Modelling Holling type II functional response in deterministic and stochastic food chain models with mass conservation

The Rosenzweig-MacArthur predator-prey model is the building block in modeling food chain, food webs and ecosystems. There are a number of hidden assumptions involved in the derivation. For instance the prey population growth is logistic without predation but also with predation. In order to reveal these we will start with modelling a resource-predator-prey system in a closed spatially homogeneous environment. This allows us to keep track of the nutrient flow. With an instantaneous remineralisation of the products excreted in the environment by the populations and dead body mass there is conservation of mass. This allows for a model dimension reduction and yields the mass balance predator-prey model. When furthermore the searching and handling processes are much faster that the population changing rates, the trophic interaction is described by a Holling type II functional response, also assumed in the Rosenzweig-MacArthur model. The derivation uses an extended deterministic model with number of searching and handling predators as model variables where the ratio of the predator/prey body masses is used as a mechanistic time-scale parameter. This extended model is also used as a starting point for the derivation of a stochastic model. We will investigate the stochastic effects of random switching between searching and handling of the predators and predator dying. Prey growth by consumption of ambient resources is still deterministic and therefore the stochastic model is hybrid. The transient dynamics is studied by numerical Monte Carlo simulations and also the quasi-equilibrium distribution for the population quantities is calculated. The body mass of the prey individual is the scaling parameter in the stochastic model formulation. This allows for a quantification of the mean-field approximation criterion for the justification of replacement of the stochastic by a deterministic model.


Introduction
Predator-prey models such as the Rosenzweig-MacArthur model (RM-model) (Rosenzweig and MacArthur, 1963), are important building block in food chain, food web and ecosystem models.A trophic interaction is described by the Holling type II functional response.A descriptive logistic growth for the self limiting prey population dynamics is assumed.The dynamics of the prey-resources is not modelled explicitly introducing the bottleneck for having mass conservation (Kooijman et al., 2007;Kooi et al., 1998).
The starting point here is a model discussed in Kooi et al. (1998Kooi et al. ( , 2002) ) for the lowest three trophic levels of an aquatic food chain whereby resource for the prey in the closed environment is modelled explicitly.This allows to keep track of the resource flow and to fulfill mass conservation for the whole resource-prey-predator system.The volume of the environment is constant and consequently total number of individuals and spatial concentration are interchangeable quantities for the size of a population.
The trophic interaction between the prey and the resources is modelled by Lotka-Volterra functional response, where the caught prey biomass is immediately converted into predator biomass.The predator prey interaction where handling prey requires time is given by the Holling type II response function.Assuming instantaneous remineralisation of the products and dead body mass for both prey and predator populations, the conservation of mass law can be used to reduce the dimension of the model by an elimination of the nutrient variable.This gives a two-dimensional mass balance predator-prey model, called the PPh-model of which the dynamics will be compared with that of the also two-dimensional predator-prey RM-model.
Table 1 gives an overview of the used notations and abbreviations for the models.For an introduction into ecosystem-modelling in the closed systems with recycling of elemental matter and using mass-balanced models we refer also to Gurney and Nisbet (1998).
To model predation and consumption explicitly the predator population is split up in searchers and handlers.The predator-prey model becomes three dimensional, called the PPd-model.When the searching/ handling of prey is much faster than the population size changes a timescale separation argument, called the quasi-steady-state assumption (QSSA), is used to justify usage of the Holling type II functional response providing a reduction of the PPd-model back to the PPh-model.
We will show that this procedure leads to an inconsistency where one term appears in expressions for both fast and slow when also the dynamics of the prey population is taken into account.It will be shown that the ratio of the body size of the prey individuals and the predator individuals can be used the time-scale parameter.In that sense we revisit the derivation of the Holling type II functional response.
Empirical results show that this ratio is smallest for the lower trophic levels of a food chain, Jonsson and Ebenman (1998) and we will show that a small ratio justifies the use of the time-scale separation argument.Therefore, usage of an algebraic expression (Holling's disk equation Gurney and Nisbet, 1998) is allowed for these lower trophic levels where the predator dynamics is again described by one ODEfor the whole predator population with the Holling type II functional response.This is studied by comparison of the results of a bifurcation analysis of the RM-model with the two mass-balance model versions: the PPh-model and the PPd-model.
A general problem with the application of deterministic models is that biological phenomena in the real world have often dynamical behaviour that is intrinsically erratic due to random processes, and thus can be better described by stochastic systems.This random behaviour is due to bio-variability in the physiological processes at the individual level together with populations that consists of small numbers of the individuals.The formulation of the stochastic PPs-model is based on internal physiological random processes.The model variables are the same as in the deterministic PPd-model where the total system biomass is fully conserved across the random events.It is an alternative for models where an external 'noise' component is added to a deterministic ODE model.The resulting model equations are called SDE stochastic differential equations, see for instance (Gardiner, 2009;Stollenwerk et al., 2017) As in the PPd-model we model the trophic interaction by searching and handling processes and split the predator population in prey-searching and prey-handling predators.Searching is modelled using the law of mass action which gives the encounter rate of one prey and one predator and handling with exponentially distributed handling times.Furthermore the proportional death rate for the predator populations is modelled as a stochastic death process.All actions are modelled as Poisson processes modelling counting the completely at random occurrences of certain events that appear to happen at a certain intensity.
The growth of the prey population remains described by its deterministic formulation.This shows that the PPs-model is hybrid (Kang and Kurtz, 2013;Kang et al., 2014) involving differential equations driven stochastic equations, that is, a Markov process with piecewise deterministic ODEs with a time dependent solution.The motivation to model the searching, handling and the populations death process stochastically, while the prey growth process remains deterministically, is that in this study the size of the nutrient particles (considered as a nutrient soup) remains relatively small with respect to the prey individuals.
The analysis of the models will be used to address the following research topics 1.In the mass-balance formulation using biomass as descriptive variable for the size of the populations, an inconsistency occurs because the same term appears in expressions for both fast (searching and handling predator) and slow prey timescales.Is it possible to overcome this problem by using besides biomass also number of individuals as a dimension?2. Are the differences in the bifurcation diagrams and consequently the long-term dynamics, for the two-dimensional mechanistic massbalance PPh-model and the descriptive RM-model large? 3. Are the differences in the long-term dynamics of the deterministic and stochastic versions in application to food chain models, with parameter values within realistic ranges, of a tri-trophic aquatic food chain model large?
Different analysis methods, such as ODE solvers and bifurcation packages are used to study the dynamics of the deterministic models.We refer to Guckenheimer and Holmes (1985), Kuznetsov (2004) for an introduction into bifurcation analysis.In case of stochastic models, Monte Carlo type numerical simulations produce realizations from which mean and variances are calculated.In Kooijman et al. (2007) the direct Poissonian process method was used and here we use the direct Gillespie method (Gillespie, 1976;2001;Stollenwerk et al., 2017) based on exponential next-event distributions for law of mass-action predator-prey encounters and mortality.Furthermore, we study a quasi-equilibrium distribution approximated by a multivariate normal distribution (Grasman and Herwaarden, 1999).
In phase-space the stochastic solution shows the characteristic counter clock limit cycles where the predator dynamics follows that of the prey around the stable equilibrium of the associated deterministic model (Gurney and Nisbet, 1998).Stochastic simulations show that this holds for spiral equilibria but also for node equilibria.
The paper is organised as follows.In Section 2 the predator-preyresources PPR-model for the lowest three trophic levels of an aquatic food chain is formulated.In Section 3 that model is reduced by one dimension using mass conservation.The Holling type II functional response describes the predator-prey trophic interaction in this PPhmodel.In Section 4 we focus on the modelling of trophic interactions.The predator population is divided into searchers and handlers of prey, called the PPd-model and one extra parameter is introduced: the ratio of the body size of the prey versus the predator individuals.A small values of the ratio gives by using a time-scale argument the PPh-model again.Section 5 gives the observed values of this ratio for real food chains depending on the trophic level.In Section 6 the dynamics of the PPdmodel and the RM-model deterministic model is analysed.In Section 7 we investigate the stochastic version of the PPd-model of the PPs-model by random simulations.This gives again one extra scaling parameter

Three-trophic level population model
This section describes the predator-prey-resources PPR-system in a batch reactor environment which is spatially homogeneous and closed for biomass.Such a system forms for instance the three lowest trophic levels in a food chain, (Gurney and Nisbet, 1998, Chapter 7).In marine systems autotrophic producers, for instance phytoplankton (algae), convert light energy to chemical energy which is stored in its organic compounds from inorganic nutrients they take up from the environment.Their compounds are the resources that become available to heterotrophic consumers.
We denote the resource concentration by x 0 , prey biomass density by x 1 and predator biomass density by x 2 expressed for instance in units Cmoles per constant environmental volume.The following physiological processes in the populations are considered: assimilation, growth and maintenance govern how populations can exist and under which conditions they survive.The efficiencies of these processes are in general not 100% and therefore products are excreted by the prey and predator organisms in the reactor.
The seven dimensional model reads (see Kooi et al., 2002) dx 0 dt The prey consumes resource at rate I 01 f 01 , a proportion μ 01 f 01 of it goes to reproduction of new prey.The assimilation efficiency is y 01 = μ 01 /I 01 .The remaining proportion (I 01 − μ 01 f 01 ) cannot be processed by the prey, and hence, is excreted to form the excreted biomass p 1 .At rate α 1 , this excreted biomass becomes a new resource that can then be used by prey and predators.Another contribution to resource comes from dead prey z 1 .We assume that prey individuals die with m 1 , and at rate β 1 dead prey are converted into the resource.Similarly predator consumes prey at rate I 12 f 12 described by a scaled Holling type II functional response, a proportion μ 12 f 12 of it goes to reproduction of new predator.The assimilation efficiency is y 12 = μ 12 /I 12 .The remaining portion (I 01 − μ 01 )f 12 cannot be processed by the predator, and hence, is excreted to form the excreted biomass p 2 .At rate α 2 , this excreted biomass becomes ads to the resource as well as dead predator z 2 .We assume predators die with rate m 2 , and at rate β 2 dead predators are converted into the resource.
The consumption of the nutrients by the prey is modelled by the Lotka-Volterra or linear functional response, f 01 , and the predator-prey interaction is modelled by the Holling type II functional response given by the scaled functional responsef 12 , hence By assuming instantaneous remineralisation of all excreted products from the assimilation, maintenance processes of the populations and dead biomass, model ( 1) becomes the PPR-model: Let C(t) be the sum of the three biomass densities, then summation of the three Eq.( 3) gives that is, the sum of the biomasses is a given constant C 0 .This is the deterministic mass-balance population model where f 12 is the Holling type II functional response Eq. ( 2) and the state values are expressed in dimension of biomass per environmental volume.Table 2 gives a list of all the parameters and variables.The model is referred to as the three dimensional PPR-model.

Deterministic predator-prey models
In this section we discuss two predator-prey models with Holling type II functional response.Using mass the conservation principle we can eliminate the variable x 0 and system Eq.( 3) is equivalent to the twodimensional predator-prey system which we call the PPh-model where the Holling type II functional response in given in Eq. ( 2).Mass conservation gives the side-condition This system resembles the RM-model where the lowest trophic level, the nutrients for the prey, is not modelled explicitly.In that formulation the prey grows logistically in absence of the predator and this term is not altered when the predator population exists, that is the − x 1 x 2 term is missing, hence When the predator population is extinct, x 2 = 0, models Eqs.(5a) and (6a) for the remaining prey population are equal because the extra term is zero.Starting with Eqs.(3a), (3b) where x 2 = 0 and following the same mass conservation principle we get a dimension reduction for the prey dynamics to the logistic growth, see also Kooi et al. (1998) and (Gurney and Nisbet, 1998, Chapter 7).
However, we show here that starting with the three level predatorprey-resources system Eq.( 3), the same principle gives the extra − x 1 x 2 term in Eq. (5a).The biological interpretation is that the prey uses the ingested resources partly for its own but also an extra amount needed as building blocks by the predators.The mass conservation law shows that in a closed reactor system the prey and predator are also (indirectly) competing for the ambient resources in the reactor and therefore, there are two negative predator-prey interaction terms in the equation for the growth of the prey population.
Using the usual parameters and notations, the RM-model is presented in the familiar expressions with intrinsic growth rate r and carrying capacity K.
In order to have r > 0 and K > 0 we assume: μ 01 C 0 /k 01 − m 1 > 0. Hence, the total biomass in the closed system, parameter C 0 is an environmental parameter similar to the carrying capacity but with a different interpretation because it is a conserved quantity and does not refer to an asymptotic solution.While in the RM-model we have for a stable preyonly equilibrium lim t→∞ x 1 = K, in the mass-balance equation it also depends on lim t→∞ x 2 .For the biological interpretation of the parameters the reader is referred to Table 2.
It is well known that the Holling type II functional response f 12 , Eq. ( 2), can be derived using mechanistically motivated time-scale separation between predator-prey physiological processes of searching and handling, see for instance (Dawes and Souza, 2013;Metz and Diekmann, 1986;Metz and Batenburg, 1985) and (Stollenwerk et al., 2017, App. A2).Using a time-scale separation argument, Eq. ( 5) with Eq. ( 2) is purely formulated at the population level.Generally, in this derivation the predator population is split up into two sub-population, searchers and handlers yielding a three-dimensional system where the time-scale parameter ε is small, and the slow time scale t is such that τ = t /ε defines the fast timescale: Where searching and handlers predators change at a fast timescale.
With ε→0, the fast time-scale dynamics of the predator system (x 2s (τ), x 2h (τ)) is described by The quasi-steady-state variables x * 2s and x * 2h are the fast equilibrium values given by where we used x 2 = x 2s + x 2h and Eq. ( 2).Adding Eqs. ( 9b) and (9c) gives, however only in the unrealistic case ε = 1, the two-dimensional model Eq. ( 6) where the conversion from prey biomass into predator mass is described by an algebraic equation, the Holling type II functional response, Eq. ( 2).The last term in Eq. (9a) and the first term in Eq. (9b) are unequal when ε < 1 and therefore, since they present the same quantity, there is an inconsistency for small ε.In the next section we will revise this derivation.

Deterministic PPd-model
In compartmental models all individuals that constitute, for both prey and predator populations, are assumed identical (also their chemical composition).As a consequence the dimension of the population can be modelled either by numbers of individuals or by the total biomass.Let n 1 ∈ N denote the number of the prey individuals per environmental volume and n 2 ∈ N the number of predator individuals per environmental volume.
Let x ε 1 denote the biomass of an single prey individual and x ε 2 the biomass of a predator individual.Since all individuals are identical we have for the prey population x 1 = n 1 x ε 1 and for the predator population x 2 = n 2 x ε 2 .We will need the ratio of the masses of the individuals of the two populations, the coefficient y ε 12 = x ε 1 /x ε 2 .For the modelling of the predator-prey interaction the predator population is divided in searchers and handlers where the number of predator individuals are denoted by n 2s = x 2s /x ε 2 and n 2h = x 2h /x ε 2 and the number of prey individuals by n 1 = x 1 /x ε 1 and of all predators n 2 = x 2 /x ε 2 .When searching predators meet prey at contact rate b they are converted into handling predators.These handling predators do not feed and they convert back into searching predators with a constant specific rate k where 1/k is the mean handling time for a prey individual.The equations for the searching and handling subpopulation using the population individuals n 1 , n 2s and n 2h are derived using the pseudo-reaction scheme The equations read for the trophic interaction We interpret now the results in the mass density formulation, that is the state variables are the population biomass densities x 1 , x 2s and x 2h .The implementation of these details requires a partitioning of the predators x 2 (t) = x 2s (t) + x 2h (t) into searching x 2s and food-handling x 2h sub-populations.
Using biomasses again instead of numbers for the state variables we get We assume that the interactions of the population individuals n 1 , n 2s and n 2h , that rate b is related to I 12 /k 12 and k to I 12 as follows Both rates b and k are inversely proportional to the ratio y ε 12 of the biomass of the prey and predator.This is motivated as follows where we assume that the predator body size is fixed.When the biomass of the prey individuals is, for instance, halved, the number of prey individuals doubles and according to the law of mass action also the attack rate b doubles.Hence, we assume the attack rate proportional to the preypredator size ratio y ε 12 and this is equivalent to the assumption made in Troost et al. (2008) that the capture time is proportional to 1/y ε 12 .Furthermore, when the size of the prey individuals is halved (and the digestion rate in the predator remains the same) also the handling time is halved and therefore the handling rate k per individual doubles.
In addition we introduce the fast time scale variable τ = t/y ε 12 besides the slow time scale t.With y ε 12 →0, the fast time-scale dynamics of the predator system (x 2s (τ), x 2h (τ)) is described by Eq. ( 10) with ε = y ε 12 .The positive fast equilibrium, quasi-steady-state (x * 2s , x * 2h ) Eq. ( 11), of the system Eq.(10) needs to be unique and stable considered as the equilibrium of the single equation: The Jacobian evaluated in the quasi-steady-state equilibrium reads and hence the quasi-steady-state equilibrium Eq. ( 11) is globally stable.
N. Stollenwerk et al.This shows that there are no spurious equilibria occurring and therefore that the application of the time-scale technique is robust.
Combining the predator-prey trophic interaction Eq. ( 10) with ε = y ε 12 the terms presenting the other processes such a prey growth and prey and predator mortality from Eq. ( 5) with x 2 (t) = x 2s (t) + x 2h (t), gives the PPd-model The three-dimensional system PPd-model ( 18) becomes equivalent with the PPh-model when limy ε 12 ↓0.To show this we add of (18b) and (18c) and we use The PPd-model is also taken as starting point for the formulation the stochastic version PPs-model in Section 7.

Fig. 1.
Phase-space analysis for prey x 1 (t) and predator x 2 (t) phase-space analysis for left column k 12 = 1/3 (stable equilibrium: '•') and right column k 12 = 1 /8 (unstable equilibrium: '∘').The initial values for the state variables prey x 1 (0) = 0.2 and predator x 2 (0) = 0.7 are indicated by '□'.The blue curve is the parabolic x 1 -nullcline.Top panel: RM-model system (6).Middle panel: system PPh-model: Eq. ( 5).Bottom panel: PPd-model system Eq.( 18).The additional parameter is y ε 12 = 1/4.For comparison the equilibrium point in the left panel indicated by the open circle ∘ on the nullcline belongs to the PPh-model is shown.Note that there is no convergence to this PPh-model equilibrium '∘' but to a somewhat higher prey value x * 1 .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) N. Stollenwerk et al.

Body mass ratios in food chains
In Jonsson and Ebenman (1998) the empirical pattern of decreasing predator-prey body mass ratios with increasing trophic height is studied.The predator-prey body mass ratio 1/y ε 12 was proposed to be the following function of the trophic position i of the consumer in the chain, and decreases toward unity with increasing consumer trophic position as observed in many real food chains The starting value for the average prey-predator body mass ratio between the primary and secondary consumer (trophic level, i = 1) in a data set analyzed by the authors is y ε 12 ≈ 1/300, see (Jonsson and Ebenman, 1998, Eq. ( 5)).For i = 2 y ε 12 is 1/110, i = 2: 1/42, i = 3: 1/16, i = 4: 1/6.5.It increases toward unity with increasing predator trophic level.So, for higher levels of long food chains, the sizes of the predator and its prey become almost equal.
Note that in Jonsson and Ebenman (1998) an additional parameter k occurs but the authors state that its value is close to 1 and therefore it is not introduced here.

Phase-space and bifurcation analysis of deterministic models
We present the results of a bifurcation analysis for the deterministic population based models.We refer to Guckenheimer and Holmes (1985), Kot (2001), Kuznetsov (2004) for an introduction into bifurcation analysis.The reference parameter values are given in Table 2.
Phase-space plots are shown in Fig. 1 for the RM-model (top panel), PPh-model (middle panel) and PPd-model (bottom panel).Two values for are used: k 12 = 1/3 (left panels) with a stable equilibrium and k 12 = 1/8 (right panels) with a stable limit cycle.
The phase-space plots for the RM-model are shown in Fig. 1a,b where we take the same initial values as for the other deterministic models, that is x 1 (0) = 0.2 and x 2 (0) = 0.7.This model in described by the set of ODEs in Eq. ( 6) together with Eq. ( 8).
Figure 1 c,d displays the numerical simulation results for the PPhmodel.In the phase-space diagrams because C 0 = 1 and x 0 > 0 the dynamics is limited to the left-lower triangle of the state space (0 ≤ x 1 ≤ 1, 0 ≤ x 2 ≤ 1) while for the RM-model there is not this restriction.The amplitude of the limit cycle is lower that that in the RM-model and especially the bottom value of the cycle is larger.
Figure 1 e,f displays the numerical simulation results for the PPdmodel ( 18).The additional parameter is y ε 12 = 1/4, so larger than 1/ 300 the empirical value for the lowest level of food chains and webs in Eq. ( 20).
In Table 3 the explicit expressions for the equilibria of the RM-model and PPh-model are given, see also Kooi and Poggiale (2018).Since these models are two-dimensional it is easy to perform a bifurcation analysis.A Hopf bifurcation occurs when the trace of the Jacobian matrix evaluated at equilibrium E 2 is zero and a transcritical bifurcation TC when its determinant is zero, see (Kot, 2001, p. 134).This TC bifurcation occurs for parameter values where x 2 = 0, hence the predator-prey equilibrium Fig. 2. One-parameter bifurcation diagram for k 12 .All other parameter values are given in Table 2. a) the RM-model Eq. ( 6), b) PPh-system Eq. ( 5), c) PPd-system Eq. ( 18).In the bottom-subpanels the prey population x 1 is depicted and in the top-subpanels the predator 2 .Red lines indicate stable equilibria or maxima and minima values of limit cycles and blue lines unstable equilibria.The prey-only equilibria E 1 , see Table 3 (above transcritical bifurcation TC) and interior predatorprey E 2 equilibria (between TC and Hopf bifurcation H) and the maximum and minimum values for the limit cycle L 2 (above H) are shown.((For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)E 2 coincides with the prey-only equilibrium E 1 .For the PPh-model this means that x 1 = C 0 and hence x 0 = 0 as a result of the assumption that m 1 = 0.This is biological realistic with instantaneous remineralisation.
The analytical expressions were obtained with Maple (2017) and shown also in Table 3 with the exception of the three dimensional PPd-model because in that case the expressions are unmanageable.
We perform a numerical bifurcation analysis with the parameter values from Table 2. One-and two-parameter diagrams are calculated with the computer package AUTO Doedel and Oldeman (2009).Continuation of the parameter k 12 the gives the bifurcation diagrams shown in Fig. 2.There are two bifurcations involved.The transcritical bifurcation TC and the Hopf bifurcation H. Below the transcritical bifurcation TC only the prey population exists in stable equilibrium E 1 .There is coexistence in stable equilibria E 2 in the interval between transcritical bifurcation TC and Hopf bifurcation H and above the Hopf bifurcation as stable limit cycles.
These numerical results for the RM-model in Fig. 2a show that the limit cycle for parameter values above the supercritical Hopf bifurcation H is stable and that the minimum values become very small for small k 12 .This phenomenon is related to the "paradox of enrichment" (Dawes and Souza, 2013;Rosenzweig, 1971), because catastrophic extinction due to stochastic fluctuations is likely, as we will see in the discussion of the PPs-model results below.
For the PPh-model by the bifurcation diagram is shown in Fig. 2a where the additional parameter is y ε 12 = 1/4.This diagram looks the same as Fig. 2b for the RM-model.The main difference is that the Hopf bifurcation H is at somewhat higher k 12 value and the predator biomass is smaller and this holds similarly for the amplitude of the limit cycle.
For the PPd-model the one-parameter bifurcation diagram is shown in Fig. 2c.This diagram looks again similar to the one in Fig. 2b for the PPh-model.The main difference is that the transcritical bifurcations TC and Hopf bifurcation H occur at lower k 12 values.Furthermore the equilibrium E 2 values and the amplitude of the limit cycle are smaller.
At the Hopf bifurcation the real parts of two conjugate eigenvalues of the Jacobian are zero while at the transcritical bifurcation one eigenvalue is zero and the other is negative.This indicates that there is a focus-node point that plays the same role as a bifurcation point.At this point two real eigenvalues are equal and negative.It marks the transition between a stable node and a focus (spiral sink).Hence the qualitative long-term dynamics of the deterministic system does not change when the bifurcation parameter crosses this point; in both cases there is convergence to the equilibrium.We will see below that this point is important when stochasticity is taken into account.
The two-parameter diagrams are shown in Fig. 3.For I 12 = 5/3 the results are shown in the one-parameter diagram Fig. 2 for parameter k 12 .There are two bifurcation curves TC and H. Below the transcritical bifurcation TC only the prey population exists in stable equilibrium E 1 .There is coexistence in stable equilibria in the region denoted by E 2 between Hopf bifurcation H and transcritical bifurcation TC curve and as stable limit cycles denoted L 2 in the region above the H curve.
The analytical expressions for the two-dimensional models: RMmodel ( 6) and PPh-model ( 5) are given in Table 4.For the TC curve is fixed by the condition that the determinant of the Jacobian matrix evaluated in the equilibrium is zero.In the case of the Hopf bifurcation it is that its trace equals zero.For the three-dimensional PPd-model the analytical expressions are not shown because they are very complicated.
The bifurcation portrait (equilibria prey population, equilibria preypredator equilibrium and limit cycles) is, with the reference parameter values from Table 2, almost the same for the three models.
Modelling mechanistic mass conservation in the PPh-model with respect to the descriptive RM-model logistic growth rate of the prey population gives limited changes in the position of the Hopf bifurcation H that separates the regions with stable coexistence in equilibria E 2 and on limit cycles L 2 , see Fig. 3b versus Fig. 3a.Introducing additionally two compartments for the predator as searchers and handlers in the PPdmodel alters the position of the transcritical bifurcation TC, see Fig. 3c versus Fig. 3b.

Stochastic population PPs-model
We start with a heuristic derivation of the PPs-model.The PPd-model Eq. ( 18) describes the dynamics of the predator-prey dynamics in the closed environment whereby the predator population consists of two compartments, searchers and handlers, hence a three-dimensional system for x = (x 1 , x 2s , x 2h ) ∈ R 3 + , the deterministic variables.In the derivation of the PPs-model the three model variables will be expressed intermediately in numbers of individuals instead of biomasses via n 1 = x 1 /x ε 1 and n 2 = x 2 /x ε 2 .Later in the final PPs-model we will again use biomasses with X = (X 1 , X 2s , X 2h ) ∈ R 3 + now random variables where X 2 = X 2s + X 2h .Mass conservation can be used to calculate X 0 = C 0 − X 1 − X 2 which is now also a stochastic variable.
However, at random events these quantities will change in finite steps with sizes of the prey x ε 1 and predator x ε 2 individuals.We adhere to

Table 2
Parameters and state variables for RM-, PPh-, PPd-, and PPs-models; : t=time, m=biomass, #=numbers, v=volume of the reactor.Reference values and dimensions are given.For the subindices we have: Carrying capacity of prey

Table 3
Equilibria of RM-model ( 6) and PPh-model ( 5) for reference values given in Table 2 where I 12 = 5/3.When the prey-predator equilibrium E 2 is unstable a stable limit-cycle denoted by L 2 exists.In the plots only the maximum and minimum values are shown.
) Prey-predator N. Stollenwerk et al. predator population goes to the origin x 1 = x 2 = 0.In a similar way the predator population can go extinct (in the last two loss processes in Table 5 x 2 becomes negative, that is zero) and the prey population goes deterministically to the origin x 1 = C 0 ,x 2 = 0, possible because m 1 = 0.This is related to the "paradox of enrichment" (Rosenzweig, 1971) phenomenon.
We compare the stochastic PPs-model realizations with their deterministic counterparts PPd-model where x ε 1 = 0.0002 and y ε 12 = 1 /4 in Fig. 5.In the upper panels of Fig. 5a,b we have k 12 = 1/3 a stable equilibrium and in lower panels of Fig. 5c,d with k 12 = 1 /8 an unstable equilibrium and stable limit cycle.
Again the transient dynamics for both models are similar but the asymptotic dynamics is not.After the relatively small stochastic fluctuations along the trajectory they are much larger when the trajectory starts circling around the equilibrium point of the PPd-model on the x 1 -nullcline.
Figure 6 gives one realization with the long-term result after transient dynamics shown in Fig. 5a.Also the quasi-equilibrium distribution approximated by the multivariate normal distribution is shown.In Appendix B the expressions for the 80% confidence ellipse of this distribution are derived.

Discussion
A well known classical predator-prey model is the deterministic RMmodel (Rosenzweig, 1971) and its dynamics has been studied intensively, we mention only (Kot, 2001) with further references there.Stochastic formulations have also been studied in the literature, we mention for instance (Allen, 2000;Dawes and Souza, 2013;Gardiner, 2009;Gurney and Nisbet, 1998;van Kampen, 1992;McKane and Newman, 2005;Renshaw, 1991;Stollenwerk et al., 2017).The elementary model formulation for a single population is by birth-death processes.This method has been extended for the predator-prey models such as the Lotka-Volterra predator-prey (McKane and Newman, 2005) and also the RM-model with Holling type II functional response in Renshaw (1991), Dawes and Souza (2013).
In this paper we discuss an alternative mass-balance PPh-model for the RM-model, representing the lowest levels of a aquatic food chain.We distinguish the individual and population level and use the 'number of individuals' besides the 'population biomass' description of the state  18).The additional parameters for the PPs-models are y ε 12 = 1/4 and x ε 1 = 0.0002.With k 12 = 1/3 (red curve) there is "convergence" to the stable equilibrium of the PPd-model and with k 12 1 /8, (green curve) to the limit cycle.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 6.Phase-space analysis for the PPs-model k 12 = 1/3 with additional parameter values y ε 12 = 1/4 and x ε 1 = 0.0002.The long-term dynamics is depicted and the green curve is the 80% confidence ellipse described by Eq. (B.6).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) N. Stollenwerk et al. variables like in Gurney and Nisbet (1998), McKane and Newman (2005), Dawes and Souza (2013).
There is in equation for the prey growth of the PPh-model an extra term that shows that the prey and the predator population are competing for the ambient resources besides the direct predator-prey interaction.In the RM-model this term − x 1 x 2 is missing because of the descriptive nature of used logistic growth, although correct for the prey-resource system also used for food chains longer than two.Similar to our result in the mean field description obtained in Dawes and Souza (2013) in the prey growth rate also this extra term exists.
Comparing the results of the PPh-model with those of the RM-model shows that the equilibrium biomass of the predator is larger in the classical RM-model, see the one parameter bifurcation diagrams Fig. 2a and b.When the system oscillates, however, the amplitude of the limit cycles is smaller and the minimum values during these cycles is larger and therefore there is less prone to extension due to stochastic effects.Qualitatively the bifurcation pattern, the position of the Hopf and transcritical bifurcation curves in the parameter space, is for the two models the same.This is in agreement with results in Dawes and Souza (2013).
In Kooi and Poggiale (2018) the PPh-and RM-model are compared in a more extreme case where the prey and predator time scales differ orders of magnitude.The RM-model shows complex dynamics when passing the Hopf bifurcation.First the originating stable limit cycle at the Hopf bifurcation leads to a stable limit cycle of which the amplitude grows when the distance from the Hopf bifurcation increases.However, suddenly at one critical point changing the bifurcation parameter further, there is a canard explosion after which the dynamics is similar to the relaxation oscillations predicted in the limiting case where the ratio of the predator/prey time scales equals zero, see also Rinaldi and Muratori (1992).On the other hand, in the PPh-model this complex dynamics is absent.
In Gurney and Nisbet (1998), Dawes and Souza (2013), Stollenwerk et al. (2017) the Holling type II functional response expression is derived by an ad hoc introduction of a fast time trophic interaction part with the small ε in system (10).In these papers the dynamics of the fast searching-handling predators is considered independent of the interaction with the prey population.In this paper the parameter y ε 12 , being the prey-predator body mass ratio is introduced for the population biomass as state variable allowing to use experimental results in Jonsson and Ebenman (1998) to get realistic values in food chains.The introduction of the new parameter y ε 12 , revises the inconsistency of the formulation by applying the time-scale separation technique with the derivation of the Holling type II functional response for the predator-prey interaction.An assumption for this application is still that the feeding rate is fast with respect to the population dynamics rates.When an ecosystem is close to or in equilibrium these rates will be small but this can be violated when this ecosystem dynamics is erratic.This is more important for larger ecosystems than for just predator/prey systems.Other assumptions related to the use of compartmental models, space-homogeneity, and unstructured instead of age-, stage-or physiologically structured are not discussed.
With the formulation of the stochastic PPs-model the introduction of the prey individual biomass x ε 1 leads to three nested model formulations: PPh-, PPd-and PPs-model with increasing complexity.The stochastic model results are obtained using the direct method of Gillespie (2001) like in Stollenwerk et al. (2017) where exponentially distributed waiting times are used.This is in contrast with the technique used in Kooijman et al. (2007) which is directly based on the Poissonian processes for searching, handling or dying.Note that in McKane and Newman (2005) also both algorithms were used and they found excellent agreement.
The one-parameter bifurcation diagrams Fig. 2 show that for low k 12 values the amplitude of the limit cycles are large and furthermore that during this oscillatory dynamics the biomasses of the prey and predator become very low.This is similar to the "paradox of enrichment" Rosenzweig (1971) phenomenon.The results for the PPs-model in the parameter region where large oscillations occur, Fig. 4c shows indeed extinction due to stochastic fluctuations often claimed in the literature on deterministic predator-prey or ecosystem models in nutrient rich situations.When x ε 1 is small the stochastic predictions of the PPs-model are close to those of the deterministic model.For the low y ε 12 = 1/40 value, the results for the PPs-model are in Fig. 4 compared with those of the PPhmodel where in the limit y ε 12 = 0.For the higher y ε 12 = 1/4 value with those of the PPd-model in Fig. 5.This indicates that in case of realistic body size ratio's the PPd-model is the preferred mass-balance predatorprey model.
We also discussed the approximation of the asymptotic probability density function PDF by a multivariate normal distribution, Grasman and Herwaarden (1999), Kooijman et al. (2007), Buckingham-Jeffery et al. (2018).Furthermore these results show that the bifurcation occuring in the deterministic models have to be re-considered when stochasticity is taken into account.

Conclusions
In this paper we re-stated the law of mass action formulation in the case of the trophic interactions in the predator-prey system.To solve an inconsistency problem in the derivation of the classical Holling type II formulation where, when biomasses are used, the same term runs at two different time-scales.The formulation of the model using both biomass and number dimension yields a natural scaling factor in the stochastic model measuring the degree of mean-filed approximation, namely the prey individual biomass where the volume of the batch reactor is unit.Finally this approach is not only applicable for the modelling and analysis of predator-prey systems but also for larger ecosystem models.
The most important results obtained and answers to the research questions posed in the introduction, are: 1.There is an extra term in the growth rate of the prey population Eq.
(5a) in the mass balance model with respect to the logistic growth description Eq. (6a) in the RM-model.The biological interpretation is that the prey uses the ingested resources partly for its own but also an extra amount used as building blocks by the predators.Conservation law for the total system shows that the prey and predator are also (indirectly) competing for the ambient resources and therefore there are two negative interaction predator-prey terms.2. With formulating the PPd-model using biomass as descriptive variable for the size of the populations, an inconsistency occurs when a time scale argument is used to derive the Holling type II functional response expression.This is solved by using the numbers of (identical) individuals of the populations as well.Only when the ratio of the biomasses of prey and predator individuals used a time-scale parameter is small enough the time-scale-separation technique can be used correctly.3.There are no important differences in bifurcation pattern of the mechanistic mass-balance PPd-model and the descriptive RM-model.Only a transcritical bifurcation and a Hopf bifurcation occur.4. The ratio of the biomasses of prey and predator individuals is the time-scale separation parameter.5.When the ratio of the biomasses of prey and predator individuals is small the PPd-model gives very similar results as the PPh-model quantitatively in the phase-space: equilibria and limit cycles.Hence the use of the Holling type II expression is best for lower trophic levels and has to be justified for higher levels.6. Prey and predator individual-biomasses are used as additional parameters in the stochastic model.The flow of the predator-prey stochastic system form sustained counter-clock oscillations in the phase-space plot and they show cycling behaviour where the deterministic version predicts a stable equilibrium.
N. Stollenwerk et al. 7. The well known paradox of enrichments phenomenon has been shown.8.In the stochastic PPs-model no similar threshold effects occur at the transcritical and Hopf bifurcation points of the companion deterministic PPh-model.

Declaration of Competing Interest
None.
So, the expected value given X at time t yields E(ΔX|X) = b(X)Δt + O (Δt 2 ) where b is a vector-valued function formed by the right-hand side of the system of ODEs Eq. ( 18) for the deterministic PPd-model.This first moments of the stochastic variable X(t) is the drift b which describes the mean-field solution.
For the second moments we derive the variance-covariance matrix of the change in X which is Cov(ΔX,ΔX T ⃒ ⃒ X) = a(X)Δt + O (Δt 2 ), where a is the matrix-valued function: We do not formulate the initial value for this function (as is done in Dawes and Souza, 2013) since we are interested only in the stationary distribution of the stochastic variable X found from ∂f/∂t = 0.It can be approximated by solving the system in which the drift and diffusion are N. Stollenwerk et al. approximated with the quasi equilibrium state values, that is near the stable equilibrium X * = x * for the deterministic system dx /dt = bx.In McKane and Newman (2005), Anderson and Kurtz (2011) this is called the Langevin approximation but see also Buckingham-Jeffery et al. (2018) where it is called a Gaussian process approximation.We take a linear drift approximation for the deterministic part evaluated at the equilibrium of the x changes in the deterministic PPd-model and a constant diffusion for the stochastic part b(X) ≃ B(X − X * ) and a(X) = a(X * ) = A , (B.3) where B is the Jacobian matrix of the model evaluated at the stable equilibrium point of the PPd-model described by Eq. ( 18): X = X * , that is: B = d b T (X * )/dX.Locally in the vicinity of this stable equilibrium the system behaves as a three-dimensional Ornstein-Uhlenbeck process (Gardiner, 2009;Grasman and Herwaarden, 1999) with a multivariate normal distribution as stationary solution as a sum of that of many random variables and in accordance with the central limit theorem.Its covariance matrix S satisfies the matrix equation where χ 2 2,0.8 = 3.219 is when the Chi-square distribution with 2 degrees of freedom equals 0.8. Figure 6 gives one realization long-term result after transient dynamics in Fig. 5a and also the approximated quasi-equilibrium distribution.

Fig. 5 .
Fig. 5. Phase-space analysis with k 12 = 1/3 a,b) and k 12 = 1/8 c,d), PPs-model a,c) and PPd-model b,d) Eq. (18).The additional parameters for the PPs-models are y ε 12 = 1/4 and x ε 1 = 0.0002.With k 12 = 1/3 (red curve) there is "convergence" to the stable equilibrium of the PPd-model and with k 12 1 /8, (green curve) to the limit cycle.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) scaling factor.The ε-domain is given by the 80% confidence ellipsoid in the X-space with the stationary state in the center given by(X − X * ) T S − 1 (X − X * ) = x ε 8 = 4.641 is the value for which the distribution function of the Chi-square distribution with 3 degrees of freedom equals 0.8.To construct the 80% confidence ellipse in the (x 1 , x 2 )-plane we carry out the transformationY = M(X − X * ) with M = ( that Y T = (x 1 , x 2 ).Eq. (B.5) then becomes Y T ( MSM T ) − 1 Y = x ε