3
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ @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. $\endgroup$ Nov 15, 2021 at 19:30
  • 2
    $\begingroup$ 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. $\endgroup$ Nov 15, 2021 at 20:07
  • $\begingroup$ @LutzLehmann I see! I was wondering why it seemed to bug after a certain point. Thank you for the insight, it was very helpful. $\endgroup$ Nov 15, 2021 at 20:19

1 Answer 1

5
$\begingroup$

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.

enter image description here

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.

$\endgroup$
1
  • 2
    $\begingroup$ I've just added your corrections @LutzLehmann with the updated plot. $\endgroup$
    – coyote
    Nov 16, 2021 at 8:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.