CV0, furent ensuite ajustés manuellement afin d'assurer que le développement du maïs simulé lors d'une année à conditions climatiques typiques s'apparente au rendement maximal observé de même qu'au calendrier de développement établi (MAGA, s. d.) pour la région du Guatemala correspondant à notre étude de cas. La simulation des cultures débuta le 1er décembre de la première année de simulation du modèle couplé, et le maïs fut semé le 1er mai suivi d'une récolte à maturité ou au plus tard après 240 jours.

Une description du traitement des données climatiques suit ci-dessous.

### 9.3.4 Connexion des modèles

Les simulations furent implémentées avec le logiciel *Tinamït*, un paquet Python qui facilite la construction et la simulation des modèles des dynamiques des systèmes couplés à un modèle physique (J. J. Malard et al. 2017). *Tinamït* offre l'avantage d'être plus flexible et rapide à implémenter que les autres approches disponibles pour coupler des modèles des dynamiques des systèmes avec des modèles physiques externes. Tandis que les méthodes utilisant des macros en Excel (Inam et al. 2017a; 2017b) ou des extensions pour Vensim (Parsons, Nicholson, Blake, Ketterings, Ramírez-Aviles, Fox, et al. 2011; Patrick Smith, Holzworth, et Robertson 2005) furent utilisées dans le même but par le passé, ces approches sont généralement peu flexibles et nécessitent la modification du code source (et possiblement une compilation) à chaque modification auprès des variables connectées. L'approche de *Tinamït*, tout au contraire, facilite une connexion automatisée des variables des deux modèles et se charge ensuite d'échanger leurs valeurs lorsqu'une simulation est amorcée. En outre, ce logiciel permet, par l'entremise du paquet de gestion des données climatiques *تقدير* (taqdir) (J. Malard, s. d.), la connexion automatisée des variables du modèle couplé avec des bases de données climatiques externes, y compris des prévisions de changement climatique fournies par MarkSim (P. G. Jones et Thornton 2000).

Le modèle des dynamiques des systèmes fut construit avec le logiciel visuel de modélisation des dynamiques des systèmes Vensim PLE (Vensim PLE (version 7.3.5) 2017), mais toute simulation fut effectuée par l'entremise de *Tinamït* et de son moteur de simulation par défaut, PySD (Houghton 2016). Tandis que *Tinamït* offre également l'option de simuler des modèles couplés avec l'engin de simulation de Vensim, ces fonctionnalités dépendent de la bibliothèque de liens dynamiques uniquement disponible dans la version payante DSS du logiciel Vensim, bibliothèque qui, en tout cas, ne fonctionne que sur le système d'exploitation
Windows. L'expérience des auteurs est aussi que l'engin PySD s'avère généralement plus rapide que l'engin de Vensim lors d'une simulation couplée.

Une « enveloppe » spécifique à l'utilisation du modèle PCSE dans cette étude fut développée selon la spécification de Tinamït afin de permettre la connexion du modèle des cultures PCSE et du modèle des dynamiques des systèmes. Cette enveloppe gère la création du modèle PCSE pour une région donnée du Guatemala et effectue les connexions nécessaires entre Tinamït et PCSE pour le transfère des données des variables climatiques et du rendement modélisé à chaque année. Le Tableau 9.1 ci-dessous résume les variables échangées entre les différentes composantes du modèle couplé.

<table><caption>Tableau 9.1: Variables échangées entre les différentes composantes du modèle couplé. PCSE indique le modèle des cultures et MDS le modèle des dynamiques des systèmes. Tinamït/téقدir indique des variables dont les valeurs sont fournies, à travers Tinamït, par le logiciel de gestion des données climatiques téقدir [taqdir].</caption><thead><tr><td>Variable</td><td>Code</td><td>Unité</td><td>Source</td><td>Utilisation</td></tr></thead><tbody><tr><td>Rendement du maïs</td><td>Rendimiento</td><td>kg/ha</td><td>PCSE</td><td>MDS</td></tr><tr><td>Potentiel d'évaporation de l'eau</td><td>E0</td><td>cm</td><td>Tinamït/téقدir</td><td>PCSE</td></tr><tr><td>Potentiel d'évaporation du sol</td><td>ES0</td><td>cm</td><td>Tinamït/téقدir</td><td>PCSE</td></tr><tr><td>Potentiel d'évapotranspiration des cultures</td><td>ET0</td><td>cm</td><td>Tinamït/téقدir</td><td>MDS, PCSE</td></tr><tr><td>Pression de vapeur</td><td>vap</td><td>kPa</td><td>Tinamït/téقدir</td><td>PCSE</td></tr><tr><td>Précipitation</td><td>precip</td><td>mm</td><td>Tinamït/téقدir</td><td>MDS, PCSE</td></tr><tr><td>Température moyenne quotidienne</td><td>temp.prom</td><td>C</td><td>Tinamït/téقدir</td><td>PCSE</td></tr><tr><td>Température minimale quotidienne</td><td>temp.mín</td><td>C</td><td>Tinamït/téقدir</td><td>PCSE</td></tr><tr><td>Température maximale quotidienne</td><td>temp.máx</td><td>C</td><td>Tinamït/téقدir</td><td>PCSE</td></tr></tbody></table>

### 9.3.5 Données et changements climatiques

Les données climatiques historiques furent obtenues de la base de données d'observations satellitaires NASA POWER (POWER Project, s. d.), tandis que les prévisions climatiques furent obtenues du Projet de comparaison de modèles couplés (PCMC phase 5) du Programme mondiale de recherche sur le climat (PMRC) par l'entremise du logiciel MarkSim (P. G. Jones et Thornton 2000). Ces sources de données étant toutes les deux
incluses en tant que sources par défaut dans le logiciel *تقدير*, l'extraction et l'intégration de leurs données dans le modèle furent donc possibles en quelques lignes de code Python.

Un défi important lors de l'application au niveau local des prévisions climatiques des modèles globaux concerne la correction du biais intrinsèque des prévisions à grande échelle en comparaison avec les observations locales. En somme, la moyenne autant que la variabilité des prévisions générées avec un scénario sans changement climatique devraient s'apparenter à celles des observations historiques correspondantes de la région. Lors de cette étude, les prévisions de température (maximale, minimale et moyenne quotidiennes) et de précipitation furent ajustées en utilisant l'année 2010 comme date limite (inclusive) pour les données historiques.

Pour la température, la méthode de correction linéaire telle que présentée par Luo et al. (2018) fut appliquée. Lors de cette correction, la température quotidienne corrigée $T_{m,j}^{corr}$ pour le jour $j$ du mois $m$ est donnée par l'équation suivante :

$$T_{m,j}^{corr} = T_{m,j} + [ \bar{T}_{obs,m} - \bar{T}_{mod,m} ]$$

Où $T_{m,j}$ est la température prévue pour le jour $j$ du mois $m$ sous un certain régime de changement climatique, $\bar{T}_{obs,m}$ est la moyenne de la température historique observée pour le mois $m$ et $\bar{T}_{mod,m}$ la moyenne de la température du mois $m$ prévue par le modèle lorsqu'il est simulé sans changement climatique (climat historique). La même approche fut utilisée pour l'évapotranspiration de référence.

La situation se complique légèrement dans le cas de la précipitation. En soi, cette variable contient deux composantes importantes, 1) si ou non il y a de la pluie au cours une journée donnée (valeur binaire), et 2) le cas échéant, la quantité de pluie (valeur positive). Les méthodes recensées par Luo et al. (2018) ne s'appliquent que lorsque les prévisions du modèle incluent un trop grand nombre de journées pluvieuses (Vrac, Noël, et Vautard 2016), tandis que le contraire fut le cas pour cette étude. Vrac, Noël, et Vautard (2016) proposèrent donc une nouvelle méthode nommée « Élimination stochastique de singularités » où chaque journée sèche se voit assignée une valeur aléatoire entre zéro et un seuil critique représentant le valeur minimum de la pluie quotidienne recensée dans les données. Des méthodes traditionnelles de correction de biais peuvent ensuite être appliquées à l'ensemble des données, à la suite duquel les jours qui ont toujours des valeurs inférieures au seuil critique se
voient réattribuées une précipitation nulle. Cette approche fut utilisée dans cette étude, combinée avec la méthode de correction linéaire suivie d'une variante de la transformation par puissance (Luo et al. 2018).

La première correction linéaire ajuste la moyenne de précipitation mensuelle :

$$P_{m,j}^{\text{corr}} = P_{m,j} * \left[ \frac{\bar{P}_{\text{obs},m}}{\bar{P}_{\text{mod},m}} \right]$$

Où $P_{m,j}$ est la précipitation modélisée pour le jour $j$ du mois $m$, $\bar{P}_{\text{obs},m}$ la moyenne du total de la précipitation historique observée pour le mois $m$ et $\bar{P}_{\text{mod},m}$ la moyenne du total de la précipitation modélisée pour le même mois.

Suivant cette première correction, une transformation de puissance fut appliquée pour corriger la variance de la précipitation. Contrairement à la description de cette transformation par Luo et al. (2018), la correction fut appliquée dans le but de minimiser la différence entre les valeurs observées et modélisées de la précipitation annuelle cumulative à chaque mois de l'année. Cette décision fut prise suite à une analyse initiale qui démontrait qu'une simple correction au niveau de chaque mois individuel générait fréquemment un mélange de mois secs et humides au cours de la même année et ne représentait donc pas de manière satisfaisante la variabilité interannuelle en précipitation, cause clef des sécheresses et des pertes de récoltes dans la région.

Selon cette approche, une mesure $s_m$ de différence après correction entre la variabilité interannuelle de la précipitation cumulative observée et modélisée est établie pour chaque mois $m$ :

$$s_m = \sigma\left(\sum_{i=0}^{m} P_{\text{obs},i}\right) - \sigma\left[\left(P_{\text{mod},m}\right)^{b_m} * \frac{\bar{P}_{\text{mod},m}}{\left(\bar{P}_{\text{mod},m}\right)^{b_m}} + \sum_{i=0}^{m-1} P_{\text{mod},i}\right]$$

Où $\sigma$ indique la fonction d'écart type interannuel, $P_{\text{obs},i}$ est l'observation de précipitation pour le mois $i$, $P_{\text{mod},m}$ est la précipitation modélisée pour le mois $m$, $b_m$ est le facteur de correction pour le mois $m$, $\bar{P}_{\text{mod},m}$ est la moyenne interannuelle du total de la précipitation modélisée pour le mois $m$ et $P_{\text{mod},i}$ est la précipitation modélisée déjà corrigée pour le mois $i$.
L'algorithme tente ensuite d'optimiser $b_m$ afin de minimiser la valeur absolue de $s_m$. La précipitation corrigée pour chaque jour $j$ du mois $m$ est enfin calculée selon la formule suivante :

$$P_{m,j}^{\text{corr}} = P_{m,j} * \left[ \frac{(P_{\text{mod},m})^{b_m}}{P_{\text{mod},m}} \right] * \left[ \frac{\bar{P}_{\text{mod},m}}{(\bar{P}_{\text{mod},m})^{b_m}} \right]$$

Où $P_{m,j}$ est la précipitation modélisée pour le jour $j$ du mois $m$ et les autres variables sont définies tel que ci-dessus.

En commençant par le premier mois de l'année (janvier), la procédure est répétée pour chacun des 12 mois afin d'ajuster les valeurs de précipitation de chaque mois de l'année. Janvier servit de premier mois de l'année parce qu'il se situe dans une période sèche et sans culture d'awan au Tz'olöj Ya' et offre donc une bonne référence de départ pour la précipitation annuelle d'un point de vue agronomique.

Tandis que les valeurs d'évaporation et de transpiration, de même que de la pression de vapeur, sont incluses dans la base de données d'observations POWER, elles ne font pas partie des prévisions générées par MarkSim. Les valeurs d'évaporation et de transpiration pour les prévisions climatiques furent donc calculées selon l'équation Penman-Monteith telle qu'implémentée dans PCSE (de Wit 2019), tandis que la pression de vapeur fut calculée à base de la température minimum journalière selon l'implémentation en Python dans le logiciel agroécologique *Tiko'n* (Malard et al. 2020; Malard 2019) de la méthode initialement implémentée en FORTRAN par le modèle des cultures DSSAT (J. W. Jones et al. 2003).

### 9.3.6 Analyses de scénario

Six différents scénarios de politique furent appliqués en combinaison avec quatre scénarios de changement climatique (aucun changement, PCMC 2,6; PCMC 4,5; PCMC 6,0 et PCMC 8,5), pour un total de 35 différentes combinaisons, y compris le scénario de base. Les politiques, dont la majorité (incluant la reforestation, l'autosuffisance pour les intrants agricoles, les achats de groupe, les jardins familiaux et l'éducation) avaient été proposées par les parties prenantes lors des ateliers de groupe, furent choisies pour représenter une grande gamme d'actions sur les plans économique, environnemental et social.

Le Tableau 9.2 ci-dessous offre un résumé de l'implémentation de chacune de ces politiques.
<table><caption>Tableau 9.2: Spécification des politiques appliquées au modèle.</caption><thead><tr><td>Politique</td><td>Description</td></tr></thead><tbody><tr><td>Reforestation</td><td>Reforestation de 50 ha par année</td></tr><tr><td>Autosuffisance intrants agricoles</td><td>Diminution des dépenses agricoles de 75%</td></tr><tr><td>Achats de groupe d'intrants agricoles</td><td>Diminution des dépenses agricoles de 15%</td></tr><tr><td>Jardins familiaux</td><td>Cible d'autosuffisance de 90%</td></tr><tr><td>Meilleurs salaires</td><td>Salaire minimum de 20 Q de l'heure</td></tr><tr><td>Éducation</td><td>Progrès scolaire parfait</td></tr></tbody></table>

### 9.3.7 Outils informatiques

Le modèle et toutes les analyses furent implémentées en Python 3,7. Tout modèle, donnée, et code d'analyse nécessaire à la reproduction des résultats présentés au cours de cet article est disponible en tant que matériel supplémentaire.

## 9.4 Résultats et discussion

### 9.4.1 Utilité d'une approche couplée

Tel que présenté par la Figure 9.2, le modèle couplé prévoit moins de pertes des cultures en moyenne, mais distingue mieux entre les différents scénarios climatiques.

<figure><img src="image_35.png" /><figcaption>Figure 9.2: Risque de perte des cultures (lissé) selon l'année de simulation et le scénario de changement climatique (TCR = Trajectoire de concentration représentatif) avec le modèle des dynamiques des systèmes (DS) seul et couplé.</figcaption></figure>

Cette différence s'expliquerait par l'inclusion de la variabilité intra-annuelle de la précipitation dans le modèle couplé. Des déficits de précipitations encourus pendant la saison
sèche affectent le bilan d'eau total, mais n'ont pas le même impact qu'une sécheresse en période de croissance critique pour la culture du maïs. Ces impacts potentiels des changements climatiques sur le système socioéconomique n'auraient pas été pris en compte sans l'utilisation d'un modèle couplé. En outre, le couplage des modèles permet d'inclure la variabilité intra-annuelle (journalière) dans l'analyse, même si le modèle des dynamiques des systèmes lui-même n'a qu'une résolution temporelle annuelle.

La Figure 9.3 ci-dessous démontre l'influence du modèle des cultures sur la partie socioéconomique du modèle, en l'occurrence, sur la progression de la croissance de la population (gauche : utilisation du modèle des dynamiques des systèmes seul; droite : utilisation du modèle couplé). Sans le modèle des cultures, les plus grandes pertes de rendement prévues selon l'implémentation de la croissance des cultures en dynamiques des systèmes se traduisent en un taux de migration plus important et, par conséquent, une croissance populationnelle moins forte que lorsque le modèle couplé est utilisé.

<figure><img src="image_36.png" /><figcaption>Figure 9.3: Influence du modèle des cultures externe sur les projections de population avec la politique d'éducation. À gauche : modèle des dynamiques des systèmes (SD) seul ; à droite : modèle couplé. La ligne noire pointillée indique le scénario de base (sans politique ni changement climatique). TCR = Trajectoire de concentration représentatif.</figcaption></figure>

Cette approche ouvre aussi la porte à des analyses plus détaillées (intrants agricoles, différentes cultures), limitées uniquement par les capacités du modèle des cultures utilisé et les données régionales disponibles. De telles analyses, tout comme celle réalisée au cours de cette recherche, n'auraient été possibles avec l'utilisation d'un modèle des dynamiques des systèmes en isolation. De l'autre part, une analyse basée uniquement sur un modèle des
cultures manquerait les composantes socioéconomiques du système et les réponses complexes de la société face aux changements climatiques et aux autres défis.

### 9.4.2 Changements climatiques et résilience

La Figure 9.4 ci-dessous démontre le changement auprès de plusieurs variables d'intérêt environnemental, économique ou social occasionné par de différentes combinaisons de politiques et de scénarios de changement climatique.

<figure><img src="image_37.png" /></figure>

**Figure 9.4:** Impact relatif de chaque combinaison de politique et scénario de changement climatique sur différentes variables clefs au cours des 70 dernières années de la simulation. Abréviations de politiques : B=Base, Rf=Reforestation, AS=Autosuffisance produits agricoles, AG=Achats de groupe, JF=Jardins familiaux, MS=Meilleurs salaires, Éd=Éducation. Les chiffres (0; 2,6; 4,5; 6,0 et 8,5) indiquent le scénario climatique TCR (Trajectoire de concentration représentatif). Toute valeur est relative à la valeur de la même variable pour le scénario B0.

Seuls les politiques de salaire minimum ou d'éducation universelle mènent à d'importants changements dans les variables d'intérêt. Elles sont même plus efficaces que la politique de reforestation envers l'augmentation de la couverture forestière. Une analyse temporelle (voir Figure 9.5 ci-dessous) démontre que tandis que ces deux politiques s'avèrent plus efficaces au long terme, la reforestation offre des résultats plus rapides. Ces résultats suggèrent donc la possibilité de combiner des politiques pour maximiser les impacts et au court terme, et au
futur plus lointain. Il est toutefois important de préciser que les politiques d'éducation et de
salaire minimum entraînent tous les deux une forte migration hors de la région d'étude, plus
forte même que l'impact graduel que des changements climatiques de plus en plus sévères
occasionnent sur les taux de migration. Figure 9.6 ci-dessous éclaircit les processus et
rétroactions responsables pour ces comportements. Lorsque l'éducation s'améliore, une plus
haute demande pour des emplois de qualité mène à une réduction de la pauvreté (cycle positif
auto-renforçant), et, si la demande est trop importante pour l'économie locale, à une exode
des villages vers les centres urbains (cycle d'équilibre).

<figure><img src="image_38.png" /><figcaption>Figure 9.5: Effet des politiques de reforestation (gauche) et de salaire minimum (droite) sur la couverture forestière. TCR : Trajectoire de concentration représentatif.</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>Figure 9.6: Diagramme simplifié des dynamiques d'emplois et de salaire dans le modèle. R indique un cycle de rétroaction renforçant, tandis que É indique un cycle d'équilibre.</figcaption></figure>

Ces résultats surlignent aussi l'importance, dans une telle analyse socio-environnementale, d'inclure une large gamme de variables d'intérêt, car les mêmes changements auprès des mêmes variables peuvent être plus ou moins désirables selon le contexte. Une diminution dans la population, par exemple, peut être le produit d'une politique robuste d'anti-pauvreté et d'une diminution du taux de fertilité, mais peut aussi être le résultat d'une migration de masse.

De manière similaire, une diminution en pauvreté est un développement positif si elle découle d'une amélioration auprès des conditions de vie et des salaires des résidents, mais pourrait même être moins positive si elle reflète une exode populationnel suivant une perte agricole, laissant que les résidents les mieux nantis derrière. Cela tient également pour les dynamiques forestières : une augmentation de la couverture est positive si elle représente une meilleure prise en charge environnementale dans la région, mais pas si positive si elle est le résultat d'un abandon total de la région suivant une catastrophe environnementale ou autre, tel fut déjà le cas au nord du Guatemala il y a plus d'un millénaire (Evans et al. 2018; Gill et al. 2007).

Une analyse des dynamiques des systèmes encourage les parties prenantes autant que les chercheurs et chercheuses externes à réfléchir explicitement à chacun de ces processus et rétroactions parfois moins qu'intuitifs et, par conséquent, à développer des politiques plus intégrées qui prennent mieux en compte les réalités socio-environnementales complexes de la région d'étude.

Cette analyse est la première à appliquer une méthodologie de modèle des dynamiques des systèmes couplé à un modèle physique à l'analyse de la sécurité alimentaire et des systèmes
socioéconomiques et environnementaux dont elle dépend. C'est par le biais d'une telle analyse intégrée que fut démontré que les politiques appliquées au système peuvent parfois générer des impacts très forts sur son comportement, y compris sur des variables très éloignées de celles sur lesquelles elles exercent leur première influence. Une analyse qui n'incorpore pas les rétroactions entre différentes composantes du système socio-environnemental n'aurait su déceler de tels impacts. L'approche couplée permet aussi de générer des résultats plus précis que ceux qui pourraient être obtenus à l'aide d'un modèle des dynamiques des systèmes seul, tout en incorporant au sein de l'analyse les processus et rétroactions plus complexes du système.

Cette méthodologie permet de représenter les interactions complexes entre les domaines environnementaux, sociaux et économiques, interactions qui doivent être prises en compte lors du développement de politiques publiques.

## 9.5 Conclusion

La modélisation des dynamiques des systèmes est de plus en plus utilisée pour représenter les interactions entre humains et environnement. Dans sa forme participative, elle permet aussi aux parties prenantes de présenter et d'inclure leurs points de vue dans le développement du modèle et dans l'analyse des politiques. Cependant, le langage de modélisation des dynamiques des systèmes reste mal adapté à la représentation des composantes physiques des systèmes environnementaux, tels la croissance des cultures ou les changements climatiques. L'étude présente utilise un modèle des dynamiques des systèmes participatif de la sécurité alimentaire et du système agroalimentaire sous-jacent de la région majoritairement Indigène de Concepción, Tz'olöj Ya' (Guatemala). Ce modèle est couplé à un modèle externe des cultures et utilisé pour simuler l'impact sur le système alimentaire des changements climatiques et de diverses propositions de politiques pour lutter contre la malnutrition. Les résultats montrent que le modèle des dynamiques des systèmes permet de modéliser des liens complexes entre les domaines environnemental, social et économique, tandis que le modèle des cultures se distingue en apportant des prévisions de la réponse du rendement du maïs au changement climatique qui s'avérèrent plus réalistes qu'aurait été le cas avec un modèle des dynamiques des systèmes seul.
## 9.6 Remerciements

Les auteurs tiennent à remercier l'appui de l'IARNA (Universidad Rafael Landívar) lors de la recherche sur le terrain, autant que toutes les parties prenantes qui ont participé à ce projet de recherche. Cette étude fut possible grâce au financement d'une bourse de découverte du CRSNG (no. RGPIN-2015-05554) octroyée au Dr Jan Adamowski, autant que d'une Bourse de 3ᵉ cycle FRQNT, d'une Bourse du CRDI aux chercheurs candidats au doctorat (no : 107759-99906075-017), d'une Bourse d'études supérieures Canada Vanier et d'un Supplément pour études à l'étranger Michael-Smith (Programme de bourses d'études supérieures du Canada) octroyées à Julien Malard.

## 9.7 Contributions des auteurs

JM développa la méthodologie, exécuta les analyses et écrivit l'article. JA, HT et HMQ supervisèrent le projet.

## 9.8 Références

Ammar, M.E., et E.G.R. Davies. 2019. « On the accuracy of crop production and water requirement calculations: Process-based crop modeling at daily, semi-weekly, and weekly time steps for integrated assessments ». Journal of Environmental Management 238: 460-72. https://doi.org/10.1016/j.jenvman.2019.03.030.

Bojórquez-Tapia, Luis A., Marco Janssen, Hallie Eakin, Andres Baeza, Fidel Serrano-Candela, Paola Gómez-Priego, et Yosune Miquelajauregui. 2019. « Spatially-explicit simulation of two-way coupling of complex socio-environmental systems: Socio-hydrological risk and decision making in Mexico City ». Socio-Environmental Systems Modelling 1 (mai). https://doi.org/10.18174/sesmo.2019a16129.

Brisson, N, C Gary, E Justes, R Roche, B Mary, D Ripoche, D Zimmer, et al. 2003. « An Overview of the Crop Model Stics ». European Journal of Agronomy, Modelling Cropping Systems: Science, Software and Applications, 18 (3): 309-32. https://doi.org/10.1016/S1161-0301(02)00110-7.

Chapman, Alexander, et Stephen Darby. 2016. « Evaluating Sustainable Adaptation Strategies for Vulnerable Mega-Deltas Using System Dynamics Modelling: Rice Agriculture in
the Mekong Delta's An Giang Province, Vietnam ». Science of The Total Environment
559 (juillet): 326-38. https://doi.org/10.1016/j.scitotenv.2016.02.162.

Chlous-Ducharme, Frédérique, et Françoise Gourmelon. 2012. « Modélisation d'accompagnement : appropriation de la démarche par différents partenaires et conséquences: Expérience originale en milieu côtier ». VertigO, no Volume 11 Numéro 3 (février). https://doi.org/10.4000/vertigo.12163.

Di Baldassarre, Giuliano, Murugesu Sivapalan, Maria Rusca, Christophe Cudennec, Margaret Garcia, Heidi Kreibich, Megan Konar, et al. 2019. « Sociohydrology: Scientific Challenges in Addressing the Sustainable Development Goals ». Water Resources Research 55 (8): 6327-55. https://doi.org/10.1029/2018WR023901.

Diallo, Alpha Issaga Pallé, Yves Baudouin, et Michel Raymond. 2014. « Modélisation de la dynamique des systèmes socio-écologiques en vue de l'élaboration d'un cadre de référence pour la durabilité environnementale et territoriale : application au massif du Fouta-Djalon (Guinée) ». Cybergeo, mai. https://doi.org/10.4000/cybergeo.26296.

Elsawah, Sondoss, Suzanne A. Pierce, Serena H. Hamilton, Hedwig van Delden, Dagmar Haase, Amgad Elmahdi, et Anthony J. Jakeman. 2017. « An overview of the system dynamics process for integrated modelling of socio-ecological systems: Lessons on good modelling practice from five case studies ». Environmental Modelling & Software 93 (juillet): 127-45. https://doi.org/10.1016/j.envsoft.2017.03.001.

Evans, Nicholas P., Thomas K. Bauska, Fernando Gázquez-Sánchez, Mark Brenner, Jason H. Curtis, et David A. Hodell. 2018. « Quantification of Drought during the Collapse of the Classic Maya Civilization ». Science 361 (6401): 498-501. https://doi.org/10.1126/science.aas9871.

Gill, Richardson B., Paul A. Mayewski, Johan Nyberg, Gerald H. Haug, et Larry C. Peterson. 2007. « DROUGHT AND THE MAYA COLLAPSE ». Ancient Mesoamerica 18 (2): 283-302. https://doi.org/10.1017/S0956536107000193.

Halbe, Johannes, Claudia Pahl-Wostl, et Jan Adamowski. 2018. « A Methodological Framework to Support the Initiation, Design and Institutionalization of Participatory
Modeling Processes in Water Resources Management ». Journal of Hydrology 556
(janvier): 701-16. https://doi.org/10.1016/j.jhydrol.2017.09.024.

Hossain, Md Sarwar, John A Dearing, Felix Eigenbrod, et Fiifi Amoako Johnson. 2017. « Operationalizing Safe Operating Space for Regional Social-Ecological Systems ». Science of The Total Environment 584-585 (avril): 673-82. https://doi.org/10.1016/j.scitotenv.2017.01.095.

Houghton, James. 2016. PySD. Python. https://pysd.readthedocs.io/en/master/.

Inam, Azhar, Jan Adamowski, Johannes Halbe, et Shiv Prasher. 2015. « Using Causal Loop Diagrams for the Initialization of Stakeholder Engagement in Soil Salinity Management in Agricultural Watersheds in Developing Countries: A Case Study in the Rechna Doab Watershed, Pakistan ». Journal of Environmental Management 152 (avril): 251-67. https://doi.org/10.1016/j.jenvman.2015.01.052.

Inam, Azhar, Jan Adamowski, Shiv Prasher, Johannes Halbe, Julien Malard, et Raffaele Albano. 2017a. « Coupling of a Distributed Stakeholder-Built System Dynamics Socio-Economic Model with SAHYSMOD for Sustainable Soil Salinity Management – Part 1: Model Development ». Journal of Hydrology 551 (août): 596-618. https://doi.org/10.1016/j.jhydrol.2017.03.039.

———. 2017b. « Coupling of a Distributed Stakeholder-Built System Dynamics Socio-Economic Model with SAHYSMOD for Sustainable Soil Salinity Management. Part 2: Model Coupling and Application ». Journal of Hydrology 551 (août): 278-99. https://doi.org/10.1016/j.jhydrol.2017.03.040.

INE. 2019. « Resultados del Censo 2018 ». XII Censo Nacional de Población y VII de Vivienda. 2019. https://www.censopoblacion.gt/mapas.

Iwanaga, Takuya, Daniel Partington, Jenifer Ticehurst, Barry F.W. Croke, et Anthony J. Jakeman. 2020. « A Socio-Environmental Model for Exploring Sustainable Water Management Futures: Participatory and Collaborative Modelling in the Lower Campaspe Catchment ». Journal of Hydrology: Regional Studies 28 (avril): 100669. https://doi.org/10.1016/j.ejrh.2020.100669.
Jones, J. W, G Hoogenboom, C. H Porter, K. J Boote, W. D Batchelor, L. A Hunt, P. W Wilkens, U Singh, A. J Gijsman, et J. T Ritchie. 2003. « The DSSAT Cropping System Model ». European Journal of Agronomy, Modelling Cropping Systems: Science, Software and Applications, 18 (3): 235-65. https://doi.org/10.1016/S1161-0301(02)00107-7.

Jones, Peter G., et Philip K. Thornton. 2000. « MarkSim: Software to Generate Daily Weather Data for Latin America and Africa ». Agronomy Journal 92 (3): 445-53. https://doi.org/10.2134/agronj2000.923445x.

Luo, Min, Tie Liu, Fanhao Meng, Yongchao Duan, Amaury Frankl, Anming Bao, et Philippe De Maeyer. 2018. « Comparing Bias Correction Methods Used in Downscaling Precipitation and Temperature from Regional Climate Models: A Case Study from the Kaidu River Basin in Western China ». Water 10 (8): 1046. https://doi.org/10.3390/w10081046.

MAGA. s. d. « Manual de Campo para la Identificación de las Etapas de Desarrollo del Frijol en Guatemala y su Equivalencia con las Fases Fenológicas del Sistema de Monitoreo de Cultivos -SMC- ». Sistema Mesoamericano de Alerta Temprana para Seguridad Alimentaria (MFEWS). www.mfews.net.

Malard, J. J., A. Inam, E. Hassanzadeh, J. Adamowski, H. A. Tuy, et H. Melgar-Quiñonez. 2017. « Development of a software tool for rapid, reproducible, and stakeholder-friendly dynamic coupling of system dynamics and physically-based models ». Environmental Modelling & Software 96 (octobre): 410-20. https://doi.org/10.1016/j.envsoft.2017.06.053.

Malard, Julien. 2019. Tiko'n. Python.

———. s. d. تقدیر. Python. https://taqdir.readthedocs.io/fr/latest/.

Malard, Julien, Jan Adamowski, Marcela Rojas Díaz, Ling Yun (王凌云) Wang, Héctor Tuy, et Hugo Melgar-Quiñonez. 2020. « Participatory system dynamics modelling for agricultural development and on food security: small-scale agriculture and food systems in Indigenous communities of Tz'olöj Ya' and K'iche' (Guatemala) ». Agricultural Systems. En considération.
Malard, Julien, Jan Adamowski, Héctor Tuy, et Hugo Melgar-Quiñonez. 2020. « Hybrid Bayesian inference-based participatory system dynamics modelling for temporally data-poor regions: small-scale Indigenous agriculture and food systems in Tz'olöj Ya' and K'iche' (Guatemala) ». Agricultural Systems. En considération.

Malard, Julien J, Jan Franklin Adamowski, Marcela Rojas Díaz, **جسيكا بو نصار** (Jessica Bou Nassar), **نلْبَ历ْسَاةٍ مِنْ أَبَنَّةِ الْمَسَاجِد** (Nallusamy Anandaraja), Héctor Tuy, Luís Andrés Arévalo-Rodriguez, et Hugo Ramiro Melgar-Quiñonez. 2020. « Development of a New Modular Software Tool for Agroecological Food Web Model Development: Tiko'n ». European Journal of Agronomy 116 (mai): 126033. https://doi.org/10.1016/j.eja.2020.126033.

Nguyễn Thành Tựu. 2018. «agalactin-1,6-lactone production by transesterification of linoleic acid with corn oil using a new transesterification system ».[«Un modèle des dynamiques des systèmes des ressources d'eau dans le delta du Mékong face au changement climatique et aux activités humaines ».] 춘천시 한림대학교길: 한림대학교 대학원 [Chuncheon, Université Hallym].

Parsons, David, Charles F. Nicholson, Robert W. Blake, Quirine M. Ketterings, Luis Ramírez-Aviles, Jerome H. Cherney, et Danny G. Fox. 2011. « Application of a Simulation Model for Assessing Integration of Smallholder Shifting Cultivation and Sheep Production in Yucatán, Mexico ». Agricultural Systems 104 (1): 13-19. https://doi.org/10.1016/j.agsy.2010.08.006.

Parsons, David, Charles F. Nicholson, Robert W. Blake, Quirine M. Ketterings, Luis Ramírez-Aviles, Danny G. Fox, Luis O. Tedeschi, et Jerome H. Cherney. 2011. « Development and Evaluation of an Integrated Simulation Model for Assessing Smallholder Crop-Livestock Production in Yucatán, Mexico ». Agricultural Systems 104 (1): 1-12. https://doi.org/10.1016/j.agsy.2010.07.006.

Patrick Smith, F., Dean P. Holzworth, et Michael J. Robertson. 2005. « Linking Icon-Based Models to Code-Based Models: A Case Study with the Agricultural Production Systems Simulator ». Agricultural Systems 83 (2): 135-51. https://doi.org/10.1016/j.agsy.2004.03.004.
Perrone, Antonio, Azhar Inam, Raffaele Albano, Jan Adamowski, et Aurelia Sole. 2020. « A Participatory System Dynamics Modeling Approach to Facilitate Collaborative Flood Risk Management: A Case Study in the Bradano River (Italy) ». Journal of Hydrology 580 (janvier): 124354. https://doi.org/10.1016/j.jhydrol.2019.124354.

POWER Project. s. d. « POWER Data ». États-Unis d'Amérique: NASA Langley Research Center.

Srinivasan, V., M. Sanderson, M. Garcia, M. Konar, G. Blöschl, et M. Sivapalan. 2018. « Moving Socio-Hydrologic Modelling Forward: Unpacking Hidden Assumptions, Values and Model Structure by Engaging with Stakeholders: Reply to “What Is the Role of the Model in Socio-Hydrology?” » Hydrological Sciences Journal 63 (9): 1444-46. https://doi.org/10.1080/02626667.2018.1499026.

Stojkovic, Milan, et Slobodan P. Simonovic. 2019. « System Dynamics Approach for Assessing the Behaviour of the Lim Reservoir System (Serbia) under Changing Climate Conditions » . Water 11 (8): 1620. https://doi.org/10.3390/w11081620.

Vensim PLE (version 7.3.5). 2017. Ventana Systems.

Vrac, Mathieu, Thomas Noël, et Robert Vautard. 2016. « Bias Correction of Precipitation through Singularity Stochastic Removal: Because Occurrences Matter: BIAS CORRECTION OF PRECIPITATION ». Journal of Geophysical Research: Atmospheres 121 (10): 5237-58. https://doi.org/10.1002/2015JD024511.

Wit, Allard de. 2019. Python Crop Simulation Environment (version 5.4.2). Python.

———. 2020. ajwdewit/ggcmi. PLpgSQL. https://github.com/ajwdewit/ggcmi.
Connexion au chapitre 10

Le présent chapitre, qui s'inscrit dans le cadre de l'objectif 1,2 de cette thèse, applique le modèle couplé développé au cours du chapitre précédent à l'analyse des interactions entre chocs d'ordre socioéconomique et différents scénarios de changement climatique d'un point de vue de la durabilité et de la sécurité alimentaire des systèmes agricoles Indigènes à petite échelle. En comparant deux systèmes alimentaires du Guatemala, l'un axé sur la production pour l'autosuffisance et l'autre sur la production commerciale, cette approche permet de mettre en valeur les différences entre les vulnérabilités socioéconomiques et environnementales de ces deux visions alternatives du développement agricole. Elle propose aussi des recommandations résilientes aux changements climatiques et aux conditions socioéconomiques futures afin de minimiser les vulnérabilités et l'insécurité alimentaire dans chaque région.

Le présent chapitre est une traduction au français de l'article intitulé « Rukusaxïk ri k'utb'äl richin kisilonem nuk'ulem roma rutojtob'enïk runuk'ulem rik'ilal chikiwäch rujaloj ruwäch q'ij pa Tz'olöj Ya' chuqa' pa K'iche', Iximulew », présentement en révision par les pairs en vue de sa publication à titre d'article dans la revue scientifique *Revista Eutopía* (J. Malard, Adamowski, Tuy, et al. 2020c). Des coauteurs, l'auteur de la thèse développa le modèle et la méthodologie, écrivit les codes informatiques, exécuta les analyses et écrivit l'article. J. Adamowski, H. Tuy et H. Melgar-Quiñonez contribuèrent en tant que superviseurs du projet de recherche au Canada et au Guatemala.
# 10 Chapitre 10 : Application d'un modèle des dynamiques des systèmes couplé avec un modèle des cultures pour comparer les systèmes alimentaires Indigènes à petite échelle face aux changements climatiques au Tz'olöj Ya' et au K'iche', Iximulew (Guatemala)

## 10.1 Résumé

L'agriculture à petite échelle est centrale à la sécurité alimentaire de nombreuses communautés. Ces systèmes sont aussi dotés d'une très grande diversité en matière d'agriculture et de système socioéconomique, ce qui rend la question de leur comparaison et de l'analyse de leur durabilité un problème particulièrement complexe. Au cours du présent article, nous présentons l'application d'un modèle des dynamiques des systèmes couplé à un modèle des cultures à l'analyse de différentes politiques pour faire face à l'intersection des changements climatiques et des chocs socioéconomiques potentiels dans deux systèmes agricoles Indigènes à petite échelle au Guatemala, l'un axé sur l'autosuffisance et l'autre sur la production commerciale. Ce modèle, qui fut développé de manière participative avec de diverses parties prenantes des deux régions, modélise les dynamiques complexes entre composantes environnementales, économiques et sociales des systèmes agricoles. Les résultats de l'analyse suggèrent que le système d'autosuffisance est plus vulnérable aux pertes d'emplois, mais plus résilient aux fluctuations des prix agricoles; en outre, une amélioration du salaire minimum occasionnerait des améliorations marquées auprès de diverses variables clefs, y compris les indicateurs environnementaux. Ces résultats n'auraient pu être obtenus qu'avec une méthodologie participative de modélisation des dynamiques des systèmes couplée avec un modèle des cultures externe.

## 10.2 Introduction

L'insécurité alimentaire est particulièrement sévère chez les petits exploitants agricoles au Guatemala (Malard et al. 2018), autant que dans le monde entier (FAO et al. 2018). Au moment que de différents projets de développement rural sont entrepris avec l'intention d'améliorer les conditions économiques, environnementales, agricoles ou sociales, la nécessité d'une méthodologie participative d'évaluation intégrative et sensible à la complexité de ces systèmes devient évidente.
Ces projets, qui poursuivent fréquemment des philosophies et cheminements de développement très différents, tentent tous néanmoins d'influencer les systèmes socioéconomiques et environnementaux locaux dans une direction positive. Ils bénéficieraient donc du développement d'une méthodologie intégrée et appliquée pour la comparaison de ces différentes méthodes de développement, le tout dans le but d'améliorer la résistance des projets existants, de même que de générer des recommandations pour diriger les politiques et projets futurs. Cependant, les projets de développement rural ou agricole se voient très fréquemment évalués selon les points de vue technique, économique et environnemental, ou bien selon leurs effets observés.

À titre d'exemple, plusieurs études analysèrent des propositions de développement de projets communautaires, tels que l'hôtellerie (Barrera Martínez 2012), la capture des eaux pluviales (Barillas Quezada 2018), la gestion agricole (Valenzuela Pos 2013) et l'adoption de semences hybrides (Valdivia Bernal et al. 2007), dans différentes régions du Guatemala selon leur faisabilité technique et potentiel économique. Sur le plan économique, d'autres études en Amérique se sont penchées sur l'estimation du rendement économique de divers projets, y compris le potentiel d'un plan de développement économique intégré pour une municipalité au Guatemala (Chubay Gallina et al. 2010), autant que de l'agriculture urbaine en Colombie (Zapata et al. 2020). Mentionnons que plusieurs de ceux-ci (Barrera Martínez 2012; Valenzuela Pos 2013; Zapata et al. 2020) incluent aussi dans leur procédure d'évaluation la considération des impacts potentiels sur l'environnement.

D'autres tentent également d'évaluer les impacts découlant de l'implémentation d'un projet, entre autres l'utilisation de la méthodologie participative des « histoires-photos » pour analyser l'impact d'un projet de renforcement des capacités communautaires dans la communauté Tseltal (Maya) mexicaine (Gomez et al. 2017) et l'évaluation d'ateliers communautaires au sujet de l'alimentation dans le cadre d'un projet d'intervention sociale (Amaya-Castellanos et al. 2019), également au Mexique.

Hors du contexte de l'évaluation des projets de développement, les systèmes agricoles se voient souvent évalués selon des indexes de performance soit environnementaux, agronomiques ou économiques, telle l'étude de Moeller et al. (2014) dans laquelle fut analysée, à l'aide d'une série d'objectifs et d'indicateurs, la durabilité des systèmes de production de blé, et le système d'analyse de durabilité DEXiPM (Pelzer et al. 2012), qui inclut 75 indicateurs représentatifs des composantes environnementales, sociales et économiques de la durabilité des systèmes agricoles. D'autres recherches ont aussi fait appel
à des modèles économiques ou à des régressions statistiques (Malard et al. 2018; Matuschke, Mishra, et Qaim 2007) pour analyser les impacts ou les facteurs explicatifs de phénomènes telles la sécurité alimentaire ou l'adoption des technologies par les communautés agricoles.

Ces méthodologies, malgré fort utiles, demeurent linéaires et ne prennent pas compte des rétroactions potentielles entre les composantes environnementales, sociales et économiques du système. Elles risquent donc de manquer de vision intégrée des dynamiques responsables pour le comportement du système que ces interventions visent justement à modifier. En outre, des relations non linéaires entre les sphères humaines et environnementales compliquent les réponses de ces systèmes aux interventions et peuvent même mener à des réponses contre-intuitives (Stephens et al. 2012; Hossain et al. 2020; Inam et al. 2017b). Tel que formulé par Packett et al. (2020), le processus d'évaluation ne peut être libre de biais inconscients de la part des chercheuses et chercheurs, et ces biais peuvent influencer de manière importante les résultats et les recommandations du projet. La participation des parties prenantes est donc incontournable si l'on désire obtenir des résultats représentatifs de la réalité autant que des recommandations appropriées pour la région d'étude.

Dans ce contexte, la modélisation participative, et en tout particulier celle axée sur la
modélisation des dynamiques des systèmes et qui est par ce fait capable de s'adresser aux
enjeux d'intégration énumérés ci-dessus, jouit d'une popularité croissante dans le cadre des
analyses socio-environnementales (Voinov et Bousquet 2010; Turner et al. 2016; Ahmadi et
Zarghami 2019; Adamowski et Halbe 2011; Halbe et Adamowski 2019), y compris l'analyse
participative des systèmes agricoles. Entre autres, Hossain et al. (2017; 2020) évaluèrent les
risques des changements climatiques sur les communautés côtières du Bangladesh à l'aide
d'un modèle participatif des dynamiques des systèmes, tandis qu'Inam et al. (2017b; 2017a;
2015) développèrent un modèle (également participatif) de la gestion de l'eau dans le
contexte agricole de la salinisation des sols au پنجاب (Penjab) au Pakistan. La modélisation
des dynamiques des systèmes fut également appliquée à l'analyse des systèmes
agropastoraux au Sahel (Rasmussen et al. 2012) et à la sécurité alimentaire et l'agriculture à
petite échelle au Kenya (Stephens et al. 2012). En Amérique, d'autres chercheuses et
chercheurs (Parsons, Nicholson, Blake, Ketterings, Ramírez-Aviles, Fox, et al. 2011; Parsons,
Nicholson, Blake, Ketterings, Ramírez-Aviles, Cherney, et al. 2011) développèrent un modèle
des dynamiques des systèmes, couplé à un modèle externe des cultures, afin de représenter
les systèmes agricoles à petite échelle au Yucatán (Mexique).
L'importance d'une vraie participation lors d'activités de modélisation environnementale visant à générer des recommandations et des politiques publiques fut récemment soulignée par Villamor et al. (2019) de même que par Packett et al. (2020), lesquels illustrent le fait que la participation est un outil important pour combiner les points de vue qui découlent de différentes connaissances et expériences de vie afin de diminuer les biais personnels des scientifiques (autant que des parties prenantes elles-mêmes) qui risqueraient autrement de se glisser dans la conceptualisation du modèle.

Dans le cas spécifique du Guatemala, Noj Costop (2014) effectua un projet de cartographie participative pour favoriser un aménagement du territoire plus inclusif au département de San Marcos. Une étude précédente par les auteurs du présent article (J. Malard, Adamowski, Tuy, et al. 2020b; J. Malard, Adamowski, Rojas Díaz, et al. 2020; J. Malard, Adamowski, Tuy, et al. 2020a) utilisa tant qu’à elle une méthodologie de modélisation participative pour développer un modèle des dynamiques des systèmes du système agroalimentaire dans la région du Tz’olöj Ya’. Ces études facilitèrent l’identification des cycles et des dynamiques clefs dans la région et des vulnérabilités et opportunités pour des politiques visant à protéger contre les impacts des changements climatiques sur l’agriculture et les systèmes socioéconomiques connexes.

Ces dernières études, toutefois, n’analysent pas les combinaisons de chocs climatiques et socioéconomiques. Ces combinaisons peuvent pourtant s’avérer importantes à la détermination de la résilience du système et de ses réponses aux changements climatiques (Hossain et al. 2017), car si c’est les changements climatiques qui stressent le système, ce sont ses composantes socioéconomiques qui détermineront la sévérité de l’impact de ces stresses et les opportunités pour l’adaptation.

Lors de cette étude, nous appliquons une méthodologie de construction de modèle des dynamiques des systèmes couplé à un modèle externe des cultures pour analyser la performance et la résilience au long terme, face aux changements climatiques et à de divers chocs socioéconomiques, de différents programmes de développement agricole, avec comme étude de cas deux régions au Guatemala. L’une, Concepción au Tz’olöj Ya’, suit une modèle de développement axé sur la production agricole commerciale, tandis que l’autre, Chiche’ au K’iche’, suit un système de développement axé de façon prioritaire sur l’autosuffisance alimentaire.
Cette recherche propose comme objectifs 1) d'utiliser une méthodologie participative de
modélisation des dynamiques des systèmes afin d'analyser la résilience de différents
systèmes agricoles Indigènes à petite échelle et de relier celle-ci à ses causes structurelles et
2) de proposer des politiques intégrées spécifiques à chaque région visant à améliorer la
résilience de chacune. Ce projet est le premier à utiliser une méthodologie des dynamiques
des systèmes participative pour comparer des projets de développement agricole rural suivant
de différentes philosophies de développement.

## 10.3 Méthodologie

Cette recherche utilise le modèle développé lors d'une étude précédente (J. Malard, Adamowski, Tuy, et al. 2020b; J. Malard, Adamowski, Rojas Díaz, et al. 2020; J. Malard, Adamowski, Tuy, et al. 2020a) par les auteurs de la présente étude. Le modèle fut développé selon la méthodologie de modélisation participative des dynamiques des systèmes de Halbe, Pahl-Wostl, et Adamowski (2018) avec la participation de diverses parties prenantes provenant de chaque étude de cas (voir la section 10.3.1). Lors d'entrevues individuelles, des diagrammes de boucles causales furent élaborés par chaque participante et participant du processus de modélisation de groupe. Ces diagrammes qualitatifs furent par la suite combinés en un modèle qualitatif représentant chaque étude de cas (J. Malard, Adamowski, Rojas Díaz, et al. 2020), puis finalement combinés en un seul modèle quantitatif (J. Malard, Adamowski, Tuy, et al. 2020b) qui fut présenté au cours de plusieurs ateliers de modélisation de groupe dans chaque région d'étude.

Ce modèle des dynamiques des systèmes fut ensuite couplé (J. Malard, Adamowski, Tuy, et al. 2020b) avec un modèle externe des cultures, PCSE (de Wit 2019; Diepen et al. 1989), afin de modéliser l'impact des changements climatiques sur l'agriculture et la sécurité alimentaire au Tz'olöj Ya', et c'est cette version, adaptée pour être également applicable au K'iche', qui fut utilisée lors de la présente étude. Cette section présente brièvement le modèle et son fonctionnement; voir (J. Malard, Adamowski, Rojas Díaz, et al. 2020; J. Malard, Adamowski, Tuy, et al. 2020b; 2020a) pour une description complète de son développement.

Les simulations et analyses présentées dans cette étude furent effectuées en Python 3,7 avec le logiciel *Tinamït* (J. J. Malard et al. 2017). Tout code, donnée et modèle nécessaire pour la reproduction de ces résultats est inclus en tant que matériel supplémentaire.
### 10.3.1 Études de cas

Cette étude fut effectuée dans deux différentes municipalités du Guatemala, Concepción (Tz'olöj Ya') et Chiche' (K'iche'). L'agriculture à petite échelle est centrale aux économies des deux régions, qui souffrent toutefois de très hauts niveaux d'insécurité alimentaire et de malnutrition (au-delà de 50%). Cependant, chacune de ces communautés suit un philosophie de développement très différente, la première (Chiche') optant pour le développement axé sur la production pour l'autosuffisance, la seconde (Concepción) optant pour l'agriculture commerciale. Ces décisions entraînent des différences marquées dans les cultures prédominantes (maïs, fèves et autres cultures du système d'awan dans le système d'autosuffisance; légumes et fruits dans le système commercial), autant que dans les conditions socioéconomiques qui encadrent chaque système (J. Malard, Adamowski, Rojas Díaz, et al. 2020).

### 10.3.2 Modèle des dynamiques des systèmes

Le modèle des dynamiques des systèmes, portant sur les aspects socioéconomiques de la région, fut développé tel que décrit lors d'études précédentes (J. Malard, Adamowski, Tuy, et al. 2020b; J. Malard, Adamowski, Rojas Díaz, et al. 2020). Il contient des modules économique, d'utilisation du territoire, d'agriculture (commerciale et de subsistance), de population et d'alimentation et de nutrition et fut construit par le biais d'un processus participatif avec plusieurs parties prenantes de chaque région d'étude. Voir la Figure 10.1 pour un diagramme explicatif de la méthodologie participative.

<figure><img src="image_40.png" /><figcaption>Figure 10.1: Procédure de modélisation participative. De J. Malard, Adamowski, Rojas Díaz, et al. (2020).</figcaption></figure>
Tout en commençant par la modélisation qualitative par l'entremise d'entrevues individuelles (Figure 10.1, étapes 1 à 3) (J. Malard, Adamowski, Rojas Díaz, et al. 2020), le projet procéda ensuite à la quantification du modèle et à son développement continu au cours de plusieurs ateliers de modélisation de groupe (Figure 10.1, étape 5) (J. Malard, Adamowski, Tuy, et al. 2020b; 2020a).

Notons toutefois que le modèle conceptuel, ayant été développé avec des parties prenantes provenant des deux études de cas, utilise une structure identique dans les deux cas et ne différencie entre les régions que par l'entremise des valeurs initiales (p. ex., populations, éducation, utilisation du territoire) et de certaines constantes (p. ex., rendement économique de l'agriculture et des petites entreprises) qui furent estimées selon les données disponibles du sondage national ENCOVI 2011.

Toute autre différence entre les régions découle donc directement de ces différences initiales. À titre d'exemple, le fait que Chiche' se dédie beaucoup moins à l'agriculture commerciale que Concepción n'est pas une propriété prédéterminée du modèle tel qu'appliqué à chaque région, mais plutôt le résultat des différences dans les opportunités économiques de l'agriculture commerciale, lesquelles se traduisent ensuite en incitatifs pour les agricultrices et agriculteurs dans le modèle. Cette approche permet donc de modéliser des changements potentiels dans les régimes socio-environnementaux présents dans chaque région en réponse aux changements climatiques ou socioéconomiques.

### 10.3.3 Modèle des cultures

Tel qu'établi auparavant (Inam et al. 2017a; J. J. Malard et al. 2017; Ammar et Davies 2019), la modélisation des dynamiques des systèmes n'est pas idéale pour la représentation des systèmes environnementaux ou physiques complexes, et il existe déjà des modèles physiques bien adaptés à la question. C'est pour cette raison que J. Malard, Adamowski, Tuy, et al. (2020a) développèrent un modèle couplé, c.-à-d., un modèle externe des cultures (PCSE (de Wit 2019)) qui fut connecté au modèle des dynamiques des systèmes par l'entremise du logiciel Tinamït (J. J. Malard et al. 2017). Ce même modèle des cultures fut calibré pour la production de maïs selon les informations disponibles au sujet de son développement dans la région d'intérêt (MAGA, s.d.) suivant une méthodologie déjà appliqué à la région du Tz'olöj Ya' (J. Malard, Adamowski, Tuy, et al. 2020a). Lors de la présente étude, le modèle fut cependant calibré de manière indépendante pour chacune des deux régions d'étude, celles-ci, malgré leurs similitudes, ayant des climats et des semences distinctes.
### 10.3.4 Changements climatiques

Les changements climatiques furent modélisés tel que proposé par J. Malard, Adamowski, Tuy, et al. (2020a). Afin de corriger le biais des prévisions climatiques envers les valeurs observées, la méthode de correction linéaire de Luo et al. (2018) fut appliquée aux données quotidiennes de température, tandis qu'une modification de la méthode nommée « Élimination stochastique de singularités » de Vrac, Noël, et Vautard (2016) fut combinée avec une transformation de puissance (Luo et al. 2018; J. Malard, Adamowski, Tuy, et al. 2020a) pour corriger de manière simultanée les biais auprès de la fréquence et de l'intensité de la précipitation. Le biais fut corrigé de manière indépendante pour chacune des deux régions de l'étude.

Les données observées furent obtenues de la base de données d'observations satellitaires NASA POWER (POWER Project, s.d.), tandis que chaque scénario climatique fut généré par l'entremise du modèle MarkSim (Jones et Thornton 2000) selon les prévisions du Projet de comparaison de modèles couplés (PCMC phase 5) orchestré par le Programme mondiale de recherche sur le climat (PMRC).

### 10.3.5 Choix de politiques

Afin faciliter la comparaison des résultats, les mêmes politiques que celles appliquées au cours d'une étude précédente (J. Malard, Adamowski, Tuy, et al. 2020a) furent utilisées dans cette étude. Ces politiques visent à améliorer la situation nutritionnelle ou agricole par de différentes voies, soit en améliorant la situation agricole, soit en offrant de meilleures opportunités économiques à la population. Voir le Tableau 10.1 ci-dessous pour une description de chaque politique considérée. La majorité des propositions de politique (reforestation, autosuffisance des produits agricoles, achats de groupe, jardins familiaux, et éducation) découlent de propositions formulées par les parties prenantes au cours des ateliers de modélisation participative.

<table><caption>Tableau 10.1: Politiques appliquées au modèle.</caption><thead><tr><td>Politique</td><td>Description</td><td>Implémentation dans le modèle</td></tr></thead><tbody><tr><td>Reforestation</td><td>Programme de reforestation au long terme</td><td>Reforestación humana = 50 ha/année</td></tr><tr><td>Autosuffisance produits agricoles</td><td>Production locale, à moindre coût, des intrants agricoles</td><td>Gastos agrícolas = gastos agrícolas * 0,75</td></tr><tr><td>Achats de</td><td>Moindre coût des intrants agricoles</td><td>Gastos agrícolas = gastos agrícolas *</td></tr></tbody></table>
<table><thead><tr><td>Politique</td><td>Description</td><td>Implémentation dans le modèle</td></tr></thead><tbody><tr><td>groupe</td><td>par l'entremise de l'achat en gros</td><td>0,85</td></tr><tr><td>Jardins familiaux</td><td>Utilisation des jardins familiaux pour subvenir aux besoins alimentaires des familles</td><td>Meta autosuficiencia = 0,9</td></tr><tr><td>Meilleurs salaires</td><td>Augmentation réelle du salaire minimum</td><td>Salario mínimo = 20 Q/h</td></tr><tr><td>Éducation</td><td>Éducation universelle</td><td>Progreso escolar = 1</td></tr></tbody></table>

### 10.3.6 Scénarios climatiques et chocs socioéconomiques

Une combinaison de 5 scénarios climatiques (0, 2,6; 4,5; 6,0; et 8,5; chaque chiffre correspondant à un scénario du PCMC phase 5 et 0 indiquant l'absence de changement climatique) et de 5 chocs socioéconomiques fut appliquée à chacune des 6 politiques proposées, pour un total de 210 simulations pour chacune des deux municipalités. Cette analyse plus approfondie permet de mesurer les impacts, auprès de la durabilité des interventions potentielles, des interactions entre des conditions hors du contrôle des parties prenantes et les scénarios climatiques futurs potentiels. Le modèle fut simulé pour une période de 105 ans, de 1988 à 2093.

Le Tableau 10.2 ci-dessous décrit les chocs socioéconomiques appliqués au modèle.

<table><caption>Tableau 10.2: Chocs socioéconomiques appliqués au modèle.</caption><thead><tr><td>Choc</td><td>Description</td><td>Implémentation dans le modèle</td></tr></thead><tbody><tr><td>Agrocommercial</td><td>Perte de marchés pour la production agricole commerciale</td><td>Rendimiento agrícola comercial = 0</td></tr><tr><td>Économique</td><td>Crise économique engendrant une diminution systématique des salaires</td><td>máx sueldo = 20 Q/h</td></tr><tr><td>Emplois</td><td>Crise économique provoquant une perte de la moitié de la capacité d'emploi</td><td>Disponibilidad empleos = 0,5</td></tr><tr><td>Inflation</td><td>Inflation des prix des aliments</td><td>Inflación precios = 1,5</td></tr><tr><td>Agrocommercial et économique</td><td>Combinaison des chocs agrocommercial et économique</td><td>Rendimiento agrícola comercial = 0<br>Disponibilidad empleos = 0,5</td></tr></tbody></table>

Ces chocs furent choisis en raison de leur importance socioéconomique aux systèmes agricoles locaux. Chacun affecte une variable à la frontière du modèle, là ou le modèle au
niveau municipal s’arrête et où il doit par conséquent se reposer sur des suppositions, implicites ou explicites, de l’état du système externe. En appliquant des chocs à ces niveaux-ci, l’analyse permet d’évaluer la résilience du système modélisé à des conditions externes que le modèle lui-même ne peut pas prédire.

En titre d'exemple, trois des chocs (économique, agrocommercial et la combinaison des deux) traitent des changements défavorables qui pourraient survenir aux marchés agricole et du travail. Le choc d'emplois, tant qu'à lui, simule une diminution dramatique de la capacité de la région de soutenir des emplois salariés, tel que pourrait être le cas suivant une crise économique d'ordre nationale ou supranationale. Finalement, le choc d'inflation simule une inflation du coût de la nourriture, éventualité probable dans le cas d'une perte des cultures généralisée dans le pays ou même en raison d'une crise économique internationale, tel que fut le cas au niveau mondial en 2008.

Cette approche facilite, de même, l'analyse comparative des deux systèmes agroalimentaires, systèmes qui sont, rappelons-le, représentés par le même modèle informatique qui ne diffère que selon leurs valeurs initiales et les définitions de certaines constantes externes (tel le rendement économique de l'agriculture commerciale) dans chaque région.

## 10.4 Résultats et discussion

Les Figures 10.4 et 10.5 à la fin de cette section résument de manière visuelle les résultats des analyses de vulnérabilité et de résilience. Dans les deux localités, les chocs économique, d'inflation, et des emplois s'avèrent les plus néfastes pour la malnutrition. Ceci peut s'expliquer par la boucle de rétroaction de renforcement entre l'éducation, la pauvreté, et la nutrition (Figure 10.2 ci-dessous). Tandis qu'une meilleure éducation peut améliorer beaucoup d'autres variables connexes et pousser le système entier vers un état stable plus favorable, ces trois chocs ont tous le potentiel d'occasionner l'impact opposé sur le système, et d'empêcher la région de sortir de la pauvreté.
<figure><img src="image_41.png" /><figcaption>Figure 10.2: Diagramme de boucles causales représentant la relation entre la malnutrition et l'économie.</figcaption></figure>

Au contraire, les politiques du salaire minimum ou de l'éducation universelle peuvent toutes les deux mitiger la malnutrition lorsqu'aucun autre choc est présent. Cependant, de ces deux, seul le salaire minimum est résilient face aux chocs d'ordre économique (incluant l'inflation et la perte d'emploi). Tel que démontré dans la Figure 10.2 ci-dessus, une intervention d'éducation reste toutefois vulnérable aux troubles d'ordre économique, tandis qu'un salaire minimum adresse le problème plus directement. L'éducation présente également un certain délai dans son mode d'action (plusieurs années doivent s'écouler entre l'éducation des jeunes d'aujourd'hui et leur entrée dans le marché du travail), tandis qu'une intervention de salaire minimum peut prendre effet immédiatement.

Les politiques environnementales et agricoles, pour leur part, s'avèrent en grande majorité vulnérables aux changements climatiques combinés avec des chocs économiques dans la région de Concepción. Lorsque comparée à la même politique appliquée à un scénario sans changement climatique, la politique de l'autosuffisance alimentaire par l'entremise des jardins familiaux se montre même plus vulnérable que le scénario sans intervention. Il semble qu'une plus grande dépendance sur l'agriculture en temps de récoltes incertaines risque effectivement d'augmenter la vulnérabilité du système.

La migration, de son côté, est vulnérable aux chocs agrocommercial et d'emplois, ces deux-ci augmentant la quantité de migration de manière dramatique. Les jardins familiaux, tant qu'à eux, augmentent la résistance de la migration face à un choc agrocommercial, vraisemblablement dû à une importance plus importante placée sur l'autosuffisance au dépens de la production commerciale. En outre, Chiche', où l'agriculture commerciale a moindre importance, se voit complètement à l'abri des impacts d'un choc des prix agrocommerciaux, tel que pourrait être prévu (voir Figure 10.5). Cependant, une moindre dépendance sur les systèmes commerciaux augmente aussi le risque de migration face à une perte d'emplois, vu
que les alternatives économiques dans le secteur agricole se voient limitées (Figure 10.3 ci-dessous).

<figure><img src="image_42.png" /><figcaption>Figure 10.3: Diagramme de boucles causales représentant les dynamiques économiques du système agroalimentaire.</figcaption></figure>

Il demeure intéressant que, à Concepción, des meilleurs salaires ou de la meilleure
éducation diminuent la migration dans le cas d'un choc d'emplois ou d'un choc
agrocommercial, mais l'augmente dans les autres cas. Cependant, ces migrations ne sont pas
une « migration forcée » par les changements climatiques ou autre crise, mais plutôt le
résultat de personnes ayant le potentiel de gagner des meilleurs revenus qui quittent les
régions rurales à la recherche d'emplois. Pour ce qui est de la pauvreté, un salaire minimum
mêne à un réduction générale, tandis que l'éducation génère un impact similaire dans la
majorité des cas mais risque d'empirer la pauvreté dans le cas d'un choc économique à
Concepción, justement parce que l'économie a transitionnée vers une système moins agricole
et plus dépendant du système économique.

Dans le case des changements climatiques, il peut être surprenant de constater que des
changements plus sévères mènent à moins de pauvreté. Cependant, cette diminution de la
pauvreté telle que mesurée n'est pas nécessairement un développement positif : au contraire,
elle peut être l'effet de pertes agricoles qui forcent les personnes à adopter un emploi, tout en
étant accompagnée d'une augmentation en malnutrition en raison de la perte de capacité de
production de la nourriture locale. Ces résultats soulignent l'importance d'évaluer plusieurs
indicateurs dans ces analyses de modélisation socioéconomique, et de toujours rester vigilant envers des résultats potentiellement trompeurs ou contreproductifs à première vue.

Du côté environnemental, la reforestation offre certainement un impact positif envers la couverture forestière, mais un salaire minimum donne des effets plus robustes envers les chocs socioéconomiques. Malgré le fait que les changements climatiques tendent à diminuer l'impact d'un salaire minimum sur la couverture forestière, notons que les impacts de celui-ci demeurent néanmoins supérieurs à ceux de la reforestation, même sous le scénario 8,5 de changement climatique maximal. De son côté, l'éducation est efficace comme politique mais est très vulnérable envers un choc économique, lorsque la population se replie vers un système économique axé sur l'agriculture.

À Chiche', là où l'agriculture commerciale est moindre, l'impact d'une intervention de salaire minimum sur toutes les variables d'intérêt est plus marqué, et, en outre, plus résilient à de différents chocs socioéconomiques. Notons qu'à Chiche', l'impact des changements climatiques était nettement mois néfaste qu'à Tz'olöj Ya'. Ceci ne correspond cependant pas aux observations à date de la région ("La sequía afectó a 50,000 familias y dañó 42,000 hectáreas en Guatemala | El Economista" s. d.; Felipe Q 2016). Plusieurs facteurs pourraient être en cause, entre eux 1) l'utilisation des données satellitaires comme point de référence historique pour l'ajustement du biais dans les données climatiques, 2) la basse résolution spatiale des ces données satellitaires, 3) la quantité limité de données satellitaires historiques et 4) la basse résolution des modèles de changement climatique globaux utilisés pour les prévisions du futur. Lors d'études suivantes, ces limitations pourraient être mitigées si des observations locales (sur le terrain) d'une longue durée historique pouvaient être obtenues, ou si un modèle de circulation régional était développé pour le Guatemala afin de mieux corriger les biais dans les modèles climatiques (Ibrahim, Karambiri, et Jan 2011). Entre temps, nous recommandons la prudence et la réflexion critique lors de l'analyse des résultats de prévision des changements climatiques.

En bref, un accès à de meilleurs salaires paraît être la solution la plus robuste face aux changements climatiques et à de différents chocs socioéconomiques afin d'améliorer la situation nutritionnelle et environnementale, et ce, dans le cas du système d'autosuffisance de même celui de commercialisation agricole. Ce dernier, cependant, se montre plus résilient face aux pertes d'emplois, au prix, toutefois, d'une plus grande vulnérabilité aux pertes de revenu agricole.
Un accès universel à l'éducation s'avère également efficace envers les enjeux environnementaux et humains, mais demeure moins résilient que l'instauration d'un salaire minimum face aux chocs socioéconomiques. Cependant, ces politiques risquent d'augmenter la migration hors des zones rurales. Tandis que cela n'indique pas nécessairement un impact négatif, ces éventualités doivent être prises en considération lors du développement de politiques intégrées pour le développement rural.

Ces résultats, qui combinent impacts environnementaux, économiques et sociaux autant que les rétroactions parfois inattendues qui les relient, n'auraient pu être décelés qu'avec une méthodologie participative de modélisation des dynamiques des systèmes couplée avec un modèle physique externe des cultures telle que fut élaborée lors de cette étude.
<figure><img src="image_43.png" /><figcaption>Figure 10.4: Interactions entre politiques (horizontal) et chocs socioéconomiques (vertical) pour différentes variables de Concepción (Tz'olöj Ya'). Première rangée : Valeurs de diverses variables au cours des 70 dernières années de la simulation sous le scénario climatique 8,5, relatives au scénario de base sans changement climatique. Deuxième rangée : Impact du scénario climatique 8,5 pour chaque combinaison de politique et de choc, relatif au résultat de la même combinaison sans changement climatique (changement relatif; 0 indique qu'il n'y a aucun changement en comparaison avec le scénario sans changement climatique). Abréviations de politique : B = scénario de base (aucune politique); Rf = reforestation; AS = autosuffisance produits agricoles; AG = achats de groupe; JF = jardins familiaux; MS = meilleurs salaires; Éd = éducation. Abréviations de choc : Agr = agrocommercial; Écon = économique; Empl = emplois; Infl = inflation; AgrÉcon = agrocommercial et économique.</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>Figure 10.5: Interactions entre politiques (horizontal) et chocs socioéconomiques (vertical) pour différentes variables de Chiche' (K'iche'). Première rangée : Valeurs de diverses variables au cours des 70 dernières années de la simulation sous le scénario climatique 8,5, relatives au scénario de base sans changement climatique. Deuxième rangée : Impact du scénario climatique 8,5 pour chaque combinaison de politique et de choc, relatif au résultat de la même combinaison sans changement climatique (changement relatif; 0 indique qu'il n'y a aucun changement en comparaison avec le scénario sans changement climatique). Abréviations de politique : B = scénario de base (aucune politique); Rf = reforestation; AS = autosuffisance produits agricoles; AG = achats de groupe; JF = jardins familiaux; MS = meilleurs salaires; Éd = éducation. Abréviations de choc : Agr = agrocommercial; Écon = économique; Empl = emplois; Infl = inflation; AgrÉcon = agrocommercial et économique.</figcaption></figure>
## 10.5 Conclusion

Lors de cette étude, un modèle des dynamiques des systèmes couplé avec un modèle externe des cultures fut appliqué à l'analyse de la résilience de différents systèmes de développement agricole au Guatemala, l'un agrocommercial et l'autre d'autosuffisance, face aux défis simultanés des changements climatiques et de divers scénarios de chocs socioéconomiques. Les résultats démontrent des interactions fortes entre les composantes économiques, environnementales et sociales du modèle couplé et suggèrent que les mesures d'égalité économique, tel que les salaires minimums et l'éducation universelle, figurent parmi les politiques les plus résistantes aux chocs climatiques et économiques.

Toute politique, cependant, demeure un compromis entre de diverses variables (telles la nutrition, la migration et la pauvreté) qui doivent être considérées avec soin par toutes les parties prenantes locales avant l'implémentation d'un projet de développement par des parties externes.

## 10.6 Remerciements

Les auteurs tiennent à remercier l'IARNA (Universidad Rafael Landívar) pour leur appui lors de la recherche sur le terrain, autant que toutes les parties prenantes qui ont participé à ce projet de recherche. Cette étude fut possible grâce au financement d'une bourse de découverte du CRSNG (no. RGPIN-2015-05554) octroyée au Dr Jan Adamowski, autant que d'une Bourse de 3e cycle FRQNT, d'une Bourse du CRDI aux chercheurs candidats au doctorat (no : 107759-99906075-017), d'une Bourse d'études supérieures Canada Vanier et d'un Supplément pour études à l'étranger Michael-Smith (Programme de bourses d'études supérieures du Canada) octroyées à Julien Malard.

## 10.7 Contributions des auteurs

JM développa la méthodologie, exécuta les analyses et écrivit l'article. JA, HT et HMQ supervisèrent le projet.

## 10.8 Références

Adamowski, Jan, et Johannes Halbe. 2011. "Participatory water resources planning and management in an Agriculturally Intensive Watershed in Quebec, Canada using Stakeholder Built System Dynamics Models". Annals of Warsaw University of Life
Sciences - SGGW. Land Reclamation 43 (1): 3-11. https://doi.org/10.2478/v10060-008-0088-0.

Ahmadi, Mohammad Hossein, et Mahdi Zarghami. 2019. "Should Water Supply for Megacities Depend on Outside Resources? A Monte-Carlo System Dynamics Simulation for Shiraz, Iran". Sustainable Cities and Society 44 (enero): 163-70. https://doi.org/10.1016/j.scs.2018.10.007.

Amaya-Castellanos, M.A., M.D.C. Morales-Ruan, R. Uribe-Carvajal, A. Jiménez-Aguilar, A.A. Salazar-Coronel, B. Martínez-Tapia, et T. Shamah-Levy. 2019. "Implementación de un modelo de capacitación multimedial para brindar orientación alimentaria a los beneficiarios de un programa de ayuda social en México". Global Health Promotion 26 (4): 118-29. https://doi.org/10.1177/1757975917751908.

Ammar, M.E., et E.G.R. Davies. 2019. "On the accuracy of crop production and water requirement calculations: Process-based crop modeling at daily, semi-weekly, and weekly time steps for integrated assessments". Journal of Environmental Management 238: 460-72. https://doi.org/10.1016/j.jenvman.2019.03.030.

Barillas Quezada, José Andrés. 2018. "Diseño de un sistema de captación de agua lluvia por medio de aljibe para el área de la aldea Xayá, municipio de Santa Lucía Cotzumalguapa, departamento de Escuintla, Guatemala". Guatemala: Universidad de San Carlos de Guatemala. http://biblioteca.usac.edu.gt/tesis/08/08_0034_A.pdf.

Barrera Martínez, Maidelyn Cristina. 2012. "Estudio de la factibilidad para la implementación de un complejo hotelero ecológico, en el parque Senderos de Alux, para el desarrollo sostenible y crecimiento turístico de la comunidad de San Lucas Sacatepéquez". Guatemala: Universidad de San Carlos de Guatemala. http://biblioteca.usac.edu.gt/tesis/08/08_2563_IN.pdf.

Chubay Gallina, Rolando Leonel, Juan Orlando González López, Eric Benjamín López López, Edgar Rolando Sutuj Estrada, William Velásquez Gálvez, Carlos Amilcar Rodas López, Víctor Oswaldo Duque Alarcón, et Adelaida Xiloj Herrera. 2010. "Diagnóstico socioeconómico, potencialidades productivas y propuestas de inversión: Municipio de
Conceptualizing Social-Ecological System Dynamics in the Bangladesh Delta". Regional Environmental Change 20 (1): 28. https://doi.org/10.1007/s10113-020-01599-5.

Ibrahim, B., H. Karambiri, et P. Jan. 2011. "Evolution of main rainy season characteristics in Burkina Faso derived from rainfall data of five regional climate models (RCM)". En, 344:82-86.

Inam, Azhar, Jan Adamowski, Johannes Halbe, et Shiv Prasher. 2015. "Using Causal Loop Diagrams for the Initialization of Stakeholder Engagement in Soil Salinity Management in Agricultural Watersheds in Developing Countries: A Case Study in the Rechna Doab Watershed, Pakistan". Journal of Environmental Management 152 (abril): 251-67. https://doi.org/10.1016/j.jenvman.2015.01.052.

Inam, Azhar, Jan Adamowski, Shiv Prasher, Johannes Halbe, Julien Malard, et Raffaele Albano. 2017a. "Coupling of a Distributed Stakeholder-Built System Dynamics Socio-Economic Model with SAHYSMOD for Sustainable Soil Salinity Management – Part 1: Model Development". Journal of Hydrology 551 (agosto): 596-618. https://doi.org/10.1016/j.jhydrol.2017.03.039.

———. 2017b. "Coupling of a Distributed Stakeholder-Built System Dynamics Socio-Economic Model with SAHYSMOD for Sustainable Soil Salinity Management. Part 2: Model Coupling and Application". Journal of Hydrology 551 (agosto): 278-99. https://doi.org/10.1016/j.jhydrol.2017.03.040.

Jones, Peter G., et Philip K. Thornton. 2000. "MarkSim: Software to Generate Daily Weather Data for Latin America and Africa". Agronomy Journal 92 (3): 445-53. https://doi.org/10.2134/agronj2000.923445x.

"La sequía afectó a 50,000 familias et dañó 42,000 hectáreas en Guatemala | El Economista". s.d.. Consultado el 2 de junio de 2020. https://www.eleconomista.net/actualidad/La-sequia-afecto-a-50000-familias-y-dano-42000-hectareas-en-Guatemala-20190826-0042.html.

Luo, Min, Tie Liu, Fanhao Meng, Yongchao Duan, Amaury Frankl, Anming Bao, et Philippe De Maeyer. 2018. "Comparing Bias Correction Methods Used in Downscaling
Precipitation and Temperature from Regional Climate Models: A Case Study from the Kaidu River Basin in Western China". Water 10 (8): 1046. https://doi.org/10.3390/w10081046.

MAGA. s.d.. "Manual de Campo para la Identificación de las Etapas de Desarrollo del Frijol en Guatemala y su Equivalencia con las Fases Fenológicas del Sistema de Monitoreo de Cultivos -SMC-". Sistema Mesoamericano de Alerta Temprana para Seguridad Alimentaria (MFEWS). www.mfews.net.

Malard, J. J., A. Inam, E. Hassanzadeh, J. Adamowski, H. A. Tuy, et H. Melgar-Quiñonez. 2017. "Development of a software tool for rapid, reproducible, and stakeholder-friendly dynamic coupling of system dynamics and physically-based models". Environmental Modelling & Software 96 (octubre): 410-20. https://doi.org/10.1016/j.envsoft.2017.06.053.

Malard, Julien, Jan Adamowski, Marcela Rojas Díaz, Ling Yun (王凌云) Wang, Héctor Tuy, et Hugo Melgar-Quiñonez. 2020. "Participatory system dynamics modelling for agricultural development and on food security: small-scale agriculture and food systems in Indigenous communities of Tz'olöj Ya' and K'iche' (Guatemala)". Sera soumis à System Dynamics Review.

Malard, Julien, Jan Adamowski, Héctor Tuy, et Hugo Melgar-Quiñonez. 2020a. "Couplage d'un modèle des dynamiques des systèmes avec un modèle agricole : développement de politiques pour lutter contre l'insécurité alimentaire face aux changements climatiques dans un système agricole Indigène à petite échelle au Tz'olöj Ya' (Guatemala)". VertigO. En considération.

———. 2020b. "Hybrid Bayesian inference-based participatory system dynamics modelling for temporally data-poor regions: small-scale Indigenous agriculture and food systems in Tz'olöj Ya' and K'iche' (Guatemala)". Sera soumis à System Dynamics Review.

Malard, Julien J, Diana Dallmann, Ahn Bui, Lillian Yin, Miguel García-Winder, Humberto Monardes, et Hugo Melgar-Quiñonez. 2018. "Relación entre la diversidad pecuaria y la seguridad alimentaria de pequeños productores agropecuarios en Guatemala". Revista Eutopía 3 (6): 3-30.
Matuschke, Ira, Ritesh R. Mishra, et Matin Qaim. 2007. "Adoption and Impact of Hybrid Wheat in India". World Development 35 (8): 1422-35. https://doi.org/10.1016/j.worlddev.2007.04.005.

Moeller, Carina, Joachim Sauerborn, Peter de Voil, Ahmad M. Manschadi, Mustafa Pala, et Holger Meinke. 2014. "Assessing the Sustainability of Wheat-Based Cropping Systems Using Simulation Modelling: Sustainability = 42?" Sustainability Science 9 (1): 1-16. https://doi.org/10.1007/s11625-013-0228-2.

Noj Costop, Ixchebel Nacxit. 2014. "Diseño de una guía metodológica y elaboración de planes indicativos de ordenamiento territorial diagnóstico y servicios en cinco comunidades de la cuenca de los ríos Coatán y Suchiate del departamento de San Marcos, Guatemala, C.A." Guatemala: Universidad de San Carlos de Guatemala. http://biblioteca.usac.edu.gt/tesis/01/01_2969.pdf.

Packett, E., N.J. Grigg, J. Wu, S.M. Cuddy, P.J. Wallbrink, et A.J. Jakeman. 2020. "Mainstreaming gender into water management modelling processes". Environmental Modelling and Software 127. https://doi.org/10.1016/j.envsoft.2020.104683.

Parsons, David, Charles F. Nicholson, Robert W. Blake, Quirine M. Ketterings, Luis Ramírez-Aviles, Jerome H. Cherney, et Danny G. Fox. 2011. "Application of a Simulation Model for Assessing Integration of Smallholder Shifting Cultivation and Sheep Production in Yucatán, Mexico". Agricultural Systems 104 (1): 13-19. https://doi.org/10.1016/j.agsy.2010.08.006.

Parsons, David, Charles F. Nicholson, Robert W. Blake, Quirine M. Ketterings, Luis Ramírez-Aviles, Danny G. Fox, Luis O. Tedeschi, et Jerome H. Cherney. 2011. "Development and Evaluation of an Integrated Simulation Model for Assessing Smallholder Crop-Livestock Production in Yucatán, Mexico". Agricultural Systems 104 (1): 1-12. https://doi.org/10.1016/j.agsy.2010.07.006.

Pelzer, Elise, Gabriele Fortino, Christian Bockstaller, Frédérique Angevin, Claire Lamine, Camilla Moonen, Vasileios Vasileiadis, et al. 2012. "Assessing Innovative Cropping Systems with DEXiPM, a Qualitative Multi-Criteria Assessment Tool Derived from
POWER Project. s.d.. "POWER Data". États-Unis d'Amérique: NASA Langley Research Center.

Rasmussen, Laura Vang, Kjeld Rasmussen, Anette Reenberg, et Simon Proud. 2012. "A System Dynamics Approach to Land Use Changes in Agro-Pastoral Systems on the Desert Margins of Sahel". Agricultural Systems 107 (marzo): 56-64.
https://doi.org/10.1016/j.agsy.2011.12.002.

Stephens, E.C., C.F. Nicholson, D.R. Brown, D. Parsons, C.B. Barrett, J. Lehmann, D. Mbugua, S. Ngoze, A.N. Pell, et S.J. Riha. 2012. "Modeling the impact of natural resource-based poverty traps on food security in Kenya: The Crops, Livestock and Soils in Smallholder Economic Systems (CLASSES) model". Food Security 4 (3): 423-39.
https://doi.org/10.1007/s12571-012-0176-1.

Turner, Benjamin, Hector Menendez, Roger Gates, Luis Tedeschi, et Alberto Atzori. 2016. "System Dynamics Modeling for Agricultural and Natural Resource Management Issues: Review of Some Past Cases and Forecasting Future Roles". Resources 5 (4): 40.
https://doi.org/10.3390/resources5040040.

Valdivia Bernal, Roberto, Francisco de Jesús Caro Velarde, Margarito Ortiz Catón, Alberto Betancourt Vallejo, Alejandro Ortega Corona, Víctor Antonio Vidal Martínez, et Alejandro Espinosa Calderón. 2007. "Desarrollo participativo de híbridos sintéticos de maíz y producción de semilla por agricultores". Agricultura Técnica en México 33 (2): 135-45.

Valenzuela Pos, Oscar Alfredo. 2013. "Análisis de la gestión colectiva de la comunidad Indígena de Palín, con énfasis en el uso y manejo de los recursos naturales en la finca comunal 'El Chilar', Palín, Escuintla, Guatemala, C.A." Guatemala: Universidad de San Carlos de Guatemala. http://biblioteca.usac.edu.gt/tesis/01/01_2856.pdf.

Villamor, Grace B., David L. Griffith, Andrew Kliskey, et Lilian Alessa. 2019. "Contrasting stakeholder and scientist conceptual models of food-energy-water systems: a case study
in Magic Valley, Southern Idaho". Socio-Environmental Systems Modelling 2 (octubre): 16312. https://doi.org/10.18174/sesmo.2020a16312.

Voinov, Alexey, et Francois Bousquet. 2010. "Modelling with Stakeholders". Environmental Modelling & Software, Thematic Issue - Modelling with Stakeholders, 25 (11): 1268-81. https://doi.org/10.1016/j.envsoft.2010.03.007.

Vrac, Mathieu, Thomas Noël, et Robert Vautard. 2016. "Bias Correction of Precipitation through Singularity Stochastic Removal: Because Occurrences Matter: BIAS CORRECTION OF PRECIPITATION". Journal of Geophysical Research: Atmospheres 121 (10): 5237-58. https://doi.org/10.1002/2015JD024511.

Wit, Allard de. 2019. Python Crop Simulation Environment (versión 5.4.2). Python.

Zapata, V.Y.H., F.J.V. Zabala, M.A.M. Martina, A.M.S. González, et M.E. Vélez. 2020. "Sociodemographic, cultural, environmental and agroecological characterization in order to adopt urban agriculture in the municipality of Tuluá, Colombia". WSEAS Transactions on Environment and Development 16: 67-83. https://doi.org/10.37394/232015.2020.16.8.
## Connexion au chapitre 11

Les chapitres précédents de la thèse présentent des méthodes pour la modélisation participative des systèmes agroalimentaires à petite échelle, tout particulièrement, de leurs aspects socioéconomiques et des dynamiques au long terme générées par les interactions de ceux-ci avec les composantes environnementales des mêmes systèmes.

Un tout autre défi envers la modélisation des systèmes agricoles à petite échelle est la représentation des systèmes de réseaux trophiques qui déterminent les dynamiques des populations des ravageurs et des auxiliaires dans les champs. Les modèles des cultures présentement disponibles, développés prioritairement en vue d'une application à l'agriculture à grande échelle, ne prennent pas compte de ces dynamiques, d'où l'importance de développer des outils informatiques pour faciliter la modélisation des réseaux agroécologiques que de différentes communautés agricoles gèrent de manières très diverses.

La deuxième partie de cette thèse s'adresse donc au développement de méthodologies et d'outils nécessaires pour la représentation, dans un cadre de modélisation agricole, des dynamiques des réseaux trophiques agroécologiques à court et moyen terme. Les équations fonctionnelles de prédation disponibles dans la littérature ne présentant pas toutes les caractéristiques théoriques nécessaires à l'utilisation dans le cadre d'un modèle appliqué, le présent chapitre présente une nouvelle équation fonctionnelle de prédation dont la validité mathématique correspond aux nécessités des deux articles suivants, où seront exposés le développement et l'application pratique des modèles de réseaux trophiques agroécologiques. Le présent chapitre s'inscrit dans le cadre de l'objectif 2,1 de cette thèse.

Le présent chapitre est une traduction au français de l'article intitulé « Modelling predation: Theoretical criteria and empirical evaluation of functional form equations for predator-prey systems », sous presse auprès de la revue *Ecological Modelling* en tant qu'article scientifique révisé par des pairs (Julien J Malard et al. 2020). Suivra, dans le chapitre 12, la version originale en anglais. Des coauteurs de l'article, l'auteur de la thèse développa l'algorithme et les critères de sélection théoriques, effectua les calculs théoriques et empiriques, compila les données et développa le code pour la section empirique, et écrivit l'article. (Jessica Bou Nassar) contribua à la vérification des calculs théoriques, tandis que J. Adamowski, நல்லுசாமி ஆனந்தராஜா (Nallusamy Anandaraja), H. Tuy, et H.
Melgar-Quiñonez contribuèrent en tant que superviseurs de la recherche au Canada, en Inde (தமிழ்நாடு) et au Guatemala.
# 11 Chapitre 11 : Modèles de prédation : Critères théoriques et évaluation empirique des équations de forme fonctionnelle pour les systèmes de prédateurs et de proies

## 11.1 Résumé

La représentation exacte des relations entre prédateurs et proies est centrale à la modélisation écologique et à la modélisation des dynamiques des populations. Toutefois, et malgré la présence de longue date d'une compétition au sein de la littérature entre les modèles dépendants du ratio et ceux dépendants de la proie (de même que quelques formes intermédiaires), la majorité des équations qui sont présentement appliquées à la représentation de ces relations ne répondent pas aux critères théoriques nécessaires pour garantir la cohérence biologique. La présente recherche propose une liste de tels critères applicables de manière universelle à toute équation de prédation et démontre par la suite que les équations les plus populaires au sein de la littérature ne répondent pas à ces mêmes critères. Suivra une proposition d'une nouvelle équation qui répond de manière satisfaisante à ces critères et qui combine les concepts de la dépendance à la proie et de la dépendance au ratio afin de générer des prévisions raisonnables même sous des conditions de densité populationnelle élevée du prédateur ou de la proie. La performance empirique de la nouvelle équation est ensuite évaluée, de même que celles des équations alternatives existantes, selon de différents ensembles de données expérimentaux provenant de la littérature. Les résultats indiquent que la nouvelle équation est non seulement plus cohérente que les équations existantes du point de vue mathématique, mais donne des résultats empiriques plus consistants, et ce, dans le contexte de diverses situations écologiques. Cette recherche représente la première proposition d'une liste systématique de critères pour l'évaluation des équations de la prédation, suivie d'une proposition d'une nouvelle équation qui est plus performante du point de vue théorique de même qu'empirique lorsqu'appliquée à des données provenant d'une grande gamme de systèmes de prédateur et de proie.

## 11.2 Introduction

La représentation des relations entre prédateurs et proies demeure un enjeu central des études de modélisation des interactions trophiques et des dynamiques populationnelles (Kratina et al. 2009). Une grande gamme d'équations de prédation est présentement
disponible auprès de la littérature écologique; ces équations calculent les taux de prédation (proie consommée par prédateur par unité de temps) selon une variable indépendante (densité de la proie, densité du prédateur, ou les deux) et un ou plusieurs paramètres empiriques (Abrams et Ginzburg 2000). Ces équations se voient fréquemment classifiées selon la variable indépendante et la forme mathématique (appelée réponse fonctionnelle) de l'équation qui relie cette variable au taux de prédation.

Selon cette première, les équations de prédation peuvent être classifiées en tant qu'équation dépendante de la proie (la population de la proie, N, étant la seule variable indépendante) ou en tant qu'équation dépendante du ratio (la variable indépendante étant le ratio (N/P) des populations de la proie et du prédateur, P). Notons également une grande gamme de formes intermédiaires, dont figurent parmi les plus communes les formes Hassell-Varley (Hassell et Varley 1969) et l'équation Beddington-DeAngelis (αN/((N+cP+b))) (Beddington 1975; DeAngelis et al. 1975). Ces premières représentent un compromis entre la dépendance envers la proie et celle envers le ratio, ce ratio (N/P) étant replacé par la forme N/P<sup>m</sup>, m étant un paramètre situé dans l'intervalle [0,1] qui contrôle une transition continue entre la dépendance envers la proie et celle envers le ratio. L'équation Beddington-DeAngelis, à son tour, incorpore N autant que P, sans toutefois contenir un terme de ratio pure entre ces variables.

Du côté des réponses fonctionnelles (Holling 1959), celles-ci peuvent être de Type I (relation linéaire; y=ax ), de Type II (croissance de forme logarithmique; y=ax l(x+b) ), de Type III (croissance de forme logistique; y=ax² l(x²+b) ) ou de Type IV (croissance vers un sommet, suivi d'une décroissance), x étant la variable indépendante, c'est-à-dire, N dans le cas de la dépendance envers la proie et N/P dans le cas de la dépendance envers le ratio. À l'exception des espèces qui s'alimentent par microphagie suspensivore, les réponses fonctionnelles de Type I sont rares dans la nature (Jeschke et al. 2004). Celles de Type II indiquent une saturation des prédateurs, tandis que le Type III combine la saturation des prédateurs à des hautes populations de proie avec une courbe de croissance sigmoïde lorsque la proie se fait rare, propriétés qui peuvent découler de la saturation des refuges limités offerts aux proies par l'environnement, ou bien être une indication de sélection de proie ou d'apprentissage de la part du prédateur. Les réponses de Type IV se retrouvent auprès des systèmes qui incluent une proie dangereuse pour le
prédateur (Líznarová et Pekár 2013), dans lequel cas des hautes densités de proie peuvent au contraire freiner la prédation.

Certaines équations (telles les équations Hassell-Varley, de dépendance envers la proie et de dépendance envers le ratio) peuvent être combinées avec tous les types de réponse, tandis que d'autres (telle l'équation Beddington-DeAngelis) ne se prêtent pas à de telles modifications mais présentent néanmoins un comportement qui s'apparente à l'une ou plusieurs des formes de Holling. Jeschke et al. (2002) offrent une revue approfondie des équations de réponse fonctionnelle.

Aucune de ces équations présentement disponibles n'est cependant entièrement satisfaisante d'un point de vue théorique. La dépendance au ratio n'est pas réaliste à de basses densités de la proie, lorsque la rareté de celle-ci forcera les prédateurs à passer la majorité de leur temps à la recherche de la proie et entraînera donc une diminution du taux de prédation par prédateur, et ce, même si les prédateurs eux-mêmes se font assez rares pour produire un ratio élevé de proie par prédateur (Abrams et Ginzburg 2000). Les équations dépendantes du ratio produisent également des résultats qui ont par le passé été décrits comme « inhabituels », notamment, une disponibilité infinie de proie par prédateur lorsque la densité des prédateurs approche le zéro (Gleeson 1994). La dépendance à la proie, de son côté, produit elle-même des résultats non réalistes à de basses valeurs du ratio de proie au prédateur, car les prédateurs devront obligatoirement se partager les proies disponibles (Abrams et Ginzburg 2000). En outre, des études expérimentales ont observé des réponses fonctionnelles se situant entre les extrêmes de la dépendance à la proie et la dépendance au ratio (Schenk et al. 2005).

Le besoin d'une équation qui incorpore les processus de la dépendance au ratio autant que de celle à la proie se fait donc ressentir. Les équations alternatives existantes qui se posent en tant qu'intermédiaires entre la dépendance à la proie et la dépendance au ratio, cependant, le garantissent pas nécessairement une résolution cohérente d'un point de vue mathématique envers les problèmes associés à chaque extrême. À titre d'exemple, la transition entre dépendance au ratio et dépendance à la proie au sein des formes Hassell-Varley dépend d'un paramètre fixe, et non des populations de la proie ou du prédateur; cependant, les considérations théoriques évoquées au cours du paragraphe précédent suggèrent que les niveaux populationnels jouent un rôle central envers l'applicabilité de l'une ou de l'autre de ces formes. En outre, d'autres ont même suggéré de changer d'équation selon la nature de l'interférence entre les prédateurs du système (Skalski et Gilliam 2001).
Dans ce contexte, une liste de critères clairs envers le comportement attendu d'une
équation de prédation utile, surtout dans le contexte des modèles à temps discret,
contribuerait au développement plus standardisé d'équations de prédation qui demeurent
fiables peu importe les valeurs de la densité populationnelle ou du ratio entre la proie et le
prédateur. Une telle liste de critères fut proposée par Berryman et al. (1995), mais la majorité
de ceux-ci s'appliquent au développement de modèles complets des populations de prédateur
et de proie et non à la fonction de prédation en soi. Figurent parmi les critères qui
s'appliquent également à la fonction de prédation l'idée que les prédateurs « doivent avoir un
appétit limité » et que le potentiel de la compétition intraspécifique doit être pris en compte,
critères qui disqualifient d'emblée les équations dépendantes de la proie. Cette liste de
critères ne s'applique en outre qu'aux taux de prédation par prédateur et n'inclut pas les
critères importants en matière de taux de prédation total sous des conditions de population
extrême (soit élevée ou basse) du prédateur ou de la proie.

Malgré le débat de longue date à ce sujet, la plupart de ces problèmes n'ont toujours pas
trouvé de solution satisfaisante (Abrams 2015), et les études théoriques et pratiques de la
prédation continuent d'utiliser des équations qui arborent les lacunes théoriques décrites
auparavant. Même les équations de dépendance à la proie et de réponse fonctionnelle de Type
I, dont plusieurs études on démontré l'infériorité en comparaison aux alternatives (quoique
imparfaites elles aussi) dans le contexte de la majorité des systèmes écologiques (Abrams et
Ginzburg 2000; Jeschke et al. 2004), trouvent toujours application fréquente auprès de
nombre de modèles écologiques et de dynamiques des populations (Chesson et Kuang 2008;
Sanchez et al. 2018; Sanders et al. 2018; Imbert et al. 2020).

Ces lacunes théoriques existent parce que la majorité des ces équations furent construites
selon des suppositions (telles des densités de proie et de prédateur non extrêmes) qui ne
tiennent fréquemment pas lors de l'application de l'équation à la modélisation à temps discret
des conditions naturelles, tel est souvent le cas lors d'études de modélisation prédictive ou
d'expériences sur le terrain. Celles-ci sont des erreurs structurelles, lesquelles ne peuvent
donc pas être contournées à l'aide du paramétrage ou de la calibration. Les prévisions
générées par ces équations sous des conditions extrêmes, tel est fréquemment le cas en
agriculture ou auprès des autres systèmes gérés par les humains, courent donc le risque d'être
invalides. À titre d'exemple, les équations qui se reposent sur des suppositions qui ne tiennent
pas sous des conditions de haute densité du prédateur (Gleeson 1994; Abrams et Ginzburg
2000) ont peu de chances d'offrir des prévisions exactes au sujet du contrôle biologique par
inondation, lorsque des grandes quantités de prédateurs ou de parasitoïdes seront relâchées dans le champ. Ces enjeux demeurent une barrière importante envers l'application des équations de proie et de prédateur aux questions de modélisation écologique et agroécologique. Tel qu'observé par Kratina et al. (2009), « le choix d'une équation fonctionnelle appropriée est central à la génération de prévisions exactes des dynamiques des réseaux trophiques. Lorsque la dépendance au ratio est incorporée auprès des modèles de prédateur et de proie, leur stabilité se voit généralement améliorée » (traduction de l'auteur de la thèse). L'utilisation d'équations plus exactes pour la représentation des réponses fonctionnelles est donc un enjeu central dans le contexte de la fiabilité du comportement des modèles écologiques.

Cet article vise à répondre à cette question ouverte auprès de la littérature et propose comme objectifs 1) la présentation d'une liste de critères standards que toutes les équations de prédation utiles envers une application à temps discret devraient observer, 2) une démonstration que la majorité des équations disponibles auprès de la littérature échouent à l'une ou à plusieurs des ces conditions et 3) la présentation d'une nouvelle équation (nommée « Kovai ») de modélisation des taux de prédation qui répond favorablement à tous ces critères. L'article conclut avec 4) une évaluation comparative de toutes les équations révisées envers une grande gamme d'ensembles de données de prédation provenant de la littérature afin de démontrer les avantages d'ordre empirique offerts par la formulation améliorée de l'équation de prédation. Cet recherche est la première à réviser de manière exhaustive les limitations des équations de prédation existantes et à offrir une amélioration à ces équations qui demeure consistante d'un point de vue mathématique tout en performant de manière favorable dans le contexte théorique aussi bien qu'empirique.

## 11.3 Critères standards envers les fonctions de prédation

Les équations de prédation peuvent être définies soit selon le taux de prédation par prédateur par unité de temps, $f(N,P)$, soit selon le taux de prédation totale par unité de temps, $g(N,P)=P*f(N,P)$, $N$ et $P$ étant, respectivement, les densités de la proie et du prédateur. Les sections qui suivent proposent une liste de critères standards pour l'évaluation des équations de prédation dans le contexte des simulations à temps discret et seront suivies de la proposition d'une équation qui répond à chacun de ces critères. Lors des sections qui suivent, $\alpha$ indique le taux de prédation maximal par prédateur, c'est-à-dire, la quantité de
proie par unité de temps qu’un prédateur individuel pourrait consommer s’il ne devait perdre aucun temps à la recherche de la proie.

À titre de comparaison, les équations de prédation populaires auprès de la littérature furent également évaluées selon ces critères. Figurent parmi celles-ci les équations dépendantes de la proie et dépendantes du ratio, de même que les formes intermédiaires de Hassell-Varley (Hassell et Varley 1969), chacune desquelles fut combinée aux formes de réponse fonctionnelle de Type I, II et III de Holling (Holling 1959; Abrams et Ginzburg 2000). Trois autres équations qui ne s’adonnent pas aux modifications selon les formes de réponse fonctionnelle de Holling, soit l’équation de Beddington-DeAngelis (Beddington 1975; DeAngelis et al. 1975), l’équation de Crowley-Martin (Crowley et Martin 1989) et celle de Royama-Rogers (Royama 1971; Rogers 1972; Kfir 1983) furent également évaluées. Le Tableau 11.2 présente une définition formelle de chaque équation.

### 11.3.1 Comportement à haute densité de la proie

Lorsque la densité de la proie augmente, la prédation par prédateur devrait, tout au plus, approcher le taux de prédation maximal propre à l'espèce du prédateur. En termes mathématiques, $0 \le \lim_{N \to \infty} f(N, P) \le \alpha$ , et, par conséquent, $0 \le \lim_{N \to \infty} g(N, P) \le \alpha P$ .

### 11.3.2 Comportement à basse densité de la proie

Lorsque la densité de la proie approche 0, la prédation totale, de même que la prédation par prédateur, devrait elle aussi approcher 0; c'est-à-dire, $\lim_{N \to 0} f(N, P) = \lim_{N \to 0} g(N, P) = 0$ .

### 11.3.3 Comportement à haute densité du prédateur

Lorsque la densité du prédateur augmente, la prédation totale devrait augmenter de manière à s'approcher de l'asymptote de la population totale de la proie (ou possiblement vers une valeur inférieure si certaines niches dans l'habitat offrent à un nombre limité de proie une protection complète contre la prédation, ou bien si les prédateurs interfèrent entre eux). De même, la disponibilité de proie par prédateur devrait diminuer selon la compétition entre les prédateurs. En termes mathématiques, $\lim_{P \to \infty} g(N, P) \le N$ , tandis que $\lim_{P \to \infty} f(N, P) = 0$ .
### 11.3.4 Comportement à basse densité du prédateur

Lorsque la densité du prédateur diminue, au contraire, la prédation totale devrait elle aussi approcher le 0; c'est-à-dire, $\lim_{P \to 0} g(N,P)=0$. Tandis que la prédation par prédateur peut varier, elle ne devrait en aucun cas dépasser le taux de prédation maximal propre à l'espèce du prédateur; c'est-à-dire, $0 \le \lim_{P \to 0} f(N,P) \le \alpha$.

### 11.3.5 Autres considérations

La majorité des équations de réponse fonctionnelle supposent que la forme de celle-ci dépendra de l'abondance de la proie (soit de sa densité, soit du ratio de la proie au prédateur) et que le taux de prédation ne sera pas une fonction constante de l'abondance de la proie mais au contraire changera de façon non linéaire selon cette dernière. Ce phénomène peut s'expliquer par la saturation des prédateurs ou bien par l'existence de niches dans l'habitat qui offrent une protection partielle contre la prédation. Deux des trois formes de réponse fonctionnelle de Holling (Type II et Type III) qui ont par la suite connu une grande popularité prennent justement compte de ces phénomènes.

Toutefois, et même si l'on suppose que les proies additionnelles ajoutées à un système seront plus vulnérables à la prédation que les proies qui y figuraient avant (possiblement en raison de la saturation des niches protégées), l'addition de nouveaux individus de proie au système ne devrait jamais augmenter le risque de prédation instantané envers les proies déjà présentes. C'est-à-dire, l'ajout d'une unité de proie de peut pas augmenter la prédation totale par une valeur supérieure à une unité de proie; en termes mathématiques, $0 \le \frac{dg}{dN} \le 1$.

### 11.3.6 L'équation Kovai

Tel que démontré dans le Tableau 11.2, la majorité des équations présentement utilisées pour la modélisation de la prédation ne répondent pas à tous les critères énoncés ci-dessus. Cette section présente une nouvelle équation, que nous proposons de nommer selon le nom de la ville où la recherche menant à sa découverte fut effectuée.

$$ 
\begin{aligned}
f(N,P) &= \alpha (1 - e^{-Nu/(\alpha P)}) \\
u &= 1 - \frac{b}{N} (1 - e^{-N/b})
\end{aligned}
$$
La première partie de l'équation Kovai considère le ratio entre le prédateur et la proie et divise la quantité de proie disponible entre les prédateurs présents. En ce sens, l'équation Kovai ressemble à une équation dépendante du ratio. Cependant, au lieu d'utiliser la densité de la proie lors du calcul de ce ratio, une mesure de disponibilité effective de la proie, u, est prise en compte. Cette mesure, située 0 et 1, représente la difficulté associée à la recherche de la proie par le prédateur lorsque la proie se fait rare ou lorsqu'elle peut se cacher dans des refuges. Ce faisant, elle offre une réponse fonctionnelle similaire à celle des équations

Holling de Type III mais garde l'avantage que la valeur de la dérivée $$ \frac{d}{dN} uN = 1 - e^{-N/b} $$ demeure dans l'intervalle $[0, 1)$ tant que b est non négatif, ce qui assure que la disponibilité relative de la proie n'augmentera jamais plus rapidement que la population de la proie elle-même. Le paramètre b peut être considéré soit comme la capacité maximale des refuges pour la proie, étant donné que $$ \lim_{N \to \infty} (N - Nu) = \lim_{N \to \infty} [N - N + b(1 - e^{-N/b})] = b $$ , soit comme la densité de la proie à laquelle $e^{-1} \approx 36,8\%$ de la population de la proie est accessible au prédateur.

Lorsque le paramètre b prend un valeur proche de 0, l'équation s'apparente à la réponse fonctionnelle Holling Type II; par conséquent, ce paramètre contrôle la forme de la réponse fonctionnelle de l'équation Kovai. En même temps, la dépendance du facteur u sur N au lieu du ratio N/P assure que l'équation Kovai prendra compte d'une pénurie de proie de manière similaire aux équations dépendantes de la proie.

## 11.4 Méthodes (étude empirique)

12 études expérimentales rapportées dans la littérature et comportant des données de réponse fonctionnelle de la prédation furent identifiées afin d'obtenir un total de 29 différents ensembles de données (Tableau 11.1). Lorsque les données étaient uniquement disponibles auprès d'une figure et non d'un tableau, celles-ci furent digitalisées à l'aide de la version en ligne du logiciel WebPlotDigitizer (Rohatgi 2019). Toutes les analyses furent effectuées en Python 3,7 (Python Software Foundation (PSF) 2019) avec l'aide des libraires Numpy (van der Walt et al. 2011), Pandas (McKinney 2010) et PyMC3 (Salvatier et al. 2016 p. 3) pour les analyses mathématiques, de même que Seaborn (Michael Waskom et al. 2018) et Matplotlib (Hunter 2007) pour la visualisation des résultats. Toute donnée et code nécessaire à la
reproduction des analyses et à la génération des figures présentées dans cet article est
disponible auprès du matériel supplémentaire.

<table><caption>Tableau 11.1: Études et bases de données utilisées pour les analyses. DP : dépendance de la proie, DR : dépendance du ratio; HV : Hassell-Varley; CM : Crowley-Martin; RR : Royama-Rogers; BD : Beddington-DeAngelis; A : autre équation pas considérée par la présente étude. Les chiffres romains indiquent la forme de la réponse fonctionnelle Holling. Seules les équations identifiées en raison de leur bon ajustement ou bien recommandées par les auteurs de l'étude correspondante sont incluses dans la colonne de droite.</caption><thead><tr><td>Référence</td><td>Bases de données</td><td>Système écologique</td><td>Équations sélectionnées</td></tr></thead><tbody><tr><td>Kratina et al. (2009)</td><td>1</td><td>Plathelminthe benthique-paramécium</td><td>BD*, HV III, DR III</td></tr><tr><td>Fernandez-Maldonado et al. (2017)</td><td>1</td><td>Cannibalisme de Nabis pseudoferus</td><td>A I</td></tr><tr><td>Líznarová et Pekár (2013)</td><td>5</td><td>Araignée-fourmis (diverses espèces)</td><td>A IV</td></tr><tr><td>Schenk et al. (2005)</td><td>1</td><td>Polistes dominulus-Cassida rubiginosa</td><td>HV III</td></tr><tr><td>Schenk et Bacher (2002)</td><td>2</td><td>Polistes dominulus-Cassida rubiginosa</td><td>A III</td></tr><tr><td>Long et al. (2012)</td><td>4</td><td>Cannibalisme du Crabe royal du Kamtchatka</td><td>DP II</td></tr><tr><td>Kfir (1983)</td><td>1</td><td>Parasitoïde d'œuf de Phthorimaea operculella</td><td>RR</td></tr><tr><td>Edwards (1961)</td><td>3</td><td>Parasitisme (diverses espèces)</td><td>Aucune</td></tr><tr><td>Eveleigh et Chant (1982)</td><td>4</td><td>Acariens prédateurs</td><td>HV</td></tr><tr><td>Mertz et Davies (1968)</td><td>1</td><td>Cannibalisme du tribolium rouge de la farine</td><td>A</td></tr><tr><td>Taylor (1988)</td><td>4</td><td>Ectoparasitisme (diverses proies)</td><td>RR **</td></tr><tr><td>Uttley (1980)***</td><td>2</td><td>Zygoptère/notonecte-cladocère</td><td>CM, HV II***</td></tr></tbody></table>

Chaque équation présentée dans le Tableau 11.2 fut calibrée envers chacun des ensembles
de données par inférence bayésienne avec distributions a priori peu informatives selon

<footnote>* Légèrement modifiée afin d'obtenir une réponse fonctionnelle de Type III.</footnote>

<footnote>** Légèrement modifiée.</footnote>

<footnote>*** Données et équations sélectionnées de la figure de Skalski et Gilliam (2001).</footnote>
l'algorithme d'échantillonnage sans demi-tour de Monte Carlo hamiltonien (Hoffman et Gelman 2014).

Dans le cas des expériences où les densités des proies furent maintenues à un niveau constant (opportunités de prédation théoriquement infinies), la réponse fonctionnelle fut modélisée selon la distribution poisson suivante :

$$y \sim \text{Poisson}(\mu * n)$$

Où y représente le nombre total d'événements de prédation observés auprès de toutes les répétitions de l'expérience, n est le nombre de répétitions, et μ est déterminé selon la spécification de chaque équation.

Dans le cas (nettement plus fréquent) où les proies ne furent pas remplacées afin de garder une densité constante au cours de l'expérience, les données furent modélisées selon une distribution binomiale :

$$y \sim \text{Binomiale}(n=N*n, p=\min(\mu/N, 0,999))$$

Où N est la quantité de proie initial et les autres paramètres gardent les mêmes définitions que lors de l'équation poisson ci-dessus. Notons que p dut être limité afin d'éviter des erreurs de valeurs hors bornes avec certaines équations fonctionnelles.

Suite à la calibration, la performance de chaque équation envers chaque ensemble de données fut ordonnée selon le Critère d'information largement applicable (CILA) (Watanabe 2013; Vehtari et al. 2017), lequel pénalise les équations qui contiennent un plus grand nombre de paramètres afin d'éviter le surapprentissage. Son utilisation fut recommandée au lieu de celle du Critère d'information d'Akaike dans le contexte des analyses bayésiennes (Gelman et al. 2014); des valeurs plus petites du CILA indiquent un meilleur ajustement du modèle.

# 11.5 Résultats et discussion

## 11.5.1 Étude théorique

Tel que présenté par le Tableau 11.2, aucune des équations présentement disponibles auprès de la littérature ne répond à tous les critères proposés dans la section 11.3. Ces résultats remettent en question la fiabilité de ces équations envers la génération de prévisions exactes dans le contexte des systèmes dont les populations de prédateur ou de proie dévient de manière significative des conditions d'équilibre statique. La nouvelle équation proposée par cette étude, cependant, donna des résultats satisfaisants face à toutes les conditions.
Seules les critères de comportement lorsque la population de la proie approche le 0 (c'est-
à-dire, que la prédation totale et la prédation par prédateur devraient approcher 0 elles aussi)
furent respectés par toutes les équations. Chacune des autres condition fut brisée par une ou
plusieurs équations.

Tel que prévu, les comportements des équations dépendantes de la proie ne sont pas
favorables dans des conditions de haute densité des prédateurs, conditions sous lesquelles les
formes de dépendance au ratio et de Hassell-Varley offrirent des meilleurs résultats. Les
réponses fonctionnelles de Type I, cependant, génèrent des comportements problématiques
sous des conditions de densité élevée de la proie, et ce, peu importe la classe de l'équation en
question.

Plus problématique, cependant, sont les cas où les taux de prédation sous des conditions de
densité élevée du prédateur soit excédent tout simplement la population totale de la proie (p.
ex., $\lim_{P \to \infty} (g)=\infty$ ), soit risquent de l'excéder selon le paramétrage de l'équation (p. ex.,

$\lim_{P \to \infty} (g) = \alpha N/b$). En fin de compte, seule la nouvelle équation proposée (Kovai) présente un comportement consistant envers le taux de prédation total sous des conditions de haute densité du prédateur. En outre, aucune des équations de la littérature ne prévoit correctement le taux de croissance de la prédation lorsqu'une proie additionnelle est ajouté au système. Dans de nombreux cas, l'ajout d'une unité de proie pourrait engendrer une augmentation beaucoup plus importante dans le taux de prédation total, ce qui indiquerait que l'ajout d'un individu de proie au système pourrait augmenter de manière importante la prédation auprès des autres individus de proie qui étaient déjà présentes dans ce système. Un tel phénomène paraît improbable dans la réalité et, lors d'études de modélisation, pourrait mener à la surestimation des taux de prédation. Seule l'équation Kovai offrit un comportement cohérent à cet égard.
<table><caption>Tableau 11.2: Comparaison du comportement de différentes équations de prédation populaires sous des conditions extrêmes. T indique la forme de réponse fonctionnelle Holling (I, II ou III). f(N,P) représente le taux de prédation par prédateur et g(N,P) = P * f(N,P) représente le taux de prédation total. α représente le taux de prédation maximal par prédateur, N et P sont les densités de la proie et du prédateur, T<sub>h</sub> est le temps de manutention par proie consommée et les autres variables sont des constantes empiriques. α, b, c ∈ [0, ∞); m ∈ [0, 1]. Les cellules ombragées indiquent des comportements indésirables. Voir la section 11.11 (Appendice A : Dérivation des limites et dérivées non triviales) pour les détails des calculs menant aux résultats présentés ci-dessous.</caption><thead><tr><th>Équation</th><th>f(N,P)</th><th>lim(f)<sub>P→0</sub></th><th>lim(f)<sub>P→∞</sub></th><th>lim(f)<sub>N→0</sub></th><th>lim(f)<sub>N→∞</sub></th><th>lim(g)<sub>P→0</sub></th><th>lim(g)<sub>P→∞</sub></th><th>lim(g)<sub>N→0</sub></th><th>lim(g)<sub>N→∞</sub></th><th>dg/dN</th><th>Référence</th></tr></thead><tbody><tr><td>Équation idéale</td><td></td><td>[0, α]</td><td>0</td><td>0</td><td>[0, α]</td><td>0</td><td>[0, N]</td><td>0</td><td>[0, αP]</td><td>[0, 1]</td><td></td></tr><tr><td rowspan="3">Dépendance proie</td><td>T I</td><td>αN</td><td>αN</td><td>αN</td><td>0</td><td>∞</td><td>0</td><td>∞</td><td>0</td><td>∞</td><td>[0, ∞)</td><td rowspan="3">(Holling 1959; Abrams et Ginzburg 2000)</td></tr><tr><td>T II</td><td>αN / (N+b)</td><td>αN / (N+b)</td><td>αN / (N+b)</td><td>0</td><td>α</td><td>0</td><td>∞</td><td>0</td><td>αP</td><td>[0, ∞)</td></tr><tr><td>T III</td><td>αN<sup>2</sup> / (N<sup>2</sup>+b)</td><td>αN<sup>2</sup> / (N<sup>2</sup>+b)</td><td>αN<sup>2</sup> / (N<sup>2</sup>+b)</td><td>0</td><td>α</td><td>0</td><td>∞</td><td>0</td><td>αP</td><td>[0, ∞)</td></tr><tr><td rowspan="3">Dépendance ratio</td><td>T I</td><td>αN/P</td><td>∞</td><td>0</td><td>0</td><td>∞</td><td>αN</td><td>αN</td><td>0</td><td>∞</td><td>[0, ∞)</td><td rowspan="3">(Holling 1959; Arditi et Ginzburg 1989)</td></tr><tr><td>T II</td><td>αN/P / (N/P+b)</td><td>α</td><td>0</td><td>0</td><td>α</td><td>0</td><td>αN / b</td><td>0</td><td>αP</td><td>[0, ∞)</td></tr><tr><td>T III</td><td>α(N/P)<sup>2</sup> / ((N/P)<sup>2</sup>+b)</td><td>α</td><td>0</td><td>0</td><td>α</td><td>0</td><td>0</td><td>0</td><td>αP</td><td>[0, ∞)</td></tr><tr><td>Beddington-DeAngelis</td><td></td><td>αN / (N+cP+b)</td><td>αN / (N+b)</td><td>0</td><td>0</td><td>α</td><td>0</td><td>αN / c</td><td>0</td><td>αP</td><td>[0, ∞)</td><td>(Beddington 1975; DeAngelis et al. 1975)</td></tr><tr><td>Crowley-Martin</td><td></td><td>αN / (1+bN)(1+c[P-1])</td><td>αN / (1+bN)*(1-c)</td><td>0</td><td>0</td><td>αP / (b+bc(P-1))</td><td>0</td><td>αN / (1+bN)c</td><td>0</td><td>αP / (b+bc(P-1))</td><td>(-∞, ∞)</td><td>(Crowley et Martin 1989)</td></tr><tr><td>Royama-Rogers</td><td></td><td>N(1-e<sup>-αP</sup> / (1+αT<sub>h</sub>N))</td><td>0</td><td>N</td><td>0</td><td>P / T<sub>h</sub></td><td>0</td><td>∞</td><td>0</td><td>P<sup>2</sup> / T<sub>h</sub></td><td>[0, ∞)</td><td>(Royama 1971; Rogers 1972; Kfir 1983)</td></tr></tbody></table>
<table><thead><tr><th>Équation</th><th>f(N,P)</th><th>lim(f)<br>P→0</th><th>lim(f)<br>P→∞</th><th>lim(f)<br>N→0</th><th>lim(f)<br>N→∞</th><th>lim(g)<br>P→0</th><th>lim(g)<br>P→∞</th><th>lim(g)<br>N→0</th><th>lim(g)<br>N→∞</th><th>dg/dN</th><th>Référence</th></tr></thead><tbody><tr><td>T I</td><td>α(N/P<sup>m</sup>)</td><td>αN, m=0<br>∞, m&gt;0</td><td>αN, m=0<br>0, m&gt;0</td><td>0</td><td>∞</td><td>0, 0≤m&lt;1<br>αN, m=1</td><td>∞, 0≤m&lt;1<br>N, m=1</td><td>0</td><td>∞</td><td>[0,∞)</td><td></td></tr><tr><td>Hassell-<br>Varley</td><td>T II<br>α(N/P<sup>m</sup>)<br>(N/P<sup>m</sup>) + b</td><td>αN/N+b, m=0<br>α, m&gt;0</td><td>αN/N+b, m=0<br>0, m&gt;0</td><td>0</td><td>α</td><td>0</td><td>αN/b, m=1</td><td>0</td><td>αP</td><td>[0,∞)</td><td>(Hassell et<br>Varley 1969)</td></tr><tr><td>T III</td><td>α(N/P<sup>m</sup>)<sup>2</sup><br>(N/P<sup>m</sup>)<sup>2</sup> + b</td><td>αN<sup>2</sup>/N<sup>2</sup>+b, m=0<br>α, m&gt;0</td><td>αN<sup>2</sup>/N<sup>2</sup>+b, m=0<br>0, m&gt;0</td><td>0</td><td>α</td><td>0</td><td>αN/b, m=0,5<br>0, 0.5&lt;m≤1</td><td>0</td><td>αP</td><td>[0,∞)</td><td></td></tr><tr><td>Kovai</td><td>α(1 - e<sup>-Nu/(αP)</sup>)<br>u=1 - b/N(1 - e<sup>-N/b</sup>)</td><td>α</td><td>α</td><td>0</td><td>α</td><td>0</td><td>Nu∈[0, N]</td><td>0</td><td>αP</td><td>[0,1)</td><td></td></tr></tbody></table>
11.5.2 Étude empirique

Les bases de données utilisées pour l'étude empirique (Tableau 11.1) représentent un
grande gamme de relations entre proie et prédateur, de même que de catégories d'équation et
de réponses fonctionnelles (malgré la prédominance des réponses de Type II et III). La Figure
11.1 illustre la versatilité de l'équation Kovai face à des données qui présentent de différentes
formes fonctionnelles. Cette équation, selon son paramétrage, peut générer une réponse de
Type II ou de Type III (mais toutefois pas de Type IV) envers la densité de la proie. Du même
trait, l'impact de la compétition entre prédateurs auprès du taux de prédation total est lui aussi
bien représenté.

<figure><img src="image_45.png" /></figure>

Figure 11.1: Illustration des différentes formes de réponse fonctionnelle qui peuvent découler de l'équation Kovai. a) Réponse de forme Holling Type III (Schenk et Bacher 2002, données de l'année 2000); b) réponse de forme Holling Type II (Fernandez-Maldonado et al. 2017); c) tel que prévu, ajustement moins qu'optimal de l'équation face aux données de forme Holling Type IV (Líznarová et Pekár 2013, données de Zodarion-Tetramorium); d) illustration du partage les proies entre prédateurs lorsque ces derniers commencent à interférer entre eux (Uttley 1980, figure 4.3.1). N : densité de la proie; P : densité du prédateur.
La Figure 11.2 présente les résultats de différentes équations lorsqu'appliquées aux données de prédation de Kratina et al. (2009), y compris trois équations recommandées par ces auteurs de même que l'équation Kovai.

À de basses valeurs de P, lorsque le ratio de proie au prédateur est suffisamment élevé pour éviter la compétition entre prédateurs, le modèle dépendant du ratio ne sut représenter l'impact de la rareté absolue des proies envers les taux de prédation et par conséquence confond les effets de différentes densités de la proie (Figure 11.2, a). Les équations de forme intermédiaire telles que Hassell-Varley Type III, Beddington-DeAngelis et Kovai (voir figure Figure 11.2, b, c et d) surent quant à elles distinguer entre la densité de la proie et celle du prédateur. Notons en passant que les modèles dépendants de la proie n'eurent aucun succès envers la modélisation de cet ensemble de données en raison de leur incapacité de différencier entre les différentes valeurs de P et ne sont donc pas incluses auprès de la figure.

<figure><img src="image_46.png" /><figcaption>Figure 11.2: Présentation des observations de la réponse fonctionnelle (proie consommée par prédateur) pour de différentes densités de prédateur (Stenostomum virginanum) et de proie (Paramecium aurelia), de même que des prévisions des modèles calibrés. La ligne représente la moyenne du modèle, et la zone ombragée l'écart type de la distribution a</figcaption></figure>
posteriori du modèle. a) Équation dépendante du ratio de Type III; b) Équation Hassell-Varley de Type III; c) Équation Beddington-DeAngelis; d) Équation Kovai. N indique la densité de la proie. Les données proviennent de Kratina et al. (2009).

La Figure 11.3 présente l'ajustement des diverses catégories d'équation envers chacune des bases de données utilisées lors de cette étude. Nous nous attendions à ce que l'équation Kovai, étant donné sa capacité de s'ajuster aux formes de réponse fonctionnelle de Type II ou de Type III, présente une performance plus consistante que celle des équations de Type II ou de Type III pures (soit dépendantes du ratio, soit de la proie) envers des ensembles de données présentant de diverses formes de réponse fonctionnelle. Le degré de ressemblance d'un ensemble de données à une réponse fonctionnelle de Type II fut estimé selon le classement moyen de l'ajustement CILA des équations de Type II (en comparaison à l'ajustement des autres équations) lorsqu'appliquées à cet ensemble de données. La Figure 11.3 compare le classement de diverses équations au classement des équations de Type II pour chacun des différents ensembles de données. Les classements affichant une valeur numérique plus petite indiquent un meilleur ajustement de l'équation d'intérêt en comparaison à l'ajustement des autres équations.

<figure><img src="image_47.png" /></figure>

Figure 11.3: Classement de diverses équations (un classement moins élevé indiquant un meilleur ajustement de l'équation en comparaison aux autres) selon le classement moyen des équations de Type II pour le même ensemble de données. Chaque point représente le classement d'une équation ou d'un groupe d'équations (Type III, BD ou KV) lorsque calibrée selon un ensemble de données expérimentales (un classement inférieur étant favorable). Les lignes indiquent le modèle de régression linéaire pour chaque équation ou groupe d'équations, les zones ombragées indiquant les intervalles de confiance de 95%. Type II inclut les équations de Hassell-Varley, de dépendance à la proie et de dépendance au ratio
combinées avec une réponse de Holling Type II; Type III inclut les mêmes équations avec une réponse de forme Holling Type III; BD signifie l'équation Beddington-DeAngelis; KV signifie l'équation Kovai.

Tel que prévu, le classement des équations de Type III s'améliore au fur et à mesure que les données s'éloignent d'une forme fonctionnelle de Type II. L'ajustement de l'équation Beddington-DeAngelis, dont la forme ressemble à celle de Holling Type II, s'améliore également auprès des données qui ressemblent plus à une réponse fonctionnelle de Type II.

L'équation Kovai, tout au contraire, présente un classement relativement consistent envers les différents ensembles de donnés, peut importe leur degré de ressemblance à une réponse fonctionnelle de Type II. Cet avantage découle vraisemblablement du paramétrage de l'équation Kovai, qui lui permet de modéliser les réponses fonctionnelles de Type II tout autant que celles de Type III. Étant donné que l'équation Kovai suppose que la disponibilité de la proie augmentera avec la densité de la même, il demeure important de mentionner que cette équation n'est pas apte à la représentation des réponses fonctionnelles de Type IV, lorsqu'une augmentation auprès de la densité d'une proie dangereuse diminue au contraire le potentiel de chasse des prédateurs.

En conclusion, l'équation Kovai présente plusieurs améliorations clefs en comparaison aux équations de prédation existantes auprès de la littérature dans le contexte d'une simulation à temps discret. 1) Sa contribution principale est une meilleure cohérence théorique, en tout particulier, sa prévision de la prédation totale sous des conditions de population élevée de la proie ou du prédateur, de même que sa représentation du taux de croissance de la prédation totale au fur et à mesure que des individus de proie sont ajoutés au système. 2) L'équation Kovai combine également le « meilleur des deux mondes » entre les concepts de la dépendance au ratio et de la dépendance à la proie, et ce, de manière cohérente du point de vue mathématique, en ajustant pour la compétition entre prédateurs (partage de la proie entre les prédateurs lorsque la population de la proie devient limitante) tout en gardant la capacité de distinguer entre des différentes densités absolues de la population de la proie même lorsque le ratio de proie au prédateur demeure constant. 3) De même, l'équation offre, selon son paramétrage, une transition lisse (et cohérente du point de vue mathématique) entre les réponses fonctionnelles de Type II et de Type III.

## 11.6 Conclusion

Le domaine de l'écologie est toujours en quête d'une équation qui représente les relations entre les proies et les prédateurs de manière cohérente du point de vue mathématique, et ce,
malgré la centralité de ces relations à toute modélisation écologique des dynamiques des populations. La présente recherche démontre que la grande majorité des équations en usage fréquent auprès des études écologiques d'ordre théorique de même qu'empirique ne répondent pas à plusieurs critères de base indispensables à la cohérence mathématique de ces équations. Nous proposons une liste de tels critères afin de faciliter l'évaluation des équations de prédation qui seront développées à l'avenir et, du même trait, proposons une nouvelle équation qui répond à chacun de ces critères. L'évaluation empirique de la nouvelle équation suggère qu'elle serait plus robuste envers les différences entre les formes de réponse fonctionnelle qui se retrouvent dans de différents contextes écologiques.

## 11.7 Remerciements

Cette recherche bénéficia du support d'une bourse Découverte du CRSNG (no. RGPIN-2015-05554) octroyée à Jan Adamowski, de même que d'une Bourse de 3e cycle du FRQNT, d'une Bourse du CRDI aux chercheurs candidats au doctorat (no. 107759-99906075-017), d'une Bourse d'études supérieures du Canada Vanier et d'un Supplément pour études à l'étranger Michael-Smith (Programme de bourses d'études supérieures du Canada) octroyées à Julien Malard.

## 11.8 Contributions des auteurs

JM effectua la recherche et écrivit l'article. JBN contribua à l'analyse mathématique. JA, HQM, NA et HTY supervisèrent le projet de recherche. Les auteurs ne déclarent aucun conflit d'intérêt.

## 11.9 Disponibilité des données

Tout code et donnée nécessaire à la reproduction des analyses et des figures présentées au cours du présent article sont disponibles auprès du matériel supplémentaire.

## 11.10 Références

Abrams PA. 2015. Why ratio dependence is (still) a bad model of predation. Biological Reviews 90: 794–814.

Abrams PA et Ginzburg LR. 2000. The nature of predation: Prey dependent, ratio dependent or neither? Trends in Ecology and Evolution 15: 337–41.
Arditi R et Ginzburg LR. 1989. Coupling in predator-prey dynamics: Ratio-Dependence. Journal of Theoretical Biology 139: 311-26.

Beddington JR. 1975. Mutual Interference Between Parasites or Predators and its Effect on Searching Efficiency. The Journal of Animal Ecology 44: 331.

Berryman AA, Gutierrez AP, et Arditi R. 1995. Credible, Parsimonious and Useful Predator-Prey Models: A Reply to Abrams, Gleeson, and Sarnelle. Ecology 76: 1980-5.

Chesson P et Kuang JJ. 2008. The interaction between predation and competition. Nature 456: 235-8.

Crowley PH et Martin EK. 1989. Functional Responses and Interference within and between Year Classes of a Dragonfly Population. Journal of the North American Benthological Society 8: 211-21.

DeAngelis DL, Goldstein RA, et O'Neill RV. 1975. A Model for Tropic Interaction. Ecology 56: 881-92.

Edwards RL. 1961. The Area of Discovery of Two Insect Parasites, Nasonia vitripennis (Walker) and Trichogramma evanescens Westwood, in an Artificial Environment. Can Entomol 93: 475-81.

Eveleigh ES et Chant DA. 1982. Experimental studies on acarine predator-prey interactions: the effects of predator density on prey consumption, predator searching efficiency, and the functional response to prey density (Acarina: Phytoseiidae). Can J Zool 60: 611-29.

Fernandez-Maldonado FJ, Gallego JR, Valencia A, et al. 2017. Cannibalism: Do risks of fighting and reprisal reduce predatory rates? Community Ecology 18: 87-96.

Gelman A, Hwang J, et Vehtari A. 2014. Understanding predictive information criteria for Bayesian models. Stat Comput 24: 997-1016.

Gleeson SK. 1994. Density Dependence is Better Than Ratio Dependence. Ecology 75: 1834-5.
Hassell MP et Varley GC. 1969. New Inductive Population Model for Insect Parasites and its Bearing on Biological Control. Nature 223: 1133-7.

Hoffman MD et Gelman A. 2014. The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research 15: 1593-623.

Holling CS. 1959. The Components of Predation as Revealed by a Study of Small-Mammal Predation of the European Pine Sawfly. Can Entomol 91: 293-320.

Hunter JD. 2007. Matplotlib: A 2D Graphics Environment. Comput Sci Eng 9: 90-5.

Imbert C, Papaïx J, Husson L, et al. 2020. Estimating population dynamics parameters of cabbage pests in temperate mixed apple tree-cabbage plots compared to control vegetable plots. Crop Protection 129: 105037.

Jeschke JM, Kopp M, et Tollrian R. 2002. Predator Functional Responses: Discriminating Between Handling and Digesting Prey. Ecological Monographs 72: 95-112.

Jeschke JM, Kopp M, et Tollrian R. 2004. Consumer-food systems: Why type I functional responses are exclusive to filter feeders. Biological Reviews of the Cambridge Philosophical Society 79: 337-49.

Kfir R. 1983. Functional response to host density by the egg parasite Trichogramma pretiosum. Entomophaga 28: 345-53.

Kratina P, Vos M, Bateman A, et Anholt BR. 2009. Functional responses modified by predator density. Oecologia 159: 425-33.

Líznarová E et Pekár S. 2013. Dangerous prey is associated with a type 4 functional response in spiders. Animal Behaviour 85: 1183-90.

Long WC, Popp J, Swiney KM, et Van Sant SB. 2012. Cannibalism in red king crab, Paralithodes camtschaticus (Tilesius, 1815): Effects of habitat type and predator density on predator functional response. Journal of Experimental Marine Biology and Ecology 422-423: 101-6.
McKinney W. 2010. Data Structures for Statistical Computing in Python. In: Proceedings of the 9th Python in Science Conference.

Mertz DB et Davies RB. 1968. Cannibalism of the Pupal Stage by Adult Flour Beetles: An Experiment and a Stochastic Model. Biometrics 24: 247.

Michael Waskom, Olga Botvinnik, Drew O'Kane, et al. 2018. seaborn. Zenodo.

Python Software Foundation (PSF). 2019. Python.

Rogers D. 1972. Random Search and Insect Population Models. Journal of Animal Ecology 41: 369-83.

Rohatgi A. 2019. WebPlotDigitizer.

Royama T. 1971. A comparative study of models for predation and parasitism. Res Popul Ecol 13: 1-91.

Salvatier J, Wiecki TV, et Fonnesbeck C. 2016. Probabilistic programming in Python using PyMC3. PeerJ Computer Science 2: e55.

Sanchez C, Gamez M, Burguillo FJ, et al. 2018. Comparison of predator-parasitoid-prey interaction models for different host plant qualities. Community Ecology 19: 125-32.

Sanders D, Thébault E, Kehoe R, et Veen FJF van. 2018. Trophic redundancy reduces vulnerability to extinction cascades. PNAS 115: 2419-24.

Schenk D et Bacher S. 2002. Functional response of a generalist insect predator to one of its prey species in the field. Journal of Animal Ecology 71: 524-31.

Schenk D, Bersier L-F, et Bacher S. 2005. An experimental test of the nature of predation: neither prey- nor ratio-dependent: Experimental test of functional response. Journal of Animal Ecology 74: 86-91.

Skalski GT et Gilliam JF. 2001. Functional responses with predator interference: viable alternatives to the Holling Type II model. Ecology 82: 3083-92.
Taylor AD. 1988. Host Effects on Functional and Ovipositional Responses of Bracon hebetor. The Journal of Animal Ecology 57: 173.

Uttley M. 1980. A laboratory study of mutual interference between freshwater invertebrate predators.

Vehtari A, Gelman A, et Gabry J. 2017. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput 27: 1413–32.

Walt S van der, Colbert SC, et Varoquaux G. 2011. The NumPy Array: A Structure for Efficient Numerical Computation. Comput Sci Eng 13: 22–30.

Watanabe S. 2013. A Widely Applicable Bayesian Information Criterion. Journal of Machine Learning Research 14: 867–97.

# 11.11 Appendice A : Dérivation des limites et dérivées non triviales

## 11.11.1 Dépendance à la proie Type II

Le modèle de dépendance à la proie de Holling Type II est spécifié comme suit (Holling 1959; Abrams et Ginzburg 2000) :

$$f(N, P) = \frac{\alpha N}{N + b}$$

$$g(N, P) = \frac{\alpha PN}{N + b}$$

Par conséquent :

$$\frac{dg}{dN} = \frac{\alpha P(N+b) - \alpha PN}{(N+b)^2} = \frac{\alpha Pb}{(N+b)^2}$$

Étant donné que toutes les variables sont des nombres non négatifs, l'expression

$\frac{Pb}{(N+b)^2}$ est elle aussi un nombre non négatif. Vue que $\alpha \in [0, \infty)$, l'expression

$\frac{dg}{dN} = \alpha \frac{Pb}{(N+b)^2}$ se situe donc elle aussi dans l'intervalle $[0, \infty)$.
11.11.2 Dépendance à la proie Type III

Le modèle de dépendance à la proie de Holling Type III se spécifie comme suit (Holling
1959; Abrams et Ginzburg 2000) :

$$
f(N, P) = \frac{\alpha N^2}{N^2 + b}
$$

$$
g(N,P) = \frac{\alpha PN^2}{N^2+b}
$$

Par conséquent :

$$
\frac{dg}{dN} = \frac{2 \alpha PN * (N^2 + b) - 2N(\alpha PN^2)}{(N^2 + b)^2} = \frac{2 \alpha PNb}{(N^2 + b)^2}
$$

Étant donné que toutes les variables sont des nombres non négatifs, l'expression

$$
2 \frac{PNb}{(N^2 + b)^2} \text{ est elle aussi un nombre non négatif. Vue que } \alpha \in [0, \infty) \text{, l'expression}
$$

$$
\frac{dg}{dN} = \alpha \frac{2 PNb}{(N^2 + b)^2} \text{ se situe donc dans l'intervalle } [0, \infty) \text{ .}
$$

**11.11.3 Dépendance au ratio Type II**

Le modèle de dépendance au ratio de Holling Type II se spécifie comme suit (Arditi et
Ginzburg 1989; Abrams et Ginzburg 2000) :

$$
f(N, P) = \frac{\alpha N/P}{N/P + b}
$$

$$
g(N, P) = \frac{\alpha N}{N/P + b}
$$

Par conséquent :

$$
\lim_{P \to 0} (f) = \lim_{P \to 0} \left[ \frac{\alpha N}{N + bP} \right] = \frac{\alpha N}{N} = \alpha
$$

$$
\lim_{P \to \infty} (f) = \lim_{P \to \infty} \left[ \frac{\alpha N}{N + bP} \right] = 0
$$

$$
\lim_{N \to 0} (f) = \frac{0}{0+b} = 0
$$
$$ \lim_{N \to \infty} (f) = \lim_{N \to \infty} \left[ \frac{\alpha / P}{1/P + b/N} \right] = \frac{\alpha / P}{1/P} = \alpha $$

$$ \lim_{P \to 0} (g) = \frac{\alpha N}{\infty} = 0 $$

$$ \lim_{P \to \infty} (g) = \frac{\alpha N}{0+b} = \frac{\alpha N}{b} $$

$$ \lim_{N \to 0} (g) = \frac{0}{0+b} = 0 $$

$$ \lim_{N \to \infty} (g) = \lim_{N \to \infty} \left[ \frac{\alpha}{1/P + b/N} \right] = \frac{\alpha}{1/P} = \alpha P $$

$$ dg/dN = \frac{\alpha[N/P + b] - \alpha N[1/P]}{[N/P + b]^2} = \frac{\alpha b}{[N/P + b]^2} $$

Étant donné que toutes les variables sont des nombres non négatifs, l'expression

$\frac{b}{[N/P + b]^2}$ représente elle aussi un nombre non négatif. Vue que $\alpha \in [0, \infty)$,

l'expression $\frac{dg}{dN} = \alpha \frac{b}{[N/P + b]^2}$ se situe donc dans l'intervalle $[0, \infty)$.

### 11.11.4 Dépendance au ratio Type III

Le modèle de dépendance au ratio de Holling Type III est spécifié comme suit (Arditi et Ginzburg 1989; Abrams et Ginzburg 2000) :

$$ f(N, P) = \frac{\alpha (N/P)^2}{(N/P)^2 + b} $$

$$ g(N, P) = \frac{\alpha N^2 / P}{(N/P)^2 + b} $$

Par conséquent :

$$ \lim_{P \to 0} (f) = \lim_{P \to 0} \left[ \frac{\alpha N^2}{N^2 + bP^2} \right] = \frac{\alpha N^2}{N^2} = \alpha $$
$$ \lim_{P \to \infty} (f) = \lim_{P \to \infty} \left[ \frac{\alpha N^2}{N^2 + bP^2} \right] = \frac{\alpha N^2}{\infty} = 0 $$

$$ \lim_{N \to 0} (f) = \frac{0}{0 + b} = 0 $$

$$ \lim_{N \to \infty} (f) = \lim_{N \to \infty} \left[ \frac{\alpha / P^2}{1 / P^2 + b / N^2} \right] = \frac{\alpha / P^2}{1 / P^2} = \alpha $$

$$ \lim_{P \to 0} (g) = \lim_{P \to 0} \frac{\alpha N^2}{[(N^2/P) + bP]} = 0 $$

$$ \lim_{P \to \infty} (g) = \lim_{P \to \infty} \frac{\alpha N^2}{P[(N/P)^2 + b]} = 0 $$

$$ \lim_{N \to 0} (g) = \frac{0}{0 + b} = 0 $$

$$ \lim_{N \to \infty} (g) = \lim_{N \to \infty} \frac{\alpha / P}{(1 / P)^2 + b / N^2} = \frac{\alpha / P}{(1 / P)^2} = \alpha P $$

$$ \frac{dg}{dN} = \frac{2 \alpha N/P [(N/P)^2 + b] - \alpha (N^2/P) 2 N/P^2}{[(N/P)^2 + b]^2} = \frac{2 \alpha b (N/P)}{[(N/P)^2 + b]^2} $$

Étant donné que toutes les variables sont des nombres non négatifs, l'expression

$$ \frac{b(N/P)}{[(N/P)^2 + b]^2} $$ représente elle aussi un nombre non négatif. Vue que $\alpha \in [0, \infty)$, l'expression

$$ \frac{dg}{dN} = 2\alpha \frac{b(N/P)}{[(N/P)^2 + b]^2} $$ se situe donc dans l'intervalle $[0, \infty)$.

### 11.11.5 Hassell-Varley Type I

Le modèle Hassell-Varley de Holling Type I est spécifié comme suit (Hassell et Varley 1969; Abrams et Ginzburg 2000) :

$$ f(N, P) = \alpha (N/P^m) $$

$$ g(N, P) = \alpha N P^{1-m} $$
Par conséquent, $\frac{dg}{dN} = \alpha P^{1-m}$, expression qui, étant donné $\alpha \in [0, \infty)$, se situe dans l'intervalle $[0, \infty)$.

### 11.11.6 Hassell-Varley Type II

Le modèle Hassell-Varley de Holling Type II est spécifié comme suit (Hassell et Varley 1969; Abrams et Ginzburg 2000) :

$$f(N, P) = \frac{\alpha(N/P^m)}{(N/P^m) + b}$$

$$g(N, P) = \frac{\alpha NP^{1-m}}{(N/P^m) + b}$$

Par conséquent :

$$\lim_{P \to 0} (f) = \lim_{P \to 0} \left[ \frac{\alpha N}{N + bP^m} \right] = \begin{cases} \frac{\alpha N}{N + b}, & m = 0 \\ \alpha, & m > 0 \end{cases}$$

$$\lim_{P \to \infty} (f) = \lim_{P \to \infty} \left[ \frac{\alpha N}{N + bP^m} \right] = \begin{cases} \frac{\alpha N}{N + b}, & m = 0 \\ 0, & m > 0 \end{cases}$$

$$\lim_{N \to 0} (f) = \frac{0}{0 + b} = 0$$

$$\lim_{N \to \infty} (f) = \lim_{N \to \infty} \left[ \frac{\alpha 1/P^m}{1/P^m + b/N} \right] = \frac{\alpha 1/P^m}{1/P^m} = \alpha$$

$$\lim_{P \to \infty} (g) = \lim_{P \to \infty} \left[ \frac{\alpha NP^{1-m}}{(N/P^m) + b} \right] = \begin{cases} \frac{\alpha N * \infty}{N + b} = \infty, & m = 0 \\ \frac{\alpha N * \infty}{N/\infty + b} = \infty, & 0 < m < 1 = \begin{cases} \infty, & 0 \le m < 1 \\ \alpha N/b, & m = 1 \end{cases} \\ \frac{\alpha N}{N/\infty + b} = \frac{\alpha N}{b}, & m = 1 \end{cases}$$

### 11.11.7 Hassell-Varley Type III

Le modèle Hassell-Varley de Holling Type III est spécifié comme suit (Hassell et Varley 1969; Abrams et Ginzburg 2000) :
$$ 
\begin{aligned}
f(N, P) &= \frac{\alpha (N/P^m)^2}{(N/P^m)^2 + b} \\
g(N, P) &= \frac{\alpha (N^2/P^{2m-1})}{(N/P^m)^2 + b}
\end{aligned}
$$

Par conséquent,

$$ \lim_{P \to 0} (f) = \lim_{P \to 0} \left[ \frac{\alpha N^2}{N^2 + b * P^{2m}} \right] = \begin{cases} \frac{\alpha N^2}{N^2 + b}, & m=0 \\ \alpha, & 0 < m \le 1 \end{cases} $$

$$ \lim_{P \to \infty} (f) = \lim_{P \to \infty} \left[ \frac{\alpha N^2}{N^2 + b + P^{2m}} \right] = \begin{cases} \frac{\alpha N^2}{N^2 + b}, & m=0 \\ 0, & 0 < m \le 1 \end{cases} $$

$$ \lim_{P \to \infty} (g) = \lim_{P \to \infty} \left[ \frac{\alpha N^2}{(N^2/P) + b * P^{2m-1}} \right] = \begin{cases} \frac{\alpha N^2}{N^2/\infty + b * 0} = \infty, & 0 \le m < 0.5 \\ \frac{\alpha N^2}{N^2/\infty + b * 1} = \frac{\alpha N^2}{b}, & m = 0.5 \\ \frac{\alpha N^2}{N^2/\infty + b * \infty} = 0, & 0.5 < m \le 1 \end{cases} $$

### 11.11.8 Beddington-DeAngelis

Le modèle de Beddington-DeAngelis est spécifié comme suit (Beddington 1975; DeAngelis et al. 1975) :

$$ f(N, P) = \frac{\alpha N}{N + cP + b} $$

$$ g(N, P) = \frac{\alpha NP}{N + cP + b} $$

Par conséquent,

$$ \lim_{N \to 0} (f) = \frac{0}{cP + b} = 0 $$

$$ \lim_{N \to \infty} (f) = \lim_{N \to \infty} \left[ \frac{\alpha}{1 + (cP + b)/N} \right] = \frac{\alpha}{1} = \alpha $$
$$ \lim_{P \to 0} (g) = \frac{0}{N+0+b} = 0 $$

$$ \lim_{P \to \infty} (g) = \lim_{P \to \infty} \left[ \frac{\alpha N}{c + (N+b)/P} \right] = \frac{\alpha N}{c} $$

$$ \lim_{N \to 0} (g) = \frac{0}{cP+b} = 0 $$

$$ \lim_{N \to \infty} (g) = \lim_{N \to \infty} \left[ \frac{\alpha P}{1 + (cP+b)/N} \right] = \frac{\alpha P}{1} = \alpha P $$

$$ \frac{dg}{dN} = \frac{\alpha P [N+cP+b] - \alpha NP}{[N+cP+b]^2} = \frac{\alpha [cP^2+bP]}{[N+cP+b]^2} $$

Étant donné que toutes les variables sont des nombres non négatifs, l'expression

$$ \frac{cP^2+bP}{[N+cP+b]^2} $$ représente elle aussi un nombre non négatif. Vue que $\alpha \in [0, \infty)$,

l'expression $\frac{dg}{dN} = \alpha \frac{cP^2+bP}{[N+cP+b]^2}$ se situe donc elle aussi dans l'intervalle $[0, \infty)$.

### 11.11.9 Crowley-Martin

Le modèle Crowley-Martin est spécifié comme suit (Crowley et Martin 1989) :

$$ f(N, P) = \frac{\alpha N}{(1+bN)(1+c[P-1])} $$

$$ g(N, P) = \frac{\alpha NP}{(1+bN)(1+c[P-1])} $$

Par conséquent,

$$ \lim_{N \to \infty} (f) = \lim_{N \to \infty} \left[ \frac{\alpha}{(1/N+b)(1+c[P-1])} \right] = \frac{\alpha}{b(1+c[P-1])} $$

$$ \lim_{P \to \infty} (g) = \lim_{P \to \infty} \left[ \frac{\alpha N}{(1+bN)(1/P+c[1-1/P])} \right] = \frac{\alpha N}{(1+bN)c} $$

$$ \lim_{N \to \infty} (g) = \lim_{N \to \infty} \left[ \frac{\alpha P}{(1/N+b)(1+c[P-1])} \right] = \frac{\alpha P}{b(1+c[P-1])} $$
$$
\frac{dg}{dN} = \frac{\alpha P}{1+c(P-1)} \left[ \frac{(1+bN)-bN}{(1+bN)^2} \right] = \frac{\alpha P}{[1+c(P-1)][1+bN]^2}
$$

Étant donné que toutes les variables sont des nombres non négatifs, l'expression

$\frac{\alpha P}{(1+bN)^2}$ représente elle aussi un nombre non négatif, tandis que $\frac{1}{[1+c(P-1)]}$ peut être positif ou négatif. Par conséquent, $\frac{dg}{dN}$ se situe dans l'intervalle $(-\infty, \infty)$.

11.11.10 Royama-Rogers

Le modèle Royama-Rogers est spécifié comme suit (Royama 1971; Rogers 1972) :

$$
\begin{align*}
f(N, P) &= N \left(1 - e^{\frac{-\alpha P}{1+\alpha T_h N}}\right) \\
g(N, P) &= NP \left(1 - e^{\frac{-\alpha P}{1+\alpha T_h N}}\right)
\end{align*}
$$

Par conséquent,

$$
\begin{align*}
\lim_{P \to 0} (f) &= N(1 - e^0) = 0 \\
\lim_{P \to \infty} (f) &= N(1 - e^{-\infty}) = N \\
\lim_{N \to 0} (f) &= 0 * (1 - e^{\frac{-\alpha P}{1+\alpha T_h N}}) = 0
\end{align*}
$$

La limite qui suit est de forme indéterminée (0 * $\infty$) :

$$
\lim_{N \to \infty} (f) = \lim_{N \to \infty} \frac{(1 - e^{\frac{-\alpha P}{1+\alpha T_h N}})}{(-1/N^2)}
$$

Par souci de simplicité, nous définissons les variables auxiliaires suivantes : $a = \alpha P$ et $b = \alpha T_h$ , de tel que

$$
(1 - e^{\frac{-a}{1+bN}})' = -e^{\frac{-a}{1+bN}} * \frac{ab}{(1+bN)^2}
$$

Par conséquent,
$$
\lim_{N \to \infty} (f) = \lim_{N \to \infty} \left[ e^{\frac{-a}{1+bN}} * \frac{ab}{(1+bN)^2} * N^2 \right] = \lim_{N \to \infty} \left[ e^{\frac{-a}{1+bN}} * \frac{ab}{1/N^2 + 2b/N + b^2} \right]
$$

$$
\lim_{N \to \infty} (f) = e^0 * \frac{ab}{0+0+b^2} = \frac{a}{b} = \frac{P}{T_h}
$$

Pour ce qui concerne la prédation totale :

$$
\lim_{P \to 0} (g) = 0 * N(1 - e^0) = 0 * N * 0 = 0
$$

$$
\lim_{P \to \infty} (g) = \infty * N(1 - e^{-\infty}) = \infty * N * 1 = \infty
$$

$$
\lim_{N \to 0} (g) = P * \lim_{N \to 0} (f) = 0
$$

$$
\lim_{N \to \infty} (g) = P * \lim_{N \to \infty} (f) = \frac{P^2}{T_h}
$$

$$
\frac{dg}{dN} = P\left(1 - e^{\frac{-\alpha P}{1 + \alpha T_h N}}\right) - PN\frac{\alpha^2 P T_h}{(1 + \alpha T_h N)^2} * e^{-\alpha P(1 + \alpha T_h N)}
$$

En définissant $v=e^{-a/(1+bN)}$ et selon les définitions de $a$ et de $b$ telles qu'établies ci-dessus,

$$
\begin{align*}
\frac{dg}{dN} &= P(1-v) - PN\left(\frac{ab}{(1+bN)^2}\right) & v &= P\left[1-v\left(1+\frac{Nab}{(1+bN)^2}\right)\right] \\
g'(N=0) &= P[1-e^{-a}] \in [0, \infty)
\end{align*}
$$

Étant donné que $\lim_{N \to \infty} \frac{dg}{dN} = 0$ et que $\frac{d^2 g}{dN^2}$ est négatif tant que $a, b, N \in [0, \infty)$, la dérivée $dg/dN$ doit elle aussi se situer dans l'intervalle $[0, \infty)$.

11.11.11 Kovai

L'équation de prédation Kovai est spécifiée comme suit :
$$ 
\begin{aligned}
f(N, P) &= \alpha (1 - e^{-Nu/(\alpha P)}) \\
g(N, P) &= \alpha P (1 - e^{-Nu/(\alpha P)}) \\
u &= 1 - \frac{b}{N} (1 - e^{-N/b})
\end{aligned}
 $$

Les limites du taux de prédation *f* sous des conditions extrêmes de *N* et de *P* se calculent comme suit :

$$ 
\begin{aligned}
\lim_{P \to 0} (f) &= \alpha [1 - \lim_{P \to 0} (e^{-Nu/(\alpha P)})] = \alpha [1 - 0] = \alpha \\
\lim_{P \to \infty} (f) &= \alpha [1 - \lim_{P \to \infty} (e^{-Nu/(\alpha P)})] = \alpha [1 - 1] = 0 \\
\lim_{N \to 0} (f) &= \alpha [1 - \lim_{N \to 0} (e^{[-N+b(1-e^{-N/b})]/(\alpha P)})] = \alpha [1 - e^0] = 0 \\
\lim_{N \to \infty} (f) &= \alpha [1 - \lim_{N \to \infty} (e^{-Nu/(\alpha P)})] = \alpha [1 - e^{-\infty}] = \alpha
\end{aligned}
 $$

Les limites du taux de prédation total sous des conditions extrêmes de *N* et de *P* se calculent comme suit :

$$ \lim_{P \to 0} (g) = \alpha * 0 * [1 - e^{-\infty}] = 0 $$

La limite de *g* tandis que la population des prédateurs approche l'infini se calcule comme suit :

$$ \lim_{P \to \infty} (g) = \lim_{P \to \infty} [\alpha P (1 - e^{-Nu/(\alpha P)})] $$

Cette limite étant de forme indéterminée $\infty*0$, nous reformulons l'expression et appliquons la règle de l'Hôpital :

$$ \lim_{P \to \infty} (g) = \lim_{P \to \infty} \frac{\alpha (1 - e^{-Nu/(\alpha P)})}{1/P} = \lim_{P \to \infty} \frac{\frac{Nu}{P^2} e^{-Nu/(\alpha P)}}{1/P^2} = \lim_{P \to \infty} [Nu * e^{-Nu/(\alpha P)}] = Nu * 1 = Nu $$

$$ \lim_{N \to \infty} (g) = \alpha P \lim_{N \to \infty} (1 - e^{-Nu/(\alpha P)}) = \alpha P (1 - 0) = \alpha P $$

La dérivée de *g* par rapport à *N* se calcule comme suit :
$$ \frac{dg}{dN} = \alpha P \frac{d}{dN} \left[ 1 - e^{-\frac{N+b(1-e^{-N/b})}{\alpha P}} \right] $$

Définissons $x=N-b(1-e^{-N/b})$ , ce qui donne :

$$ \frac{dg}{dN} = \alpha P \frac{d}{dN} \left[ 1 - e^{-x/(\alpha P)} \right] = \alpha P \frac{d}{dx} \left[ 1 - e^{-x/(\alpha P)} \right] \frac{dx}{dN} $$

$$ \frac{dx}{dN} = 1 - b \left( \frac{1}{b} e^{-N/b} \right) = 1 - e^{-N/b} $$

Par conséquent :

$$ \frac{dg}{dN} = \alpha P \left[ \frac{e^{-x/(\alpha P)}}{\alpha P} \right] \left[ 1 - e^{-N/b} \right] = \left[ 1 - e^{-N/b} \right] * e^{-\frac{N+b(1-e^{-N/b})}{\alpha P}} $$

Lorsque N et b se situent dans l'intervalle $[0, \infty)$ , l'expression $1-e^{-N/b}$ se situera dans l'intervalle $[0, 1)$ .

La deuxième partie de la dérivée est également bornée de manière similaire, tel que démontré ci-dessous :

$$ e^{-\frac{N+b(1-e^{-N/b})}{\alpha P}} = e^{-x/(\alpha P)} $$

En premier lieu, confirmons que $x$ se situe dans l'intervalle $[0, \infty)$ tant que $N, b \in [0, \infty)$ :

$$ x(N=0) = 0 - b(1 - e^0) = 0 $$

Étant donné que $\frac{dx}{dN} = 1 - e^{-N/b}$ se situe dans l'intervalle $[0, 1)$ , $x$ doit être un nombre non négatif tant que $N, b \in [0, \infty)$ . Par conséquent, $e^{-x/(\alpha P)}$ se situe également dans l'intervalle $[0, 1]$ et $\frac{dg}{dN}$ dans l'intervalle $[0, 1)$ .
## 11.12 Appendice B : Tableau supplémentaire

<table><caption>Tableau 11.3: CILA, soit « Critère d'information largement applicable » (en premier) et classement (en parenthèses) de chaque équation lorsqu'appliquée aux différents ensembles de données. Un classement ou une valeur de CILA plus proche de 0 indique un meilleure ajustement du modèle. Les noms des ensembles de données sont reproduits tel qu'ils apparaissent dans le code du matériel supplémentaire. KV : Kovai; DP : Dépendance à la proie; DR : Dépendance au ratio; HV : Hassell-Varley; BD : Beddington-DeAngelis; CM : Crowley-Martin; RR : Royama-Rogers. Les chiffres romains suivant les équations DP, DR et HV indiquent la forme de réponse fonctionnelle de Holling.</caption><thead><tr><th></th><th>KV</th><th>DPI</th><th>DPII</th><th>DPIII</th><th>DRI</th><th>DRII</th><th>DRIII</th><th>HVI</th><th>HVII</th><th>HVIII</th><th>BD</th><th>CM</th><th>RR</th></tr></thead><tbody><tr><td>Kratina et al, 2009, F</td><td>544 (5)</td><td>708 (10)</td><td>704 (9)</td><td>715 (11)</td><td>648 (7)</td><td>514 (4)</td><td>663 (8)</td><td>473 (1)</td><td>472 (0)</td><td>570 (6)</td><td>506 (3)</td><td>475 (2)</td><td>2408 (12)</td></tr><tr><td>Schenk et al, 2005, F</td><td>22,4 (2)</td><td>32,7 (8)</td><td>33,0 (10)</td><td>29,8 (7)</td><td>27,2 (3)</td><td>27,4 (4)</td><td>20,3 (0)</td><td>28,9 (5)</td><td>28,9 (6)</td><td>21,7 (1)</td><td>32,9 (9)</td><td>41,4 (11)</td><td>53,2 (12)</td></tr><tr><td>Long et al, 2012: Exp 2, F : Sand</td><td>53,2 (0)</td><td>145 (9)</td><td>69,1 (5)</td><td>1,56e6 (12)</td><td>177 (11)</td><td>54,1 (1)</td><td>71,0 (6)</td><td>155,1 (10)</td><td>56,4 (2)</td><td>71,7 (7)</td><td>56,6 (3)</td><td>66,9 (4)</td><td>123,3 (8)</td></tr><tr><td>Kfir 1983, F</td><td>142 (2)</td><td>768 (9)</td><td>170 (6)</td><td>273 (8)</td><td>1,31e3 (11)</td><td>143 (4)</td><td>183,5 (7)</td><td>815,1 (10)</td><td>128,6 (0)</td><td>143,1 (3)</td><td>135,1 (1)</td><td>148,9 (5)</td><td>1,35e3 (12)</td></tr><tr><td>Edwards 1961: Nasonia-Musca, F</td><td>231 (2)</td><td>394 (10)</td><td>318 (9)</td><td>317 (8)</td><td>400 (11)</td><td>232 (3)</td><td>228 (1)</td><td>295 (7)</td><td>233 (5)</td><td>227 (0)</td><td>232 (4)</td><td>245 (6)</td><td>907 (12)</td></tr><tr><td>Edwards 1961: Trichogramma-Sitotroga, F</td><td>255 (4)</td><td>959 (11)</td><td>787 (10)</td><td>768 (9)</td><td>532 (8)</td><td>227 (0)</td><td>299 (6)</td><td>412 (7)</td><td>228 (1)</td><td>274 (5)</td><td>239 (2)</td><td>243 (3)</td><td>3,13e3 (12)</td></tr><tr><td>Edwards 1961: Trichogramma-Sitotroga 48-54, F</td><td>68,8 (11)</td><td>55,0 (10)</td><td>53,5 (9)</td><td>52,2 (8)</td><td>45,9 (7)</td><td>42,2 (5)</td><td>40,8 (2)</td><td>41,4 (3)</td><td>40,3 (1)</td><td>39,0 (0)</td><td>42,1 (4)</td><td>42,9 (6)</td><td>136 (12)</td></tr><tr><td>Eveleigh & Chant</td><td>131</td><td>673</td><td>120</td><td>120</td><td>182</td><td>122</td><td>122</td><td>176</td><td>113</td><td>115</td><td>111</td><td>119</td><td>665</td></tr></tbody></table>
<table><thead><tr><th></th><th>KV</th><th>DPI</th><th>DPII</th><th>DPIII</th><th>DRI</th><th>DRII</th><th>DRIII</th><th>HVI</th><th>HVII</th><th>HVIII</th><th>BD</th><th>CM</th><th>RR</th></tr></thead><tbody><tr><td>1982: Phytoseiulus proto, F</td><td>(8)</td><td>(12)</td><td>(4)</td><td>(5)</td><td>(10)</td><td>(6)</td><td>(7)</td><td>(9)</td><td>(1)</td><td>(2)</td><td>(0)</td><td>(3)</td><td>(11)</td></tr><tr><td>Eveleigh & Chant 1982: Amblyseius proto, F</td><td>156 (4)</td><td>1,20e3 (11)</td><td>295 (10)</td><td>291 (9)</td><td>185 (8)</td><td>153 (3)</td><td>160 (5)</td><td>165 (7)</td><td>138 (1)</td><td>152 (2)</td><td>137 (0)</td><td>163 (6)</td><td>1,47e3 (12)</td></tr><tr><td>Eveleigh & Chant 1982: Phytoseiulus deuto, F</td><td>120 (5)</td><td>798 (11)</td><td>160 (9)</td><td>162 (10)</td><td>145 (8)</td><td>114 (2)</td><td>118 (4)</td><td>139 (7)</td><td>108 (1)</td><td>115 (3)</td><td>107 (0)</td><td>134 (6)</td><td>855 (12)</td></tr><tr><td>Eveleigh & Chant 1982: Amblyseius deuto, F</td><td>123 (2)</td><td>1,38e3 (11)</td><td>407 (9)</td><td>408 (10)</td><td>129 (3)</td><td>122 (0)</td><td>154 (8)</td><td>130 (4)</td><td>122 (1)</td><td>148 (7)</td><td>135 (6)</td><td>131 (5)</td><td>1,69e3 (12)</td></tr><tr><td>Mertz & Davies 1968, F</td><td>1,16e4 (12)</td><td>1,45e3 (1)</td><td>1,57e3 (4)</td><td>2,71e3 (10)</td><td>2,60e3 (8)</td><td>2,67e3 (9)</td><td>1,94e3 (7)</td><td>1,51e3 (2)</td><td>1,61e3 (5)</td><td>579 (0)</td><td>1,68e3 (6)</td><td>1,56e3 (3)</td><td>3,12e3 (11)</td></tr><tr><td>Taylor 1988, F: Medium A, kühniella</td><td>79 (4)</td><td>164 (11)</td><td>110 (8)</td><td>136 (10)</td><td>114 (9)</td><td>72 (2)</td><td>91 (5)</td><td>99 (7)</td><td>71 (1)</td><td>93 (6)</td><td>70 (0)</td><td>75 (3)</td><td>325 (12)</td></tr><tr><td>Taylor 1988, F: Small A, Kühniella</td><td>53 (0)</td><td>102 (11)</td><td>66 (7)</td><td>76 (9)</td><td>83 (10)</td><td>57 (6)</td><td>56 (4)</td><td>75 (8)</td><td>56 (5)</td><td>55 (3)</td><td>54 (1)</td><td>55 (2)</td><td>159 (12)</td></tr><tr><td>Taylor 1988, F: Medium P, interpunctella</td><td>71 (3)</td><td>217 (11)</td><td>157 (9)</td><td>188 (10)</td><td>102 (7)</td><td>76 (4)</td><td>107 (8)</td><td>79 (5)</td><td>68 (2)</td><td>98 (6)</td><td>67 (1)</td><td>67 (0)</td><td>520 (12)</td></tr><tr><td>Taylor 1988, F:</td><td>69</td><td>193</td><td>131</td><td>167</td><td>119</td><td>76</td><td>93</td><td>91</td><td>68</td><td>87</td><td>66</td><td>74</td><td>441</td></tr></tbody></table>
<table><thead><tr><th></th><th>KV</th><th>DPI</th><th>DPII</th><th>DPIII</th><th>DRI</th><th>DRII</th><th>DRIII</th><th>HVI</th><th>HVII</th><th>HVIII</th><th>BD</th><th>CM</th><th>RR</th></tr></thead><tbody><tr><td>Small P, interpunctella</td><td>(2)</td><td>(11)</td><td>(9)</td><td>(10)</td><td>(8)</td><td>(4)</td><td>(7)</td><td>(6)</td><td>(1)</td><td>(5)</td><td>(0)</td><td>(3)</td><td>(12)</td></tr><tr><td>Uttley 1980: Fig 4,3,1, F</td><td>120<br>(0)</td><td>237<br>(11)</td><td>134<br>(7)</td><td>136<br>(8)</td><td>171<br>(10)</td><td>121<br>(3)</td><td>124<br>(5)</td><td>159<br>(9)</td><td>121<br>(2)</td><td>121<br>(4)</td><td>120<br>(1)</td><td>126<br>(6)</td><td>317<br>(12)</td></tr><tr><td>Uttley 1980: Fig 8,3,2, F</td><td>114<br>(3)</td><td>168<br>(9)</td><td>114<br>(4)</td><td>118<br>(6)</td><td>181<br>(10)</td><td>118<br>(7)</td><td>115<br>(5)</td><td>2,27e6<br>(12)</td><td>112<br>(0)</td><td>113<br>(2)</td><td>112<br>(1)</td><td>161<br>(8)</td><td>211<br>(11)</td></tr></tbody></table>
# 12Chapitre 12 : Modelling predation: Theoretical criteria and empirical evaluation of functional form equations for predator-prey systems

## 12.1 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.

## 12.2 Introduction

Correctly representing the relationship between predators and prey is a key challenge to all attempts at modelling trophic interactions and population dynamics (Kratina 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 and 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.
Considering 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 and Varley 1969) and the Beddington-DeAngelis equation (aN/((N+cP+b)))
(Beddington 1975; DeAngelis et al. 1975). (Beddington 1975; DeAngelis et al. 1975). The
former represent a compromise between prey and ratio dependence, where the N/P
prey-predator ratio is replaced by N/P^m, where m is a parameter in the range [0,1] that
allows the function to transition smoothly between either prey or ratio 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=ax ), Type II (logarithmic-type growth; y=ax/(x+b) ), Type III
(sigmoidal type growth; y=ax²/(x²+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 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íznarová and 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. For a comprehensive in-depth review of functional response
equations, the reader is referred to (Jeschke et al. 2002).

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 and Ginzburg 2000). Ratio dependent equations also produce 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 and Ginzburg 2000), and studies have found experimental data partially between ratio and prey dependence (Schenk et al. 2005).

There is therefore a need for an equation that encompasses both ratio and prey dependence 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 ratio dependence in the Hassell-Varley forms is dependent on a parameter, not on prey or predator populations, while the theoretical considerations mentioned in the preceding paragraph suggest that population levels are important in determining when either relationship is more appropriate. Switching between equations based on the nature of predator interference in the system has also been proposed (Skalski and Gilliam 2001).

There is therefore a need for clear criteria and expected behaviours for all useful predation equations, especially in the context of discrete time models, 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 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 and Ginzburg 2000; Jeschke et al.
2004) have been recently used as part of larger ecological and population dynamics models (Chesson and Kuang 2008; Sanchez et al. 2018; Sanders et al. 2018; Imbert et al. 2020).

These shortcomings occur because most of these equations are built upon assumptions (mainly non-extreme prey and predator population densities) that are unlikely to hold true under natural conditions when applied to discrete time steps (as they most often are in prediction modelling and in many field studies). 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 (Gleeson 1994; Abrams and Ginzburg 2000) 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 et al. (2009) have observed, “choosing appropriate functional responses is crucial for adequate predictions of food web dynamics. When ratio dependence is incorporated into predator–prey models their stability is usually enhanced,” highlighting the importance of improved 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 when applied to discrete time steps, 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.

## 12.3 Standard criteria for predation functions

Predation equations may be defined as the predation rate per predator per unit time, $f(N,P)$, or as the total predation rate per unit time, $g(N,P)=P* f(N,P)$ where N and P are the prey and predator densities, respectively. The following sections present a
proposed set of standard criteria for predation equations when applied to discrete time step
simulations, 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 Hassell-
Varley equation (Hassell and Varley 1969), each combined with Type I, II and III Holling
functional response forms (Holling 1959; Abrams and Ginzburg 2000). Three other equations
which do not lend themselves to combination with distinct Holling functional response forms,
the Beddington-DeAngelis equation (Beddington 1975; DeAngelis et al. 1975), the Crowley-
Martin equation (Crowley and Martin 1989), and the Royama-Rogers equation (Royama
1971; Rogers 1972; Kfir 1983) were also tested. See Table 12.2 for a full definition of each
equation.

**12.3.1 Behaviour at high prey density**

As prey density increases, predation per predator should, at most, asymptotically approach
the maximum feeding rate for the predator. In mathematical terms, $\displaystyle 0 \le \lim_{N \to \infty} f(N, P) \le \alpha$,
and, by extension, $\displaystyle 0 \le \lim_{N \to \infty} g(N, P) \le \alpha P$.

**12.3.2 Behaviour at low prey density**

As prey density approaches 0, both total as well as per predator predation should tend
towards 0, that is, $\lim_{N \to 0} f(N,P) = \lim_{N \to 0} g(N,P) = 0$.

**12.3.3 Behaviour at high predator density**

As predator density increases, overall predation should increase asymptotically towards
the total prey population available, or, potentially, to a lesser value. The latter could be caused
by the presence of some habitat niches that offer complete predation protection to a limited
number of prey, or else due to interference between predator individuals. At the same time,
prey availability per predator should decrease as predators begin competing for prey. In
mathematical terms, $\lim_{P \to \infty} g(N,P) \le N$, though $\lim_{P \to \infty} f(N,P) = 0$.
### 12.3.4 Behaviour at low predator density

On the other hand, as predator density decreases, the total predation should also approach $0: \lim_{P \to 0} g(N, P)=0$. While predation per predator may vary, in no case should it be higher than the maximum feeding rate per predator; that is, $0 \le \lim_{P \to 0} f(N, P) \le \alpha$.

### 12.3.5 Other considerations

Finally, many equations consider the functional response of predation rates to be dependent of 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 \le \frac{dg}{dN} \le 1.$$

### 12.3.6 The Kovai equation

As can be later seen in Table 12.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.

$$f(N, P) = \alpha (1 - e^{-Nu/(\alpha P)}) \\ u = 1 - \frac{b}{N} (1 - e^{-Nb})$$
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, 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}{dN}uN=1-e^{-N/b}$, which is in the range [0, 1) for non-negative b, thereby ensuring that the effective prey availability will never increase faster than the prey density itself does. b may be interpreted as the maximum capacity of prey refuges

$$ \lim_{N \to \infty} (N - Nu) = \lim_{N \to \infty} [N - N + b(1 - e^{-N/b})] = 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 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 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.

## 12.4 Methods (empirical study)

12 experimental studies with predation (functional response) data were taken from the literature, for a total of 29 distinct datasets (Table 12.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 in Python 3.7 (Python Software Foundation (PSF) 2019), using the libraries Numpy (van der Walt et al. 2011), Pandas (McKinney 2010) and PyMC3 (Salvatier et al. 2016 p. 3) for mathematical analyses; graphs were plotted using Seaborn (Michael Waskom et al. 2018) and Matplotlib (Hunter 2007). All data and code used to run the analyses and plot the figures present in this article are included in the supplementary material.
<table><caption>Table 12.1: Studies and datasets used in the analyses. PD indicates prey dependence, RD ratio dependence, HV Hassell-Varley, CM Crowley-Martin, RR Royama-Rogers, 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.</caption><thead><tr><th>Reference</th><th>No. datasets</th><th>System studied</th><th>Equations used</th></tr></thead><tbody><tr><td>Kratina et al. (2009)</td><td>1</td><td>Benthic flatworm-paramecium</td><td>BD*, HV III, RD III</td></tr><tr><td>Fernandez-Maldonado et al. (2017)</td><td>1</td><td>Cannibalism in Nabis pseudoferus</td><td>O I</td></tr><tr><td>Líznarová and Pekár (2013)</td><td>5</td><td>Spider-ant (various species)</td><td>O IV</td></tr><tr><td>Schenk et al. (2005)</td><td>1</td><td>Paper wasp-shield beetle</td><td>HV III</td></tr><tr><td>Schenk and Bacher (2002)</td><td>2</td><td>Paper wasp-shield beetle</td><td>O III</td></tr><tr><td>Long et al. (2012)</td><td>4</td><td>Cannibalism in red king crab</td><td>PD II</td></tr><tr><td>Kfir (1983)</td><td>1</td><td>Egg parasitoid of Phthorimaea operculella</td><td>RR</td></tr><tr><td>Edwards (1961)</td><td>3</td><td>Parasitism (various species)</td><td>None</td></tr><tr><td>Eveleigh and Chant (1982)</td><td>4</td><td>Predatory mites</td><td>HV</td></tr><tr><td>Mertz and Davies (1968)</td><td>1</td><td>Cannibalism in flour beetles</td><td>O</td></tr><tr><td>Taylor (1988)</td><td>4</td><td>Ectoparasitism (various prey)</td><td>RR **</td></tr><tr><td>Uttley (1980)***</td><td>2</td><td>Damselfly/back swimmer-cladoceran</td><td>CM, HV II***</td></tr></tbody></table>

* Slightly modified to give a Type III functional response.

** Slightly modified.

*** Data and best-fit equations from figure in Skalski and Gilliam (2001).

Each of the equations in Table 12.2 was calibrated using Bayesian inference (Hamiltonian Monte Carlo No U-Turn Sampler algorithm (Hoffman and 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 functional response was modelled with a Poisson distribution as follows:

$$y \sim \text{Poisson}(\mu * n)$$
where *y* is the total number of predation events observed (across all replicates), *n* is the number of replicates, and *μ* 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:

$$y \sim \text{Binomiale}(n=N*n, p=\min(\mu/N, 0,999))$$

with the same parameter definitions as above, and *N* being the initial number of prey. Note that *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 et al. 2017), which takes the number of parameters into account (and penalises equations with more parameters in order to avoid overfitting) and was recommended over the Akaike Information Criterion (AIC) for use in Bayesian analysis (Gelman et al. 2014); lower scores indicate better model fit.

# 12.5 Results and discussion

## 12.5.1 Theoretical study

As indicated in Table 12.2, none of the equations currently available in the literature fulfills the criteria proposed in Section 11.3, 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.

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., $\lim_{P \to \infty} (g) = \infty$) or else could exceed it under certain parametrisations of the equations (e.g., $\lim_{P \to \infty} (g) = \alpha N/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 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.
<table><caption>Table 12.2: Comparison of critical extreme condition properties of common predation equations. T indicates Holling response type (I, II or III). f(N,P) represents the predation rate per predator and g(N,P)=P*f(N,P) represents the total predation rate, where α is the theoretical maximum predation rate per predator, N and P are the prey and predator densities, respectively, Th is the handling time per prey and other variables are empirical constants. α, b, c ∈ [0, ∞); m ∈ [0, 1]. Shaded squares indicate undesirable behaviours as compared to the ideal equation. See Appendix A: Calculation of non-trivial limits and derivatives for details on the calculation of the results presented here.</caption><thead><tr><th>Equation</th><th>f(N,P)</th><th>lim(f)<sub>P→0</sub></th><th>lim(f)<sub>P→∞</sub></th><th>lim(f)<sub>N→0</sub></th><th>lim(f)<sub>N→∞</sub></th><th>lim(g)<sub>P→0</sub></th><th>lim(g)<sub>P→∞</sub></th><th>lim(g)<sub>N→0</sub></th><th>lim(g)<sub>N→∞</sub></th><th>dg/dN</th><th>Reference</th></tr></thead><tbody><tr><td>Ideal equation</td><td></td><td>[0, α]</td><td>0</td><td>0</td><td>[0, α]</td><td>0</td><td>[0, N]</td><td>0</td><td>[0, αP]</td><td>[0, 1]</td><td></td></tr><tr><td>T I</td><td>αN</td><td>αN</td><td>αN</td><td>0</td><td>∞</td><td>0</td><td>∞</td><td>0</td><td>∞</td><td>[0, ∞)</td><td>(Holling 1959; Abrams and Ginzburg 2000)</td></tr><tr><td rowspan="2">Prey dependence</td><td>T II</td><td>αN / (N+b)</td><td>αN / (N+b)</td><td>αN / (N+b)</td><td>0</td><td>α</td><td>0</td><td>∞</td><td>0</td><td>αP</td><td rowspan="2"></td></tr><tr><td>T III</td><td>αN<sup>2</sup> / (N<sup>2</sup>+b)</td><td>αN<sup>2</sup> / (N<sup>2</sup>+b)</td><td>αN<sup>2</sup> / (N<sup>2</sup>+b)</td><td>0</td><td>α</td><td>0</td><td>∞</td><td>0</td><td>αP</td></tr><tr><td>T I</td><td>αN/P</td><td>∞</td><td>0</td><td>0</td><td>∞</td><td>αN</td><td>αN</td><td>0</td><td>∞</td><td>[0, ∞)</td><td></td></tr><tr><td rowspan="2">Ratio dependence</td><td>T II</td><td>αN/P / (N/P+b)</td><td>α</td><td>0</td><td>0</td><td>α</td><td>0</td><td>αN / b</td><td>0</td><td>αP</td><td rowspan="2"></td></tr><tr><td>T III</td><td>α(N/P)<sup>2</sup> / ((N/P)<sup>2</sup>+b)</td><td>α</td><td>0</td><td>0</td><td>α</td><td>0</td><td>0</td><td>0</td><td>αP</td></tr><tr><td>Beddington-DeAngelis</td><td>αN / (N+cP+b)</td><td>αN / (N+b)</td><td>0</td><td>0</td><td>α</td><td>0</td><td>αN / c</td><td>0</td><td>αP</td><td>[0, ∞)</td><td>(Beddington 1975; DeAngelis et al. 1975)</td></tr><tr><td>Crowley-Martin</td><td>αN / (1+bN)(1+c[P-1])</td><td>αN / (1+bN)*(1-c)</td><td>0</td><td>0</td><td>αP / (b+bc(P-1))</td><td>0</td><td>αN / (1+bN)c</td><td>0</td><td>αP / (b+bc(P-1))</td><td>(-∞, ∞)</td><td>(Crowley and Martin 1989)</td></tr><tr><td>Royama-Rogers</td><td>N(1-e<sup>-αP / (1+αT<sub>n</sub>N)</sup>)</td><td>0</td><td>N</td><td>0</td><td>P / T<sub>h</sub></td><td>0</td><td>∞</td><td>0</td><td>P<sup>2</sup> / T<sub>h</sub></td><td>[0, ∞)</td><td>(Royama 1971; Rogers 1972; Kfir 1983)</td></tr></tbody></table>
<table><thead><tr><th>Equation</th><th>f(N,P)</th><th>lim(f)<sub>P→0</sub></th><th>lim(f)<sub>P→∞</sub></th><th>lim(f)<sub>N→0</sub></th><th>lim(f)<sub>N→∞</sub></th><th>lim(g)<sub>P→0</sub></th><th>lim(g)<sub>P→∞</sub></th><th>lim(g)<sub>N→0</sub></th><th>lim(g)<sub>N→∞</sub></th><th>dg/dN</th><th>Reference</th></tr></thead><tbody><tr><td>T I</td><td>α(N/P<sup>m</sup>)</td><td>αN, m=0<br>∞, m&gt;0</td><td>αN, m=0<br>0, m&gt;0</td><td>0</td><td>∞</td><td>0, 0≤m&lt;1<br>αN, m=1</td><td>∞, 0≤m&lt;1<br>N, m=1</td><td>0</td><td>∞</td><td>[0,∞)</td><td></td></tr><tr><td>Hassell-<br>Varley</td><td>T II<br>α(N/P<sup>m</sup>) / (N/P<sup>m</sup>) + b</td><td>αN / (N+b), m=0<br>α, m&gt;0</td><td>αN / (N+b), m=0<br>0, m&gt;0</td><td>0</td><td>α</td><td>0</td><td>∞, 0≤m&lt;1<br>αN / b, m=1</td><td>0</td><td>αP</td><td>[0,∞)</td><td>(Hassell and<br>Varley 1969)</td></tr><tr><td>T III</td><td>α(N/P<sup>m</sup>)<sup>2</sup> / (N/P<sup>m</sup>)<sup>2</sup> + b</td><td>αN<sup>2</sup> / (N<sup>2</sup>+b), m=0<br>α, m&gt;0</td><td>αN<sup>2</sup> / (N<sup>2</sup>+b), m=0<br>0, m&gt;0</td><td>0</td><td>α</td><td>0</td><td>∞, 0≤m&lt;0.5<br>αN / b, m=0.5<br>0, 0.5&lt;m≤1</td><td>0</td><td>αP</td><td>[0,∞)</td><td></td></tr><tr><td>Kovai</td><td>α(1 - e<sup>-Nu/(αP)</sup>)<br>u=1 - b/N(1 - e<sup>-N/b</sup>)</td><td>α</td><td>α</td><td>0</td><td>α</td><td>0</td><td>Nu∈[0, N]</td><td>0</td><td>αP</td><td>[0,1)</td><td></td></tr></tbody></table>
12.5.2 Empirical study

As can be seen in Table 12.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 12.1 shows how the Kovai
equation, depending on its parametrisation, can generate different Holling response type
forms for different datasets. In addition to its accurate representation of both Type II and Type
III functional responses with respect to prey density (though, as expected, not Type IV), the
Kovai equation can also represent the impact of predator competition for prey on predation
rates.

<figure><img src="image_48.png" /></figure>

Figure 12.1: Example of different functional responses possible with the Kovai equation. a) Type III Holling response form (Schenk and Bacher 2002, experiment year 2000); b) Type II Holling response form (Fernandez-Maldonado et al. 2017); c) suboptimal fit, as expected, against data with a Type IV functional response (Líznarová and Pekár 2013, Zodarion-Tetramorium data); d) example of ability of Kovai equation to partition prey between predators when predators begin interfering amongst themselves (Uttley 1980, figure 4.3.1).
N: Prey density; P: Predator density.
Figure 12.2 shows the results of different models applied to predation data from Kratina 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 12.2, a). On the other hand, intermediate form equations (Hassell-Varley Type III, Beddington-DeAngelis, and Kovai; see Figure 12.2 b, c and d, respectively), were successful in distinguishing the effects of prey and predator densities. (Prey-dependent models were overall unsuccessful at modelling the data, as these models cannot distinguish between different values of P at all, and are not shown in the figure.)

<figure><img src="image_49.png" /><figcaption>Figure 12.2: 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</figcaption></figure>
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 12.3 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’ WAIC fit (as compared to the fit of equations with different functional responses) when fitted to the data. Figure 12.3 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.

<figure><img src="image_50.png" /><figcaption>Figure 12.3: 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 (lower is better). 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.</figcaption></figure>

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 for application to discrete time step simulation. 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 are 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.

## 12.6 Conclusion

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.
## 12.7 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.

## 12.8 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'é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.

## 12.9 Data availability

All data and code used to run the analyses and plot the figures present in this article are included in the supplementary material.

## 12.10 References

- Abrams PA. 2015. Why ratio dependence is (still) a bad model of predation. Biological Reviews 90: 794–814.

- Abrams PA and Ginzburg LR. 2000. The nature of predation: Prey dependent, ratio dependent or neither? Trends in Ecology and Evolution 15: 337–41.

- Arditi R and Ginzburg LR. 1989. Coupling in predator-prey dynamics: Ratio-Dependence. Journal of Theoretical Biology 139: 311–26.

- Beddington JR. 1975. Mutual Interference Between Parasites or Predators and its Effect on Searching Efficiency. The Journal of Animal Ecology 44: 331.

- Berryman AA, Gutierrez AP, and Arditi R. 1995. Credible, Parsimonious and Useful Predator-Prey Models: A Reply to Abrams, Gleeson, and Sarnelle. Ecology 76: 1980–5.

- Chesson P and Kuang JJ. 2008. The interaction between predation and competition. Nature 456: 235–8.
Crowley PH and Martin EK. 1989. Functional Responses and Interference within and between Year Classes of a Dragonfly Population. Journal of the North American Benthological Society 8: 211-21.

DeAngelis DL, Goldstein RA, and O'Neill RV. 1975. A Model for Tropic Interaction. Ecology 56: 881-92.

Edwards RL. 1961. The Area of Discovery of Two Insect Parasites, Nasonia vitripennis (Walker) and Trichogramma evanescens Westwood, in an Artificial Environment. Can Entomol 93: 475-81.

Eveleigh ES and Chant DA. 1982. Experimental studies on acarine predator-prey interactions: the effects of predator density on prey consumption, predator searching efficiency, and the functional response to prey density (Acarina: Phytoseiidae). Can J Zool 60: 611-29.

Fernandez-Maldonado FJ, Gallego JR, Valencia A, et al. 2017. Cannibalism: Do risks of fighting and reprisal reduce predatory rates? Community Ecology 18: 87-96.

Gelman A, Hwang J, and Vehtari A. 2014. Understanding predictive information criteria for Bayesian models. Stat Comput 24: 997-1016.

Gleeson SK. 1994. Density Dependence is Better Than Ratio Dependence. Ecology 75: 1834-5.

Hassell MP and Varley GC. 1969. New Inductive Population Model for Insect Parasites and its Bearing on Biological Control. Nature 223: 1133-7.

Hoffman MD and Gelman A. 2014. The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research 15: 1593-623.

Holling CS. 1959. The Components of Predation as Revealed by a Study of Small-Mammal Predation of the European Pine Sawfly. Can Entomol 91: 293-320.

Hunter JD. 2007. Matplotlib: A 2D Graphics Environment. Comput Sci Eng 9: 90-5.
Imbert C, Papaïx J, Husson L, et al. 2020. Estimating population dynamics parameters of cabbage pests in temperate mixed apple tree-cabbage plots compared to control vegetable plots. Crop Protection 129: 105037.

Jeschke JM, Kopp M, and Tollrian R. 2002. Predator Functional Responses: Discriminating Between Handling and Digesting Prey. Ecological Monographs 72: 95–112.

Jeschke JM, Kopp M, and Tollrian R. 2004. Consumer-food systems: Why type I functional responses are exclusive to filter feeders. Biological Reviews of the Cambridge Philosophical Society 79: 337–49.

Kfir R. 1983. Functional response to host density by the egg parasite Trichogramma pretiosum. Entomophaga 28: 345–53.

Kratina P, Vos M, Bateman A, and Anholt BR. 2009. Functional responses modified by predator density. Oecologia 159: 425–33.

Líznarová E et Pekár S. 2013. Dangerous prey is associated with a type 4 functional response in spiders. Animal Behaviour 85: 1183–90.

Long WC, Popp J, Swiney KM, and Van Sant SB. 2012. Cannibalism in red king crab, Paralithodes camtschaticus (Tilesius, 1815): Effects of habitat type and predator density on predator functional response. Journal of Experimental Marine Biology and Ecology 422–423: 101–6.

McKinney W. 2010. Data Structures for Statistical Computing in Python. In: Proceedings of the 9th Python in Science Conference.

Mertz DB and Davies RB. 1968. Cannibalism of the Pupal Stage by Adult Flour Beetles: An Experiment and a Stochastic Model. Biometrics 24: 247.

Michael Waskom, Olga Botvinnik, Drew O’Kane, et al. 2018. seaborn. Zenodo.

Python Software Foundation (PSF). 2019. Python.

Rogers D. 1972. Random Search and Insect Population Models. Journal of Animal Ecology 41: 369–83.
Rohatgi A. 2019. WebPlotDigitizer.

Royama T. 1971. A comparative study of models for predation and parasitism. Res Popul Ecol 13: 1-91.

Salvatier J, Wiecki TV, and Fonnesbeck C. 2016. Probabilistic programming in Python using PyMC3. PeerJ Computer Science 2: e55.

Sanchez C, Gamez M, Burguillo FJ, et al. 2018. Comparison of predator-parasitoid-prey interaction models for different host plant qualities. Community Ecology 19: 125-32.

Sanders D, Thébault E, Kehoe R, and Veen FJF van. 2018. Trophic redundancy reduces vulnerability to extinction cascades. PNAS 115: 2419-24.

Schenk D and Bacher S. 2002. Functional response of a generalist insect predator to one of its prey species in the field. Journal of Animal Ecology 71: 524-31.

Schenk D, Bersier L-F, and Bacher S. 2005. An experimental test of the nature of predation: neither prey- nor ratio-dependent: Experimental test of functional response. Journal of Animal Ecology 74: 86-91.

Skalski GT and Gilliam JF. 2001. Functional responses with predator interference: viable alternatives to the Holling Type II model. Ecology 82: 3083-92.

Taylor AD. 1988. Host Effects on Functional and Ovipositional Responses of Bracon hebetor. The Journal of Animal Ecology 57: 173.

Uttley M. 1980. A laboratory study of mutual interference between freshwater invertebrate predators.

Vehtari A, Gelman A, and Gabry J. 2017. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput 27: 1413-32.

Walt S van der, Colbert SC, and Varoquaux G. 2011. The NumPy Array: A Structure for Efficient Numerical Computation. Comput Sci Eng 13: 22-30.
Watanabe S. 2013. A Widely Applicable Bayesian Information Criterion. Journal of Machine
Learning Research 14: 867–97.

# 12.11 Appendix A: Calculation of non-trivial limits and derivatives

## 12.11.1 Prey dependent Type II

The prey-dependent Holling Type II model is specified as (Holling 1959; Abrams and
Ginzburg 2000):

$$f(N, P) = \frac{\alpha N}{N + b}$$

$$g(N, P) = \frac{\alpha PN}{N + b}$$

Therefore:

$$\frac{dg}{dN} = \frac{\alpha P(N+b) - \alpha PN}{(N+b)^2} = \frac{\alpha Pb}{(N+b)^2}$$

Since all variables are nonnegative numbers, $\frac{Pb}{(N+b)^2}$ is also a nonnegative number.

Given $\alpha \in [0, \infty)$, $\frac{dg}{dN} = \alpha \frac{Pb}{(N+b)^2}$ is therefore in $[0, \infty)$.

## 12.11.2 Prey dependent Type III

The prey-dependent Holling Type III model is specified as (Holling 1959; Abrams and
Ginzburg 2000):

$$f(N, P) = \frac{\alpha N^2}{N^2 + b}$$

$$g(N, P) = \frac{\alpha PN^2}{N^2 + b}$$

Therefore:

$$\frac{dg}{dN} = \frac{2 \alpha PN * (N^2 + b) - 2N(\alpha PN^2)}{(N^2 + b)^2} = \frac{2 \alpha PNb}{(N^2 + b)^2}$$
Since all variables are nonnegative numbers, $2 \frac{PNb}{(N^2+b)^2}$ is also a nonnegative number.

Given $\alpha \in [0, \infty)$, $\frac{dg}{dN} = \alpha \frac{2PNb}{(N^2+b)^2}$ is therefore in $[0, \infty)$.

### 12.11.3 Ratio dependent Type II

The ratio-dependent Holling Type II model is specified as (Arditi and Ginzburg 1989; Abrams and Ginzburg 2000):

$$f(N, P) = \frac{\alpha N/P}{N/P + b}$$

$$g(N, P) = \frac{\alpha N}{N/P + b}$$

Therefore:

$$\lim_{P \to 0} (f) = \lim_{P \to 0} \left[ \frac{\alpha N}{N + bP} \right] = \frac{\alpha N}{N} = \alpha$$

$$\lim_{P \to \infty} (f) = \lim_{P \to \infty} \left[ \frac{\alpha N}{N + bP} \right] = 0$$

$$\lim_{N \to 0} (f) = \frac{0}{0 + b} = 0$$

$$\lim_{N \to \infty} (f) = \lim_{N \to \infty} \left[ \frac{\alpha/P}{1/P + b/N} \right] = \frac{\alpha/P}{1/P} = \alpha$$

$$\lim_{P \to 0} (g) = \frac{\alpha N}{\infty} = 0$$

$$\lim_{P \to \infty} (g) = \frac{\alpha N}{0 + b} = \frac{\alpha N}{b}$$

$$\lim_{N \to 0} (g) = \frac{0}{0 + b} = 0$$

$$\lim_{N \to \infty} (g) = \lim_{N \to \infty} \left[ \frac{\alpha}{1/P + b/N} \right] = \frac{\alpha}{1/P} = \alpha P$$
$$dg/dN = \frac{\alpha[N/P+b] - \alpha N[1/P]}{[N/P+b]^2} = \frac{\alpha b}{[N/P+b]^2}$$

Since all variables are nonnegative numbers, $\frac{b}{[N/P+b]^2}$ is also a nonnegative number.

Given $\alpha \in [0, \infty)$, $\frac{dg}{dN} = \alpha \frac{b}{[N/P+b]^2}$ is therefore in $[0, \infty)$.

### 12.11.4 Ratio dependent Type III

The ratio-dependent Holling Type III model is specified as (Arditi and Ginzburg 1989; Abrams and Ginzburg 2000):

$$f(N, P) = \frac{\alpha (N/P)^2}{(N/P)^2 + b}$$

$$g(N, P) = \frac{\alpha N^2 P}{(N/P)^2 + b}$$

Therefore:

$$\lim_{P \to 0} (f) = \lim_{P \to 0} \left[ \frac{\alpha N^2}{N^2 + bP^2} \right] = \frac{\alpha N^2}{N^2} = \alpha$$

$$\lim_{P \to \infty} (f) = \lim_{P \to \infty} \left[ \frac{\alpha N^2}{N^2 + bP^2} \right] = \frac{\alpha N^2}{\infty} = 0$$

$$\lim_{N \to 0} (f) = \frac{0}{0+b} = 0$$

$$\lim_{N \to \infty} (f) = \lim_{N \to \infty} \left[ \frac{\alpha/P^2}{1/P^2 + b/N^2} \right] = \frac{\alpha/P^2}{1/P^2} = \alpha$$

$$\lim_{P \to 0} (g) = \lim_{P \to 0} \frac{\alpha N^2}{[(N^2/P) + bP]} = 0$$

$$\lim_{P \to \infty} (g) = \lim_{P \to \infty} \frac{\alpha N^2}{P[(N/P)^2 + b]} = 0$$

$$\lim_{N \to 0} (g) = \frac{0}{0+b} = 0$$