Let's start with a simple model, that should contain the main mechanisms involved:
$$ \begin{aligned} \frac{dH}{dt} &= H (\rho(H) - \alpha D) - \frac{(1-q)H}{1+qahD + (1-q)a h_P H} aM \\ \frac{dM}{dt} &= \frac{qD + (1-q)H}{1+qahD + (1-q)a h_P H} aM - (\mu + g(D)) M \\ \frac{dD}{dt} &= i_D H- sD - \frac{qD}{1+qahD + (1-q)a h_P H} aM \end{aligned}$$The host population (or biomass) $H$ grows with a per capita birth rate of $\rho$, which decreases in the presence of detritus ($D$). It is also affected by parasitism. Mutualists ($M$) feed on detritus, and have a constant mortality rate besides a death rate imposed by host cleaning ($c$), which depends on the amount of mutualists. They can also feed on the host itself, and the parameter $q$ determines the preference for detritus or host, so that $q=1$ corresponds to a "pure" mutualist that never attacks the host. We assume no other trade-offs between "benign" mutualists and those that parasitize the host as well- in particular, the host rubs them off with the same efficiency. Finally, detritus are continually depositing in the hosts, in a rate proportional to the number of hosts, while being consumed by the mutualists. The detritus may (or may not) have a rate of decomposition $s$.
We can initially choose $\rho = 1-H/K$ (logistic growth) and $g(D) = \beta/D$. With these considerations, the basic model is in place, and we can ask the following questions:
Necessary conditions for coexistence:
$$ \begin{aligned} \mu h &< 1 \\ \frac{\alpha}{r} + s &< 1 \\ \frac{\alpha i_D}{r s} &<1 \end{aligned}$$
In [13]:
%matplotlib inline
from scipy.integrate import ode
import numpy as np
from matplotlib import pyplot
def model(t, x, a, q, h, h_p, alpha, beta, mu, s, i_d, k):
xt = np.zeros(3)
rho = 1 - x[0]/k
g = beta/x[2]
theta = 1 + q*a*h*x[2] + (1 - q)*a*h_p*x[0]
xt[0] = x[0]*(rho - alpha*x[2]) - (1 - q)*(x[0]/theta)*a*x[1]
xt[1] = x[1]*a*(q*x[2] + (1 - q)*x[0])/theta - (mu + g)*x[1]
xt[2] = i_d*x[0] - s*x[2] - a*x[1]*q*x[2]/theta
return xt
def mutate_model(t, x, a, q, q_m, h, h_p, alpha, beta, mu, s, i_d, k):
xt = np.zeros(4)
rho = 1 - x[0]/k
g = beta/x[2]
theta = 1 + q*a*h*x[2] + (1 - q)*a*h_p*x[0]
theta_m = 1 + q_m*a*h*x[2] + (1 - q_m)*a*h_p*x[0]
xt[0] = x[0]*(rho - alpha*x[2]) - (1 - q)*(x[0]/theta)*a*x[1] - (1 - q_m)*(x[0]/theta_m)*a*x[3]
xt[1] = x[1]*a*(q*x[2] + (1 - q)*x[0])/theta - (mu + g)*x[1]
xt[2] = i_d*x[0] - s*x[2] - a*x[1]*q*x[2]/theta - a*x[3]*q_m*x[2]/theta_m
xt[3] = x[3]*a*(q_m*x[2] + (1 - q_m)*x[0])/theta_m - (mu + g)*x[3]
return xt
In [14]:
def integrate(func, x0, t0, tf, dt, *params):
r = ode(func).set_integrator('lsoda')
r.set_initial_value(x0, t0).set_f_params(*params)
series = []
while r.successful() and r.t < tf:
r.integrate(r.t + dt)
series.append(r.y)
return np.array(series)
In [30]:
xo = np.array([0.6, 0.3, 0.2])
# a, q, h, h_p, alpha, beta, mu, s, i_d, k
a = 1
q = 1
h = 0.9
h_p = 1
alpha = 1
beta = 0.1
mu = 0.1
s = .01
i_d = 1
k = 1
x = integrate(model, xo, 0, 40, 0.01, a, q, h, h_p, alpha, beta, mu, s, i_d, k)
In [31]:
pyplot.plot(x[:, 0], label='H')
pyplot.plot(x[:, 1], label='M')
pyplot.plot(x[:, 2], label='D')
pyplot.legend()
pyplot.xlabel('Time')
pyplot.ylabel('Population')
pyplot.savefig('equ.png', dpi=600)
In [ ]:
In [ ]:
In [ ]:
In [105]:
xo = np.array([0.7, 0.5, 0.2, 0.1])
# a, q, h, h_p, alpha, beta, mu, s, i_d, k
a = 1.
q = 1
q_m = 0.95
h = .5
h_p = 1.1
alpha = 1
beta = 0.1
mu = 0.1
s = .2
i_d = 1
k = 1
x = integrate(mutate_model, xo, 0, 20000, 1, a, q, q_m, h, h_p, alpha, beta, mu, s, i_d, k)
In [106]:
pyplot.plot(x[:, 0], label='$H$')
pyplot.plot(x[:, 1], label='$M$')
pyplot.plot(x[:, 2], label='$D$')
pyplot.plot(x[:, 3], label='$M\'$')
pyplot.legend()
pyplot.xlabel('Time')
pyplot.ylabel('Population')
pyplot.savefig('coexist.png', dpi=600)
In [40]:
import IPython.html.widgets as widgets
def plot_interactive(**kwargs):
x0 = np.array([kwargs['h0'], kwargs['m0'], kwargs['d0'], kwargs['m_m0']])
a = 1
q = 1
q_m = 0.1
h = 0.5
alpha = 1
beta = 0.1
mu = 0.1
s = .01
i_d = 1
k = 1
t0 = 0
tf = 20
dt = .1
x = integrate(mutate_model, x0, t0, tf, dt, a, kwargs['q'], kwargs['q_m'], h, kwargs['h_p'], alpha, beta, mu, kwargs['s'], i_d, k)
pyplot.plot(x[:, 0], label='H')
pyplot.plot(x[:, 1], label='M')
pyplot.plot(x[:, 2], label='D')
pyplot.plot(x[:, 3], label='M_m')
pyplot.legend()
In [41]:
widgets.interact(plot_interactive,
h0 = widgets.FloatSliderWidget(description="H", min=0, max=1, step=0.05, value=0.6),
m0 = widgets.FloatSliderWidget(min=0, max=1, step=0.05, value=0.5),
d0 = widgets.FloatSliderWidget(min=0, max=1, step=0.05, value=0.2),
m_m0 = widgets.FloatSliderWidget(min=0, max=1, step=0.05, value=0.2),
s = widgets.FloatSliderWidget(min=0, max=1, step=0.05, value=.1),
q = widgets.FloatSliderWidget(min=0, max=1, step=0.05, value=1.0),
q_m = widgets.FloatSliderWidget(min=0, max=1, step=0.05, value=0.9),
h_p = widgets.FloatSliderWidget(description="h_P", min=0.1, max=10., step=0.1, value=1.1))
In [ ]: