# Solving and Plotting Mutualism Model in Python

I am a beginner in programming. I need to program a mutualism model of two species in python that would solve and graph using the following equations:

$$\frac{dN_1}{dt} = N_1(r_1 - e_1N_1 + \alpha _{12} N_2)$$ and $$\frac{dN_2}{dt} = N_2(r_1 - e_1N_2 + \alpha _{21} N_1)$$

$$N_1$$ is the population of Species 1. The values of $$\alpha_{12}$$ is the effect of Species 1 on Species 2 (vice versa) and would vary from 0-2. Meanwhile, $$r_1$$, $$r_2$$, $$e_1$$, and $$e_2$$ have fixed values.

I tried to make a similar code as the Lotka-Volterra Model, but it didn't work.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# parameters
r1 = 1.0  #fixed
r2 = 0.5  #fixed
e1 = 1.0  #fixed
e2 = 0.75 #fixed

alpha12 = 1.5 #vary from 0-2
alpha21 = 1 #vary from 0-2

N1_0 = 1
N2_0 = 1

# store initial values in an array
X0 = [N1_0,N2_0]

# The two equations are contained in dX
def mutualism(X,t):
N1, N2 = X
dX = np.zeros(2)
dX[0] = N1 * (r1 - (e1 * N1) + (alpha12 * N2)) # equation for dN1/dt
dX[1] = N2 * (r2 - (e2 * N2) + (alpha21 * N1)) # equation for dN2/dt
return dX

# set time length
t = np.linspace(0,100,300)

X = odeint(mutualism,X0,t)

N1 = X[:,0]; N2 = X[:,1]

plt.plot(t,N1, color='blue', lw=3)
plt.plot(t,N2, color='red', lw=3)
plt.show()


Can anyone advise me on what I'm doing wrong and how I could improve the code?

• @LutzLehmann hello! r1 and r2 denote the respective growth rates of the populations. These equations were given to me, so I need to work around it unfortunately. Nov 15, 2021 at 19:30
• The alpha terms are a positive feed-back to the growth of the other species, together this gives super-linear growth in both, with a divergence to infinity shortly after t=1.5. There is nothing wrong with the code, just with the model. Nov 15, 2021 at 20:07
• @LutzLehmann I see! I was wondering why it seemed to bug after a certain point. Thank you for the insight, it was very helpful. Nov 15, 2021 at 20:19

Looking at the documentation at scipy, they recommend to use scipy.integrate.solve_ivp rather than scipy.integrate.odeint as can be seen here.

I changed hence the function and rearranged it a little bit.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# parameters
r1 = 1.0  #fixed
r2 = 0.5  #fixed
e1 = 1.0  #fixed
e2 = 0.75 #fixed

alpha12 = 1.5 #vary from 0-2
alpha21 = 1 #vary from 0-2

N1_0 = 1
N2_0 = 1

# store initial values in an array
X0 = [N1_0,N2_0]

# The two equations are contained in dX
def mutualism(t,X):
N1, N2 = X
dX = np.zeros(2)
dX[0] = N1 * (r1 - (e1 * N1) + (alpha12 * N2)) # equation for dN1/dt
dX[1] = N2 * (r2 - (e2 * N2) + (alpha21 * N1)) # equation for dN2/dt
return dX

# set time length
t = np.linspace(0,1.6,100)

X = solve_ivp(mutualism,t[[0,-1]],X0,t_eval=t)

N1 = X.y[0]
N2 = X.y[1]
t_fin = X.t

plt.plot(t_fin,N1, color='blue', lw=3)
plt.plot(t_fin,N2, color='red', lw=3)
plt.show()


This gives me the following plot. (Not sure if this is what you wanted since I dunnot really know the mutualism model.

As stated in the comments, we obtain a singularity at t~=1.5 which tends to infinity. Therefore the code only goes a little bit above that value.

• This worked, thank you! Nov 15, 2021 at 16:46
• I've just added your corrections @LutzLehmann with the updated plot. Nov 16, 2021 at 8:50