Lotka Volterra Competition Model in Python - Part I

I have seen that exist many websites where Lotka-Volterra Predator-Prey model is explained detailed, and even with simulations in Python. However, it is less common to see sites explaining the Competition model. Then, I will try to make a very (very) simple explanation of this model as well as a simulation in Python (of course). You won't see much mathematical rigour (I am very sorry for that) as I intend to show this mainly to biologists.

Lotka-Volterra Competition model is the most famous mathematical approach of the interaction phenomenon of ecological competition. The generalised version of the model is defined as follows:

$\frac{dN_{i}}{dt} = r_{i}N_{i} (\frac{K_{i} - N_{i} - \sum_{j \neq i}^{n} \alpha_{ij}N_{j}}{K_{i}})  (1)$

Where:

$N$: Population size of i, j

$r$: Intrinsic growth rate of population i

$K$: Carrying capacity of population i

$\alpha_{ij}$: Competition effect of population j on population i

It could be a bit difficult to understand if you are not used to dealing with this kind of problems. We are going to simplify into a scenario where only two populations are competing, so we have now the following model:

$\frac{dN_{1}}{dt} = r_{1}N_{1}(1 - \frac{ [N_{1} + \alpha N_{2}]}{K_{1}}) (2)$

$\frac{dN_{2}}{dt} = r_{2}N_{2}(1 - \frac{[N_{2} + \beta N_{1}]}{K_{2}}) (3)$

For a biological context, we are going to assume that populations are always positive, so $N_{1}, N_{2} \in \mathbb{R}_{0}^{+}$. Also, we assume that populations will growth until a non-trivial carrying capacity, so $K_{1}, K_{2} \in \mathbb{R^{+}}$ and $r_{1}, r_{2} \in \mathbb{R^{+}}$. Finally, we are going to assume that competition between populations exist, then $\alpha, \beta \in \mathbb{R^{+}}$. Now, we can start our analysis.

Firstly, you have to understand is that the above equations are a system of ordinary differential equations. There is an extensive framework on this topic that I do not pretend to cover here, but if you are interested in this field, I recommend you to check this book:

Zill DG. (2013). A first course in differential equations with modelling applications. 10th edition. Brooks/Cole Cengage Learning. 384pp.

The system of differential equations in (2) and (3) represent how population size of competitors 1 and 2 changes over time. Both competitors have a logistic growth, but there is something else. If you focus on equation (2) you can see that the variable $N_{2}$ is present, so it means that the population size of competitor 1 will depend on the population size of competitor 2, and the other way around. Biologically that makes sense.

There exist three basic steps for analysing a system of ODEs:

1) Finding critical points. It consists of finding the conditions for which the system of ODEs remains constant. In other words, when the whole system of differential equations is equal to zero. Thus each equation must be equal to zero, for which we obtain the following conditions:

For equation (3)

$\frac{dN_{1}}{dt} = r_{1}N_{1}(1 - \frac{ [N_{1} + \alpha N_{2}]}{K_{1}}) = 0$

$r_{1}N_{1}(1 - \frac{ [N_{1} + \alpha N_{2}]}{K_{1}}) = 0$

$r_{1}N_{1} = 0$ or $(1 - \frac{ [N_{1} + \alpha N_{2}]}{K_{1}}) = 0$

$N_{1} = 0$ or $K_{1} - N_{1} - \alpha N_{2} = 0$

For equation (4)

$\frac{dN_{2}}{dt} = r_{2}N_{2}(1 - \frac{[N_{2} + \beta N_{1}]}{K_{2}})  = 0$

$r_{2}N_{2} (1 - \frac{[N_{2} + \beta N_{1}]}{K_{2}})  = 0$

$r_{2}N_{2} = 0$ or $(1 - \frac{[N_{2} + \beta N_{1}]}{K_{2}})  = 0$

$N_{2} = 0$ or $K_{2} - N_{2} - \beta N_{1}  = 0$

Using the above conditions, we can obtain the critical points of our system:

$P1 = (N_{1} = 0, N_{2} = 0) = (0, 0)$ Trivial

$P2 = (N_{1} = 0, K_{2} - N_{2} - \beta N_{1} = 0)$

$P2 = (0, K_{2} - N_{2} - \beta (0) = 0) = (0, K_{2} - N_{2} = 0)$

$P2 = (0, K_{2})$

$P3 = (K_{1} - N_{1} - \alpha N_{2} = 0, N_{2} = 0)$

$P3 = (K_{1} - N_{1} - \alpha (0) = 0) = (K_{1} - N_{1} = 0, 0)$

$P3 = (K_{1}, 0)$

$P4 = (K_{1} - N_{1} - \alpha N_{2} = 0, K_{2} - N_{2} - \beta N_{1} = 0)$

$P4 = (N_{1} = K_{1} - \alpha N_{2}, N_{2} = K_{2} - \beta N_{1})$

We need some extra arithmetic operation to obtain the P4

$N_{1} = K_{1} - \alpha [K_{2} - \beta N_{1}]$

$N_{1} = K_{1} - \alpha K_{2} + \alpha \beta N_{1}$

$N_{1} - \alpha \beta N_{1} = K_{1} - \alpha K_{2}$

$N_{1} (1 - \alpha \beta) = K_{1} - \alpha K_{2}$

$N_{1} = \frac{K_{1} - \alpha K_{2}}{1 - \alpha \beta}$

$N_{2} = K_{2} - \beta [K_{1} - \alpha N_{2}]$

$N_{2} = K_{2} - \beta K_{1} + \beta \alpha N_{2}$

$N_{2} - \beta \alpha N_{2} = K_{2} - \beta K_{1}$

$N_{2}(1 - \beta \alpha) = K_{2} - \beta K_{1}$

$N_{2} = \frac{K_{2} - \beta K_{1}}{1 - \beta \alpha}$

Thus, the critical points of the system are

$P1 = (0, 0)$ Trivial

$P2 = (K_{1}, 0)$ Competitor 1 survive, Competitor 2 is extinct

$P3 = (0, K_{2})$ Competitor 1 is extinct, Competitor 2 survive

$P4 = (\frac{K_{1} - \alpha K_{2}}{1 - \alpha \beta}, \frac{K_{2} - \beta K_{1}}{1 - \beta \alpha})$ Coexistence

At this point is a good idea to code. We are going to start by obtaining a solution to the model using arbitrary parameters and initial conditions. For example, I will simulate the model using:

Initial conditions: 

$N1_{0} = 6$

$N2_{0} = 5$

Parameters (arbitrary just for demonstration)

$r_{1} = 0.1$

$r_{1} = 0.1$

$\alpha = 0.4$

$\beta = 0.66$

$K_{1} = 100$

$K_{2} = 100$

	# 1. import libraries
	import numpy as np
	import matplotlib.pyplot as plt
	import scipy as sci
	from scipy import integrate
	#2. Creating the model function
	def LVC(N, t, r1, r2, alpha, beta, K1, K2):
	"""Lotka-Volterra Competition model"""
	dN_1 = r1*N[0] * (1 - (N[0] + alpha * N[1])/ K1)
	dN_2 = r2*N[1] * (1 - (N[1] + beta * N[0])/ K2)
	return(np.array([dN_1, dN_2]))
	#3. Set up initial conditions and parameters
	r1 = 0.1
	r2 = 0.1
	alpha = 0.4
	beta = 0.66
	K1 = 100
	K2 = 100
	NI = np.array([6, 5]) # This is a temporary variable which store N1 and N2 values
	#4. Time vector
	t = np.linspace(0, 150, 1000)
	#5. Solving the ODE
	sol = sci.integrate.odeint(LVC, NI, t, args=(r1, r2, alpha, beta, K1, K2))
	#6 Creating the graph
	plt.plot(t, sol[:, 0], color = "blue", label = "Competitor 1")
	plt.plot(t, sol[:, 1], color = "red", label = "Competitor 2")
	plt.legend(loc = "upper right")
	plt.xlabel("Time")
	plt.ylabel("Population Size")
	plt.show()

view raw gistfile1.txt hosted with ❤ by GitHub

Using this arbitrary initial condition and parameters, we obtain that the two species will coexist (i.e. both will reach stable non zero population in time) what is cool. However, it is even more interesting to determine the exact condition which determines the competitors coexist or one of them will get extinct. To do that, we will proceed to study the phase diagram.

2) Studying the Phase Diagram.

A phase diagram could be understood as a representation of the variables under study, the critical points and the trajectories of the solutions as vectors. This diagram also includes the nullclines of the system. Nullclines are trajectories where the derivate of a variable is equal to zero.

For example, for the Competitor $N_{1}$ the nullcline will be:

$N_{1} = 0$ a constante function (a line)

$K_{1} - N_{1} - \alpha N_{2} = 0$ a straight line

Analogously, the nullclines for Competitor $N_{2}$ will be:

$N_{2} = 0$ a constante function (a line)

$K_{2} - N_{2} - \beta N_{1} = 0$ a straight line

If we set up the X axis as $N_{1}$ and Y axis as $N_{2}$. Thus, the nullclines of Competitor $N_{1}$ would be a horizontal line ($N_{1}=0$) and a straight line ($K_{1} - N_{1} - \alpha N_{2} = 0$). Focusing in the latter we can determine where this line will cross the axes:

Making $N_{1} =0$, it follows:

$K_{1} - N_{1} - \alpha N_{2} = 0$

$K_{1} - 0 - \alpha N_{2} = 0$

$K_{1} =  \alpha N_{2}$

$N_{2} = \frac{K_{1}}{\alpha}$

Thus, when $N_{1} = 0$ the nullciline of competitor 1 will cut the $N_{2}$-axis in the value $N_{2} = \frac{K_{1}}{{\alpha}}$

Likewise, making $N_{2}=0$, it follows:

$K_{1} - N_{1} - \alpha N_{2} = 0$

$K_{1} - N_{1} - \alpha (0) = 0$

$N_{1} = K_{1}$

Thus, when $N_{2} = 0$ the nullcline of competitor 1 will cut the $N_{1}$-axis in the value $N_{1} = K$. Using these points we can define the straight line of the nullcline of competitor 1:

As you can see it is a straight line which cut the N_{1}-axis in $N_{1} = K = 100$ and the N_{2}-axis in the value $N_{2} = \frac{K_{1}}{{\alpha}} = 100/04 = 250$

Now, we must estimate the behaviour of the solution on the two regions separated for the nullcline. We can do a simple analysis for this:

If we only focus our attention in the $N_{1}$, we can see that $K_{1}$ marks the point where the nullcline divide the area in two regions (right and left). Thus, we can choose arbitrary values on each side of $K_{1}$.

- For instance, consider a value $N_{1}^{+} > K_{1}$. Replacing this value on (2), we can determine the derivate and then, the tendency of the solution:

$\frac{dN_{1}}{dt} = r_{1}N_{1}(1 - \frac{ [N_{1} + \alpha N_{2}]}{K_{1}})$

$\frac{dN_{1}^{+}}{dt} = r_{1}N_{1}^{+}(1 - \frac{ [N_{1}^{+} + \alpha N_{2}]}{K_{1}})$

As we are working on the $N_{1}-axis$, then $N_{2} = 0$

$\frac{dN_{1}^{+}}{dt} = r_{1}N_{1}^{+}(1 - \frac{ [N_{1}^{+} + \alpha (0)]}{K_{1}})$

$\frac{dN_{1}^{+}}{dt} = r_{1}N_{1}^{+}(1 - \frac{ N_{1}^{+}}{K_{1}})$

$\frac{dN_{1}^{+}}{dt} = r_{1}N_{1}^{+}(\frac{K_{1} - N_{1}^{+}}{K_{1}})$

We know that $r_{1}, K_{1},  N_{1}^{+} > 0$. Then:

$\frac{dN_{1}^{+}}{dt} =  (+) . (+) (\frac{K_{1} - N_{1}^{+}}{+})$

As $N_{1}^{+} > K_{1}$, it follows:

$\frac{dN_{1}^{+}}{dt} =  (+) . (+) (\frac{-}{+})$

$\frac{dN_{1}^{+}}{dt} =  (-)$

Thus, the derivate of the point N_{1}^{+} which is on the right side of the K_{1} will be negative, indicating the population will decrease. It can be represented as a vector pointing backwards. We can do a similar analysis to determine the behaviour on a point $N_{1}^{-} < K_{1}$, which is in the left side of the carrying capacity. The procedure will be similar:

$\frac{dN_{1}^{-}}{dt} =  (+) . (+) (\frac{K_{1} - N_{1}^{-}}{+})$

$\frac{dN_{1}^{-}}{dt} =  (+) . (+) (\frac{+}{+})$

$\frac{dN_{1}^{-}}{dt} =  (+)$

And here, the derivate is positive, which we can interpret that in this region population tend to increase. It can be represented as a vector pointing forward. This analysis can be replicated to all the points in the phase diagram, obtaining a vectorial field like this:

We can do an analogous analysis for the competitor $N_{2}$, obtaining a vectorial field like this:

But remember that these nullclines and victual fields are representations of each competitor in the absence of the another. However, in the competition model, we need to consider these figure simultaneously in the same figure. For the particular parameters and initial conditions we choose at the begging, we will obtain:

The above chart represents the Phase Diagram of the Competition model. The lime line indicate a particular solution for the initial conditions $N_{1}(0) = 6$, $N_{2}(0) = 5$. Pay attention to how all the vector point to a single point, which is not trivial. It means that under these combinations of values of parameters, all solutions always will end in the coexistence point. Precisely, it happens with our initial conditions $N_{1}(0) = 6$, $N_{2}(0) = 5$, look how the lime line ends up in the coexistence point, which values coincide with our simulations in the first figure of this post.

Continuing our analysis, you should realise the behaviour of the solutions in the phase diagram will depend on the combinations of parameters we use. If we vary the parameters we can obtain the classic predictions of the competition showed in many books of ecology:

	# 1. Import libraries
	import numpy as np
	import matplotlib.pyplot as plt
	import scipy as sci
	from scipy import integrate
	#2. Define the ODE function
	def LVC(N, t, r1, r2, alpha, beta, K1, K2):
	"""Lotka-Volterra Competition Model"""
	dN_1 = r1*N[0] * (1 - (N[0] + alpha * N[1])/ K1)
	dN_2 = r2*N[1] * (1 - (N[1] + beta * N[0])/ K2)
	return(np.array([dN_1, dN_2]))
	#3. Parameters values and initial conditions
	r1 = 0.1
	r2 = 0.1
	alpha = 0.4 # Change this parameter to obtain different scenarios of LV model
	beta = 0.66 # Change this parameter to obtain different scenarios of LV model
	K1 = 100
	K2 = 100
	NI = np.array([6, 5]) # Initial conditions (Initial population sizes)
	#4. Time vector
	t = np.linspace(0, 150, 1000)
	#5. Solving the system of differential equations
	sol = sci.integrate.odeint(LVC, NI, t, args=(r1, r2, alpha, beta, K1, K2))
	#6 Create the chart
	#plt.plot(t, sol[:, 0], color = "blue", label = "Competitor 1")
	#plt.plot(t, sol[:, 1], color = "red", label = "Competitor 2")
	#plt.legend(loc = "upper right")
	#plt.xlabel("Time")
	#plt.ylabel("Population size")
	#plt.show()
	# Phase Diagram. Multiples scenarios
	N_1 = np.linspace(0, np.max([K1, K2/beta])+5, 100)
	N_2 = np.linspace(0, np.max([K2, K1/alpha])+5, 100)
	Null_N1 = (K1 - N_1)/alpha
	Null_N2 = K2 - beta * N_1
	N1_max = np.max(N_1) # Max value of competitor 1
	N2_max = np.max(N_2) # Max value of competitor 2
	x = np.linspace(0.1, N1_max, 20) # Obtaining a vector with "X (here 20, see third argument)" values from 0 to N1_max"
	y = np.linspace(0.1, N2_max, 20) # Obtaining a vector with "X (here 20, see third argument)" values from 0 to N2_max"
	xx, yy = np.meshgrid(x, y) # Creating a meshgrid with points where the vector will be drawn.
	NN1, NN2 = LVC((xx, yy), 0, r1, r2, alpha, beta, K1, K2) # Derivate each point of the Phase Diagram (it will create the vector field)
	norm = np.sqrt(NN1**2 + NN2**2)
	NN1 = NN1 / norm
	NN2 = NN2 / norm
	plt.quiver(xx, yy, NN1, NN2, cmap=plt.cm.gray) #
	plt.plot(N_1, Null_N1, color = "blue")
	plt.plot(N_1, Null_N2, color = "red")
	plt.plot(sol[:, 0], sol[:, 1], color = "lime")
	plt.xlim(0, N1_max)
	plt.ylim(0, N2_max)
	plt.title("Phase Diagram of Competition Model")
	plt.xlabel("Competitor 1")
	plt.ylabel("Competitor 2")
	plt.show()

view raw Phase Diagram hosted with ❤ by GitHub

Case 1:

$r_{1} = 0.1$

$r_{2} = 0.1$

$\alpha = 0.8$

$\beta = 1.2$

$K_{1} = 100$

$K_{2} = 100$

Case 2:

$r_{1} = 0.1$

$r_{2} = 0.1$

$\alpha = 1.3$

$\beta = 0.75$

$K_{1} = 100$

$K_{2} = 100$

Case 3:

$r_{1} = 0.1$

$r_{2} = 0.1$

$\alpha = 1.2$

$\beta = 1.3$

$K_{1} = 100$

$K_{2} = 100$

Case 4:

$r_{1} = 0.1$

$r_{2} = 0.1$

$\alpha = 0.8$

$\beta = 0.75$

$K_{1} = 100$

$K_{2} = 100$

We obtained nice figures, but what they mean? First of all, be aware that in all cases, the parameters are very similar; we are only changing the values of $\alpha$ and $\beta$. And you can see that changing these two parameters is enough to modify the qualitative behaviour of the system completely. Recalling the definition of $\alpha$ and $\beta$, they are:

$\alpha$ the effect of the competitor $N_{2}$ over N_{1}$

$\beta$ the effect of the competitor $N_{1}$ over N_{2}$

For the Case 1 ($\alpha = 0.8$, $\beta = 1.2$). We can see that $\beta > \alpha$, thereby the competitor $N_{1}$ have a major impact on competitor $N_{2}$ than $N_{2}$ over $N_{1}$. Biologically, it can be interpreted that $N_{1}$ is a better competitor than $N_{2}$. So, intuitively (thinking like an ecologist), you can conclude that $N_{1}$ will defeat $N_{2}$ in the ecological competition. And it is what happens in the respective figure (Figure Case 1). You can see that the vector (and the particular solution in lime) end up in $K_{1}$, meaning that under this conditions $N_{1}$ will reach its carrying capacity, and $N_{2}$ will be zero, which means that it will become extinct!

For the Case 2 ($\alpha = 1.3$, $\beta = 0.75$). We can see the opposite situation that Case 1. Here, $\beta > \alpha$ so the $N_{2}$ is the best competitor. Therefore, this population will reach its carrying capacity and $N_{1}$ will get extinct.

The Case 3 ($\alpha = 1.2$, $\beta = 1.3$) is more interesting. Here, both $N_{1}$ and $N_{2}$ are good competitors ($alpha$, $beta > 1$, it will be more clear in the next part), so will create it an unstable coexistence point; you can see it in the intersection of the nullclines. But, what population will win in the competitions? in this situation, both could win. What will determine the winner are the initial conditions (N_{1}(0), and N_{2}(0)). Once again, look at the vector, and you can see that depending on where you start the initial point, you can end in the $N_{1}$-axis or$N_{2}$-axis. In the particular initial contusion we choose: $N_{1}(0) = 6$, and $N_{2}(0) = 5$, we can see that the $N_{1}$ is the winner and $N_{2}$ get extinct. But, what would had happened if we choose initial conditions like $N_{1}(0) = 20$, and $N_{2}(0) = 40$? You can see the vector would lead us to a situation where $N_{2}$ would be the winner and $N_{2}$ would get extinct.

 

Finally, the Case 4 ($\alpha = 0.8$, $\beta = 0.75$). Here, both $N_{1}$ and $N_{2}$ are bad competitor ($alpha$, $beta < 1$), and it is curios to note that this produce an stable equilibria in the intersection of the nullclines.

Now, some final points. Remember the equilibrium point of the system:

$P1 = (0, 0)$ Trivial

$P2 = (K_{1}, 0)$ Competitor 1 survive, Competitor 2 is extinct

$P3 = (0, K_{2})$ Competitor 1 is extinct, Competitor 2 survive

$P4 = (\frac{K_{1} - \alpha K_{2}}{1 - \alpha \beta}, \frac{K_{2} - \beta K_{1}}{1 - \beta \alpha})$ Coexistence

- P1, P2, and P3 are present in all Cases studied, but P4 only appear in Case 3 and 4.

- Vector never get closer to the P1, so P1 is a "repulsor" fixed-point. In contrast, vectors can go nearer ("attractors") the fixed-points P2, P3 and P4, or getaway of them ("repulsors"), which will depend on the parameters' values of $\alpha$ and $\beta$. We can determine exactly the conditions in which this fixed-points will be attractors or repulsors throughout a qualitative analysis of the system, and precisely it is what we do in the next part.

3) Determining the stability of the system.

As this post is becoming too long, I will develop this analysis in a future post. It will involve the analysis of the Jacobian matrix, so we will need to discuss some topics of linear algebra.