# Example 9: Lotka-Volterra Dynamic Population Model

## Example 9: Prob-C Code for the Lotka-Volterra Dynamic Population Model

In their overview article about probabilistic programming the authors present a series of paradigmatical examples (Gordon, Henzinger, Nori & Rajamani, Probabilistic Programming, 2014). Example 9 deals with the classical predator-prey population model known as Lotka-Volterra model. The convential model is formulated as system of deterministic differential equations. Here the authors offer a stochastic solution in form of a probabilistic C-like computer code (Prob-C). The underlying mathematical model is a Continous-Time Markov Chain (CTMC).

"The Lotka-Volterra predator-prey model (Lotka, 1925; Volterra, 1926) is a population model which describes how the population of a predator and prey species evolves over time. It is specified using a system of so called 'stoichiometric' reactions as follows:

$$G \longrightarrow 2G$$

$$G+T \longrightarrow 2T$$

$$T \longrightarrow 0$$

We consider an ecosystem with goats (denoted by G) and tigers (denoted by T). The first reaction models goat reproduction. The second reaction models consumption of goats by tigers and consequent increase in tigers. The third reaction models death of tigers due to natural causes.

It turns out that this system can equivalently modelled as a Continuous Markov Chain (CTMC) whose state is an ordered pair (G, T) consisting of the number of goats G and the number of tigers T. The first reaction can be thought of as a transition in the CTMC  from state (G, T) to (G + 1, T) and this happens with a rate equal to $$c_1 \cdot G$$ where c1 is some rate constant, and is enabled only when G > 0. Next, the second reaction can be thought of as a transition in the CTMC from state (G, T) to (G-1, T+1) and this happens with a rate equal to $$c_2 \cdot G \cdot T$$, where c2 is some rate constant, and is enabled only when G > 0 and T > 0. Finally , the last reaction can be thought of as a transition in the CTMC from state (G, T) to (G, T-1) and this happens with a rate equal to $$c_3 \cdot T$$ where c3 is some rate constant, and is enabled only when T > 0.

Using a process called uniformization, such a CTMC can be modeled using an embedded Discrete-Time Markov Chain (DTMC), and encoded as a probabilistic program, shown in Figure 1. The program starts with an initial population of goats and tigers and executes the transitions of the Lotka-Volterra model until a prescribed time limit is reached, and returns the resulting population of goats and tigers.

Since the executions are probabilistic, the program models the output distribution of the population of goats and tigers. The body of the while loop is divided into 3 conditions:

1. The first condition models the situation when both goats and tigers exist, and models the situation when all 3 reactions are possible.
2. The second condition models the situation when only goats exist, which is an extreme case, where only reproduction of goats is possible.
3. The third condition models the situation when only tigers exist, which is another extreme case, where only death of tigers is possible." (Gordon et al., 2014, p.6)

## Example 9: Formal Reconstruction of Gordon's et al. Prob-C Code

We think that the so called 'stoichiometric' reactions are misleading because they are rather different from chemical stoichiometric 'equations' or 'patterns'. According to Gordon et al. agent 'molecules' can appear out of nothing or disappear into nothing. This is not allowed in classical stoichiometry (Erban & Jonathan, 2020).

So we prefer to model the Lotka-Volterra-Scenario directly as a competitive Continuous Time Markov Chain (CTMC) with countable state space along its definition (e.g. Pishro-Nik, 2014, p.666; Ross, 2010, 10/e, p.373) with two subprocesses:

1. a holding-time process, when in state i
2. a jump process from state i to a different state j

Besides their standard use in natural (Erban & Chapman, 2020) and computational (Gordon, 2014; Mitzenmacher & Upfahl, 2018; Cassandras & Lafortune, 2008) sciences CTMCs are a non-standard modelling tool in psychology (Wickens, 1982, McGill, 1967; Diederich & Mallahi-Karai, 2020) and sociology (Coleman, 1964; Bartholomew, 1967; Tuma & Hannan, 1984).

The holding-time Poisson counting process is waiting for an event of type birth (of prey) or fatal attack (of predator), or death (of predator). In this kind of process it is assumed that holding or waiting times are exponential distributed with rates:

$$\lambda(G,T) = \lambda_1(G) +\lambda_2(G,T) + \lambda_3(T).$$

$$\lambda_1(G) = c_1 \cdot G \text{ ; birth process with birth events}$$

$$\lambda_2(G, T) = c_2 \cdot G \cdot T \text{ ; birth-death process with attack events}$$

$$\lambda_3(T) = c_3 \cdot T \text{ ; death process with death events}$$

Because rates are state-dependent the Poisson counting processes are inhomogeneous. The total rate is a sum of the rates of the three individual and independent Poisson counting subprocesses (birth, birth-death, and death) running parallel and in competion to each other. Poisson race processes are widely used in such diverse fields as cognitive reaction time (RT), choice processes or reliability studies (Ibe, 2013, ch. 2.7.5; 2014, ch. 12.5.9). In independent Poisson race processes the total rate is the sum of the rates of the subprocesses and the probability of winning this race can be determined by the ratio:

$$p_{j} = \frac{\lambda_{j}}{\lambda} \text{ ; j = 1,...,3} \text{ ; probability of j is winning the process race}.$$

As the interarrival times of a Poisson counting process have an exponential distribution we can sample holding times ('dwellTime') of the merged Poisson process

$$\Delta t(G, T) \sim DExp(\lambda(G,T)) = DExp('rate')$$

add them up to the current time ('curTime'), and compare the current time with our TIMELIMIT.

At the end of the holding-time interval the CTMC has to determine which of the three subprocesses generated the event. The probability of winning the race between the three Poisson counter processes and the  jump-to probability is the simple standardized rate parameter (Ibe, 2013, ch. 2.7.5; 2014, ch. 12.5.9):

$$p_{ij}(G,T) = \frac{\lambda_{ij}(G,T)}{\lambda_i(G,T)} \text{ ; i = 1, 2, ..., ; j = 1,...,3}.$$

Since the rates and probabilities for all i are constant ("uniform"), we suppress the state index i from now on.

$$p_{j}(G,T) = \frac{\lambda_{j}(G,T)}{\lambda(G,T)} \text{ ; j = 1,...,3} \text{ ; probability that j is winning the process race}.$$

The jump-to index j is sampled from discrete distribution with the jump probabilities as parameter vector

$$j(G,T) \sim Discrete(p_1(G,T), p_2(G,T), p_3(G,T)) \text{ ; j = 1, ..., 3}$$

According to the selected jump-to index j the three processes generate at most two simultaneous actions. According to Gordon et al. these actions are:

$$(G, T) \stackrel{c_1 \cdot G}{\longrightarrow} (G+1, T) \text{ ; if j = 1, then birth process}$$

$$(G, T) \stackrel{c_2 \cdot G \cdot T}{\longrightarrow} (G-1, T+1) \text{ ; if j = 2, then birth-death process}$$

$$(G, T) \stackrel{c_3 \cdot T}{\longrightarrow} (G, T-1) \text{ ; if j = 3, then death process}$$

The simultaneous actions of the birth-death process (j = 2) are definitional for this kind of process, though they seem to be a bit unnatural in the Lotka-Volterra context. It is a bit unrealistic that during an exponential distributed attack episode a goat's disappearance is accompanied in the same episode by a birth of a tiger.

These transitions reach every state i in two-dimensional state-space S.

$$S \subseteq G \times T$$ $$G, T \subset \mathbb{N}$$

Transitions of the Lotka-Volterra CTMC are visualized in a transition-rate diagram (Fig. 2). This 2-dimensional transition space is typical for a competition process (Reuter, 1961; Iglehart, 1964) . It deviates from the conventional one-dimensional transition space typically used in CTMC definitions.

Empirical research has shown that the maximum of the predator frequency distribution appears most times later than the maximum of the prey distribution. This is not true for Gordon et al's model, when we use their rate-constants c as we demonstrate with our WebPPL-coded reconstruction (see below).

The classical Lotka-Volterra model is valuable for educational purposes but too too simple for empirical predictions and as an agent-based model it is not embedded in an ecological context (Charles J. Krebs, Rudy Boonstra, Stan Boutin, A.R.E. Sinclair, 2001). This is the fundamental deficit.

## Lotka-Volterra Simulation Run with WebPPL Script and Mateos' Parameters

Mateos' parameters were obtained from his 2018 lecture slides (visited 2020/07/27)

$$c_1 = 1.0$$

$$c_2 = 0.005$$

$$c_3 = 0.6$$

$$Goats(0) = 50$$

$$Tigers(0) = 100 .$$

The WebPPL-Code was developed independently. Results are similar though some of our simulations runs end in an ex- or implosion of the state space.

## Lotka-Volterra Simulation Run with WebPPL Script and Gordon et al.'s Parameters and Rates

Gordon et al.'s parameters and rates were obtained from the authors 2014 publication (visited 2020/07/28).

$$c_1 = 1.0$$

$$c_2 = 5.0$$

$$c_3 = 1.0$$

$$Goats(0) = 100$$

$$Tigers(0) = 4 .$$

The pure functional WebPPL-Code is an independently developed reimplementation of the imperative Prob-C code . Results demonstrate that the model does not show the typical Lotka-Volterra behavior with its cyclic trajectories. A model run with Mateos' start parameters and rate are much more typical for the Lotka-Volterra model trajectories (see above).

## References

Bartholomew, D.J., Stochastic Models for Social Processes, John Wiley, 1/e 1967, 3/e 1982

Cassandras Ch.G. & Lafortune, St., Introduction to Discrete Event Systems, 2/e, Springer, 2008

Diederich, A. & Mallahi-Karai, K., Stochastic Methods for Modelling Decision-Making, p.1-70, in: Batchelder, W.H., Colonius, H. & Dzhafarov, E.N., New Handbook of Mathematical Psychology, Vol.2: Modeling and Measurement, Cambrige University Press, 2018

Erban, R. & Chapman, J., Stochastic Modelling of Reaction-Diffusion Processes, Cambrige University Press, 2020

Gordon, A. D., et al. Probabilistic programming. In: Proceedings of the on Future of Software Engineering. 2014. S. 167-181.

Ibe, O., Markov Processes for Stochastic Modeling, 2/e, Elsevier, 2013

Ibe, O., Fundamentals of Applied Probability and Random Processes, Academic Press, 2014

Iglehart, D. L., et al. Multivariate competition processes. The Annals of Mathematical Statistics, 1964, 35. Jg., Nr. 1, S. 350-361.

Krebs, Ch. J., et al. What Drives the 10-year Cycle of Snowshoe Hares? The ten-year cycle of snowshoe hares—one of the most striking features of the boreal forest—is a product of the interaction between predation and food supplies, as large-scale experiments in the Yukon have demonstrated. BioScience, 2001, 51. Jg., Nr. 1, S. 25-35.

Lotka, A., Elements of physical biology, Williams & Wilkins company, Baltimore, 1925.

Mateos, G., Predator-Prey Population Dynamics, Lecture Slides, 2018/11/13, Dept. of ECE and Goergen Institute for Data Science, University of Rochester (visited 2020/07/28)

McGill, W.J., Stochastic Latency Mechanisms, p.309-360, in: Luce, R.D., Bush, R.R. & Galanter, E. (eds), Handbook of Mathematical Psychology, Vol. I, John Wiley & Sons, 2/e, 1967

Mitzenmacher, M. & Upfahl, E., Probability and Computing; Randomization and Probabilistic Techniques in Algorithms and Data Analysis, 2/e, 2017

Pishro-Nik, H., Probability, Statistics, and Random Processes, Kappa-Research LLC, 2014

Reuter, G. E. H. Competition processes. In: Proc. 4th Berkeley Symp. Math. Statist. Prob. 1961. S. 421-430.

Ross, S.M, Introduction To Probability Models, Elsevier, 2010, 10/e

Ruan, Sh., Poisson Race Models: Theory and Application in Conjoint Choice Analysis. 2007. PhD-Thesis. The Ohio State University.

Tuma, N.B. & Hannan, M.T., Social Dynamics: Models and Methods, Academic Press, 1984

Volterra, V., Fluctuations in the abundance of a species considered mathematically, Nature, 118: 558–560, 1926.

Wickens, Th.D., Models For Behavior: Stochastic Processes in Psychology, Freeman & Co., 1982

(Changed: 2020-08-05)