The effect of non-linear competitive interactions on quantifying niche and fitness differences

The niche and fitness differences of modern coexistence theory separate mechanisms into stabilizing and equalizing components. Although this decomposition can help us predict and understand species coexistence, the extent to which mechanistic inference is sensitive to the method used to partition niche and fitness differences remains unclear. We apply two alternative methods to assess niche and fitness differences to four well-known community models. We show that because standard methods based on linear approximations do not capture the full community dynamics, they can sometimes lead to incorrect predictions of coexistence and misleading interpretations of stabilizing and equalizing mechanisms. Specifically, they fail when both species occupy the same niche or in the presence of positive frequency dependence. Conversely, a more recently developed method to decompose niche and fitness differences, which accounts for the full non-linear dynamics of competition, consistently identifies the correct contribution of stabilizing and equalizing components. This approach further reveals that when the true complexity of the system is taken into account, essentially all mechanisms comprise both stabilizing and equalizing components and that local maxima and minima of stabilizing and equalizing mechanisms exist. Amidst growing interest in the role of non-additive and higher order interactions in regulating species coexistence, we propose that the effective decomposition of niche and fitness differences will become increasingly reliant on methods that account for the inherent non-linearity of community dynamics.


Introduction
Niche and fitness differences are two central concepts in coexistence theory which help us to understand species coexistence. Whether or not species coexist depends on how much species limit each other's growth compared to their own growth. From a mechanistic standpoint, various factors, such as resources (Letten et al. 2017), predators (Chesson and Kuang 2008;Holt et al. 1994) or mutualists (Johnson 2021), can limit population growth. Population sizes will in turn influence these factors (Meszena et al. 2006) (e.g., by depleting resources or by boosting the population size of their predators), as such creating the regulating feedback underpinning growth limitation. Thus, mechanistically, niche differences (1-niche overlap) measure how independent the feedback loops of two species are. If the feedback loops are completely independent, then niche differences are 1, conversely, if all feedback loops, intra-and interspecific, are equivalent, niche differences are 0. Fitness differences measure the relative strength of the feedback loops. Species will coexist when niche differences overcome fitness differences (Adler et al. 2007;Chesson 2013).
Unfortunately, computing niche and fitness differences taking into account the full complexity of limiting factors is not a straightforward task. Some available methods therefore apply to simple models only (Spaak et al. 2022b), where the details of population regulation are omitted. One often applied method is to reformulate the equilibrium conditions of a two-species community model in linear terms, equivalent to the equilibrium conditions of a linear Lotka-Volterra model (Letten et al. 2017;Godoy et al. 2014;Johnson 2021;Spaak et al. 2022b). This writes the equilibrium competitor's densities, which then act as limiting factors, in lieu of the actual factors that underpin species interactions. Often, this includes various assumptions on the time scales at which the limiting factors vary, as well as on the persistence of these factors (Letten et al. 2017;Johnson 2021;Chesson and Kuang 2008). While these simplifying assumptions will by design be violated in certain conditions, we do not know the implications of such violations for our capacity to understand coexistence. Such knowledge is important for application of coexistence theory, as is choosing the appropriate definition of niche and fitness difference (Godwin et al. 2020).
Here, we examine to what extent including non-linear characteristics of community models increases our understanding of coexistence. To this end, we assemble four community models, three of which explicitly model limiting factors, and one which considers direct non-linear effects of population density on per-capita growth. We rewrite the equilibrium conditions of all four models with linear functions and compare how the resulting niche and fitness differences differ from those of the full models including limiting factors and direct non-linear interactions. Based on this comparison, we first develop criteria to test whether the approximated models affect our ability to understand coexistence. Then, we graphically represent the assumptions made by the approximation method and show that these simplifications can result in large deviations in community dynamics, especially away from equilibrium. Importantly, we show that linear approximations of non-linear dynamics can sometimes result in the misunderstanding of community dynamics, e.g., positive niche differences even when both species occupy the same niche, and negative frequency dependence interpreted as positive frequency dependence and vice versa. Finally, we also show that accounting for the full non-linear dynamics of competition can yield new insights into the effects of mortality and resource availability on niche and fitness differences. Specifically, we find that changes in resource availability generally affect both niche and fitness differences, but additionally that local maxima and minima of niche and fitness differences exist.

Methods
Niche and fitness differences of four community models have so far been computed by analyzing the linear equilibrium conditions: Species competing for substitutable resources with a Holling type 1 response (Chesson 1990) and a Holling type 2 response (Letten et al. 2017), species competing for essential resources with a Holling type 2 response (Letten et al. 2017) and the annual plant model (Godoy and Levine 2014), see Table 1. For each of these models we used the niche and fitness differences ( N A and A A ) from the literature (Chesson 1990(Chesson , 2013(Chesson , 2018. The superscript A denotes their derivation from approximated models. To compute niche and fitness differences of the full model ( N C and F C ) we use the method outlined in Spaak and De Laender (2020). Spaak and De Laender (2020) have proposed this method based on theoretical considerations, which have since then been successfully used in phenomenological community models (Buche et al. 2022;Spaak et al. 2022b), but we have little understanding whether this method gives conceptual correct interpretation for a mechanistic community model (Spaak et al. 2022a). The superscript C denotes their derivation from the full, non-linearized (complex) model. We will use N and F to denote niche and fitness differences based on any mathematical definition. All simulations and computations were done in Python 3.7.1 with NumPy version 1.15.4.

Niche and fitness differences based on approximated models
The computation of the approximated niche and fitness differences N A and F A is based on rearranging the equilibrium condition of a community model, given by 1 into a form that is equivalent to the equilibrium condition of a Lotka-Volterra community model, given by: where N i is the density of species i, r i is the intrinsic growth rate, ii and ij are the intra-and interspecific competition coefficients. This involves writing the parameters r i , ij and ii of Eq. (1) as a function of the model parameters. Generally speaking, the linear approximations are chosen to match the zero net g rowth isoclines (ZNGI), i.e., Specifically, to compute the parameters of the approximation method one first sets r i = f i (0, 0) to be the intrinsic growth rate. Second, to match the growth rates of each species at their monoculture equilibrium, i.e., f i (N C i , 0) = 0 one sets ii = 1∕N C i , where N C i is the equilibrium density of species i in monoculture. Third, to match the growth rates of species i at the twospecies equilibrium, i.e., where N * i and N * j are the equilibrium densities of species i and j in the two-species community. By this procedure, the growth rates of both species are correct when no species are present, i.e., f i (0, 0) , at the (1) and at their own monoculture equilibrium, i.e., f i (N C i , 0). However, the ZNGI may not be linear, in which case such a linear representation of the ZNGI will necessarily be an approximation. Importantly, the invasion growth rate of species i, its growth rate when the resident species j is at its monoculture equilibrium density, i.e., f i (0,N C j ) , may not be correctly approximated, as we only fit the growth rate at its own monoculture equilibrium density. We show below, that the approximations do sometimes incorrectly approximate the equilibrium densities. Table 1 summarizes the linear representation for all these models. Given the interaction coefficients ij we can compute N A and F A as (Chesson 2018): Note that we slightly change the definition of F A i by defin- This change is purely aesthetic to ensure consistency between N A i and F C i , it does not affect the properties of N A . Specifically, this ensures that the approximate model and the full model, discussed below, have the same coexistence condition. Specifically, a species is assumed to coexist if

Niche and fitness differences based on the full model
We computed N C and F C for the full model based on a recently developed method (Spaak and De Laender 2020). Briefly, this method ensures that species with independent feedback loops have N C = 1 , while species with equivalent feedback loops have N C = 0 . Thus, the method critically depends on whether invasion analysis can predict coexistence for the model it is applied to (Spaak and De Laender 2020; Barabas et al. 2018;Chesson 1994;Pande et al. 2019). More precisely, for a model given by the per-capita growth rate Spaak and De Laender (2020) define N C and F C as Table 1 We consider four community models of which the full equations are given ("Model"), as well as their linear representation of the equilibrium dynamics as reported before (Chesson 2013;Godoy et al. 2014;Letten et al. 2017). The linearized versions of these models are reconstructed by replacing r i and ij in Eq. (1)by the expressions in the corresponding columns. The first three models are resource explicit models, in the main text we focus on two-species communities competing for two resources. The resource growth rates may be either in per-capita (as in the first model) or in total (as in the second and third models), this does not qualitatively change our results.
We chose resource growth rates as used in Chesson (2013) and Letten et al. (2017). For these models w il is the conversion of resource to biomass, u li is the utilization of resource l by species i, m i is the mortality rate, S l is the resource supply and k i respectively k il are halfsaturation constants. Consistently, subscripts i and l stand for species and resources, respectively. Subscript L i is the index of the more limiting resource of species i in monoculture. In the annual plant model g i is the germination rate, m i is the seed mortality rate (traditionally the model uses the survival rate s i = 1 − m i ), i is the net-production rate and ′ ij is the intra-or interspecific interaction Name Model where N C j is, as above, the monoculture equilibrium density of species j. A conversion factor, c j , which converts densities of species j to densities of species i, ensures that both species have the same dependence on limiting factors, i.e., equally strong feedback loops. c j is the solution of the equation ) to the hypothetical invasion growth rate when the two species had independent feedback loops ( f i (0, 0) ) and to another hypothetical invasion growth rate when the two species had identical feedback loops ( f i (c jN C j , 0) ). Two species coexist if the niche differences of both species are strong enough to overcome fitness differ-

Predicting and understanding coexistence
For N A and F A or, N C and F C , to predict coexistence they need to correctly predict the outcome of species interactions. For a two-species community there are three such outcomes (Tilman 1982). A first outcome is priority effects, where the outcome of competition depends on the starting conditions, this is only possible with positive frequency dependence, i.e., negative niche differences (Ke and Letten 2018). A second outcome is coexistence, where both species persist indefinitely, this is only possible with negative frequency dependence, i.e., positive niche differences. A third outcome is competitively exclusion, where the competitive dominant species excludes the other species. This is consistent with both negative or positive niche differences, and depends on the fitness differences ( 1−N C N C F C < −1 ). When N C and F C or, N A and A A predict these three outcomes for a given model, we consider these to correctly predict coexistence.
For N A and F A or, N C and F C to be used to understand coexistence, they need to reflect the mechanisms driving the outcome of species interactions. N measures the difference between the niches of two species. When the two species do not interact with each other, because they occupy completely different niches, then N should reflect total niche differences, i.e., N = 1 (Spaak and De Laender 2020; . Conversely, if two species occupy the same niche, then N should be 0. Testing whether two species occupy the same niche is difficult in general, as one may not know whether the species do differ in any niche axis not tested for. However, in the special cases where there is only one limiting factor, e.g., because species compete for only one resource or resources go extinct, species must occupy the same niche.
In addition, to understand species coexistence, N A and A A or, N C and F C , need to categorize a change in a model's parameter as a stabilizing and/or equalizing mechanism. That is, changes that increase niche differences are called stabilizing, while they are equalizing when they decrease fitness differences. For simplicity we will use the terms stabilizing and equalizing for both increases and decreases in niche differences or fitness differences, respectively. Importantly, mechanisms can be both stabilizing and equalizing. We will call any mechanism that only affects niche or fitness differences purely stabilizing or equalizing, respectively. We tested for all model parameters whether they are stabilizing or equalizing according to N A and F A or, N C and F C .

Results
Implicitly, N A and F A are linear approximations of the percapita growth rates, i.e., f i (N i (Fig. 1). Through this approximation we often lose information about community dynamics, most notably higher order interactions (Grilli et al. 2017;Levine et al. 2017), such as resource extinctions or explicit information about interactions with limiting resources (Letten and Stouffer 2019). Specifically, the linearization does not account for the possibility of substitutable resources going extinct (Chesson and Kuang 2008;Letten et al. 2017), or for the fact that the identity of the essential resource that is limiting to a species can change during competition, independent of changes in the resource supply rates (Letten et al. 2017). Importantly, the community model is approximated near the equilibrium of the system; however, it is evaluated away from the equilibrium at the invasion growth rates of the species. This difference between where the approximation is performed and where it is evaluated leads to a poor approximation of the growth rates, and of the invasion growth rate especially, for all four models (Fig. 1).
When species drive resources to extinction the approximation given in Table 1 may lead to an incorrect prediction of coexistence. That is two species may have F A > N A 1−N A but coexist nonetheless, or they may have F A < N A 1−N A but not coexist. However, in this case the approximation would also lead to different equilibrium densities, i.e., 1 − ∑ j ij N * j ≠ 0 , we therefore do not focus further on these cases. However, we caution that one must verify that the linear approximation indeed holds for the equilibrium condition (Appendices S1 and S2).

Accounting for complexity enhances understanding of coexistence
Both the approximation and the full niche and fitness differences predict the case where the two species do not interact and must have N = 1 in all four community models. For the other cases, N C and F C improve our understanding of coexistence.
For certain resource supply rates both species will drive one resource to extinction in monoculture ( Fig. 2 above blue dotted line or below green dotted line). In this case, the ZNGI of the species is non-linear and a linear representation will necessarily lose precision, irrespective of how exactly the linear representation of the model is chosen. When one resource is driven to extinction, both species effectively compete only for one resource and must have identical feedback loops and no niche differences. For all resource explicit models, N C can identify this case (i.e., N C = 0 ). However, N A ≠ 0 for these resource supply rates for the substitutional resource models ( Fig. 2A, B). Additionally, under competition for essential resources, the linearization method predicts N A = 0 for a community with priority effects (Appendix S2), which contradicts prior findings (Ke and Letten 2018;Mordecai 2011;Spaak and De Laender 2020). Both the approximated N A and A A Table 2 How N A and F A or N C and F C interpret coexistence across the four investigated models (columns). A ✓ indicates correct prediction or interpretation, an X indicates wrong prediction or interpretation. "No interaction": If the two species do not interact niche differences should be 1. "Same niche": If the two species occupy the same niche, e.g., in the case of competition for only one resource, then niche differences should be 0. "Stabilizing mechanism": A change in the model parameters can potentially affect both stabilizing and equalizing mechanisms. If a change in model parameters alters the outcome of competition from priority effects to coexistence (or vice versa), then the change in model parameters must affect stabilizing mechanisms. The annual plant model does not explicitly consider limiting factors, so there is no clear criteria for asserting that both species occupy the same niche, the two methods however agree perfectly, i.e., N C = 0 ⇔ N A = 0 , which we indicate with a " = " (see Appendix S5). Similarly, there is no clear criteria for asserting that a mechanism is stabilizing, and we found that the two definitions do not agree on this (indicated with " ≠")

Criteria
Substit-resources Holling type 1 Substit-resources Holling type 2 Essential-resources Holling type 2 Annual plant A: B: C: D: Fig. 1 The method to compute N A and F A implicitly approximates the actual community model f i (N i , N j ) with a linear per-capita growth rate r i (1 − ii N i − ij N j ) . e show this approximation explicitly as the scaled per-capita growth rates of the focal species (y-axis) as a function of the resident species density (x-axis) for the four full models (solid line = f i (0, N j )∕f i (0, 0) ) and their linear approximations (dashed = 1 − a ij N j ). The linear approximations do not capture the full complexity of the models and deviate substantially from the full model. Two main differences between the approximation and the full model stand out: First (A and B), the approximated equilibrium density of the resident N A j is smaller than the exact resident's equilibrium density N C j (indicated by different labels for N A j and N C j ). Second (A,B and C), the density for zero net growth is not equal for the approximation and the full model (intersection with gray line). These two differences occur because the linearization does not account for resources that go extinct (black triangles, A,B,C) or changes in which essential resource is limiting (black square, C). The full model eventually transforms into a horizontal line (A,B,C) when all resources are extinct and growth rates reduce to the density-independent mortality rates. However, species j would not naturally reach densities this high (i.e., above N C j ). These differences can lead to wrong predictions about the outcome of competition (see Table 2). Chosen parameters (A,B,C) are equivalent to Fig. 2 1 3 and the full N C and F C give exactly the same predictions for the absence of niche difference in the annual plant model, i.e., N C = 0 ⇔ N A = 0 , which we interpret as both methods correctly predicting if the two species occupy the same niche.
Any mechanism that is purely equalizing cannot alter the outcome of competition from coexistence to priority effects or vice versa (Ke and Letten 2018). We found that changes in mortality or resource supply rates can alter the outcome of competition from priority effects to coexistence in certain conditions, and therefore must be stabilizing (Appendix S1). N C and F C identify these mechanisms, as they affect N C i (Fig. 3). Conversely, the linear approximations classify any mechanism in the resource competition models that does not affect resource utilization u li or the conversion efficiency w il as purely equalizing; mechanisms that affect u li or w il are both equalizing and stabilizing (Table 2).
Alternatively, we can see that changes in resource supply rates can be stabilizing in the following scenario. Given two species, the resource supply rates determine whether both species are limited by the same resources (and consequently must have N = 0 ) or whether the species coexist (and consequently must have N > 0 ). Therefore, changes in resource supply rates can be stabilizing mechanisms. In general, N C and F C predict that essentially all mechanisms are both stabilizing and equalizing (Song et al. 2019).
Importantly, the new interpretation that changes in mortality or resources supply rates are stabilizing is not based on the specificities of N C and F C . Rather, we draw this argument based on intuitive knowledge of niche differences only. It should therefore apply to any reasonable definition of niche and fitness differences (Appendix S1).
The interpretation of stabilizing and equalizing mechanisms in the annual plant model depends on the method to assess niche and fitness differences. Using N A and A A , we would interpret that only mechanisms affecting ij can be 2 N (color) as a function of the resource supply rates for the three different resource explicit community models. N A and N C are not defined when one species cannot survive in monoculture (white regions). A,B: N A is independent of resource supply rates for substitutable resources. Moreover, the functional response to resource concentrations does not affect N A (given u li and w il ). In the region where both species are limited by only one resource (above blue dotted and below green dotted line) N A is non-zero, contrary to the assumption that species competing for one resource must have N = 0 . C: N A is only affected by changes in the resource supply rates if the limiting resource changes for a species (crossing a dotted line). This change from N A ≠ 0 to N A = 0 , however, is discontinuous. D,E,F: N C continuously depends on changes in resource supply rates and features local maxima (D-F) and minima (D,E). N C identifies cases where N C = 0 in the region where both species are limited by the same resource. Green and blue solid lines are the ZNGI, dashed lines delimit the coexistence region and dotted lines delimit the region where both species are limited by the same resource. Black dot represents the resource supply rates taken for Figs. 1 and 3 stabilizing, all other mechanisms are purely equalizing. Conversely, using N C and F C we would interpret that essentially all mechanisms are both stabilizing and equalizing, including changes in intrinsic growth rates, seed survival and seed mortality rates (Fig. 3). However, we do not know which of these two interpretations is more useful, as the arguments used for the resource explicit models do not apply for the annual plant model.

Accounting for complexity reveals new phenomena
A first new insight is the existence of local maxima and minima for N C . Importantly, the location of these local maxima and minima can be constructed geometrically, similar to how the coexistence region can be constructed (Appendix S3). The existence of these maxima and minima is unexpected, but not counter-intuitive. Imagine the scenario where we keep all community parameters fixed except S 1 , the supply rate of the first resource. For S 1 much smaller than S 2 , both species are limited by R 1 only, hence we must have zero niche differences, when S 1 and S 2 are comparable the species coexist, hence we expect positive niche differences, for much larger S 1 than S 2 the species are limited by R 2 , hence we expect again zero niche differences. Consequently, N should depend in a non-monotonic way on S 1 and will therefore exhibit a maximum. We should therefore expect that N depends in a non-monotonic way on both S l , as it indeed does. Another new insight from N C and F C is that competition for substitutable resources can lead to negative N C i , i.e., positive frequency dependence, for certain resource supply rates. However, even for these resource supply rates, the community will not exhibit priority effects. This highlights that negative niche differences (and thus positive frequency dependence) are necessary but not sufficient conditions for priority effects. If fitness differences are too strong, one species will exclude the other, regardless of the initial condition, which is the case in this scenario (Ke and Letten 2018).
Positive frequency dependence arises when species consume more of the resource that limits their competitor most (Ke and Letten 2018;Tilman 1982). For the resource supply rates chosen in Fig. 2D the blue species drives resource 2 to extinction in monoculture. The green species, as an invader, will therefore only consume resource 1 and have a horizontal utilization vector. The green species therefore Fig. 3 N C (color) as a function of mortality for the four investigated models. Mortality has a non-monotonic effect on N C and may lead to local maxima and minima of N C (panel B), similar to changes in the resource supply rates. Blank areas (Panel B) indicate too high mortality for species persistence. Parameters (except mortality) are chosen such as in Fig. 2, the black dot represents the same mortality as in Fig. 2

A B C D
consumes only the favored resource of its competitor (blue species), which leads to positive frequency dependence (Ke and Letten 2018). A condition for N C < 0 is therefore that resources go extinct, which is why N C < 0 does not occur for essential resources as species in monoculture cannot drive essential resources to extinction.

Discussion
An often used approach to compute niche and fitness differences is to represent the equilibrium dynamics of the community model with a linear model (Godwin et al. 2020;Godoy and Levine 2014;Letten et al. 2017;Chesson 2013;Johnson 2021). We have found that including non-linear species interactions improves our ability to understand coexistence, and yields three important new insights compared to linearization-based approaches. First, changes in mortality rates can affect both niche and fitness differences, contrary to earlier results (Chesson 2000;Godoy et al. 2014;Barabas et al. 2018;Petry et al. 2018;Letten et al. 2017). This confirms the earlier findings from Song et al. (2019) who found that niche and fitness differences are interdependent. Except in some special cases, changing any parameter always affected both niche and fitness differences. We found no mechanism that acts solely as equalizing or stabilizing.
Second, niche differences depend non-monotonically on the resource supply rates and environmental conditions. Again, this result is independent of the specific properties of N C and F C , but is based on intuitive properties of niche and fitness differences (see Appendix S1). We also found maxima and minima of niche differences as a function of resource supply. While we were able to show the existence and specific location of such extremes, we know little about their consequences. For example, niche differences are assumed to increase ecosystem function (Carroll et al. 2011;Striebel et al. 2009) or stability with respect to perturbation (Adler et al. 2007). The extent to which these supply rates translate to higher functioning and greater stability remains to be verified.
Third, niche differences based on non-linear community models highlight resource supply rates that result in negative N C (Fig. 2). At these resource supply rates, species have positive frequency dependence, which was unknown for this well-known community model. However, despite the positive frequency dependence and negative niche differences, the system is not driven by priority effects, as fitness differences are too strong (Ke and Letten 2018).
Niche and fitness differences should facilitate understanding of species coexistence by disentangling stabilizing from equalizing mechanisms. We have shown that essentially all mechanisms are both stabilizing and equalizing, as they can alter the outcome of competition from priority effects to coexistence. This is perhaps not surprising when the full non-linear dynamics are taken into account. Nevertheless, we would still anticipate that any given shift of traits or environment conditions will often be more strongly stabilizing or equalizing. It will be valuable in future work to identify the conditions under which equalizing and stabilizing mechanisms exhibit weaker or stronger interdependence. Decomposition methods that account for non-linear interactions will be critical to achieving this objective.

Future directions of niche and fitness differences
Recent research has shown that non-linear or higher order species interactions can affect community composition, community stability and coexistence (Grilli et al. 2017;Bairey et al. 2016;Mayfield and Stouffer 2017;Singh and Baruah 2019). We have shown that niche and fitness differences can be used to understand why non-linear and higher order species interactions, such as resource extinctions or changes in limiting resources, can affect community dynamics; however, our results remain limited to simple two-species communities. Applying niche and fitness differences, which include non-linear and higher order interactions, to more complex multi-species communities is an important next step to understand the effects of non-linear or higher order effects on population dynamics Letten and Stouffer 2019).
Modern coexistence theory and niche and fitness differences in particular should help us to predict and understand species coexistence. We have shown that two different methods to assess niche and fitness differences result in different interpretations of the underlying processes. However, there are a total of thirteen different methods to assess niche and fitness differences, with little consensus about which to use (Godwin et al. 2020;Spaak et al. 2022b). Some of these methods can account for non-linear species interactions, some of them cannot. Our results imply that these different methods can potentially lead to different interpretations of the underlying processes. For example, the methods of Carroll et al. (2011);Zhao et al. (2016) and Carmel et al. (2017) interpret essentially all processes as both equalizing and stabilizing, but to different extents (Appendix S4). Similarly, niche differences of Carroll et al. (2011) andZhao et al. (2016) depend non-monotonically on resource supply rates, but niche differences of Carmel et al. (2017) do not. Additionally, Song and Saavedra (2020) have shown that a mechanism increasing structural niche differences (Saavedra et al. 2017) may decrease niche differences as measured by the linearization approach. In short, our interpretation of stabilizing mechanisms depends on the method to assess niche differences (Spaak et al. 2022b).
Author contributions All authors conceived the study. J.W.S. and R.M. wrote the computer code. J.W.S. created the figure and performed the analysis. J.W.S. wrote the first draft. All authors discussed content and contributed to revisions.
Funding Open Access funding enabled and organized by Projekt DEAL. F.D.L. received support from grants of the University of Namur (FSR Impulsionnel 48454E1) and the Fund for Scientific Research, FNRS (PDR T.0048.16).

Data availability
The study used no empirical data.
Data accessibility All computer code used in this manuscript will be made publicly available on figshare.
Code availability All code used in the project will be made publicly available. The code is all written in Python, a publicly available openaccess programming language.

Consent to participate Not applicable
Consent for publication All authors agreed to be listed as authors.

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.