%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
import seaborn as sns
from rakott.mpl import legend_out
def ribeck_ode(x, t, K, r, y, h):
N1, N2, R1, R2 = x
r11, r12, r21, r22 = r
K11, K12, K21, K22 = K
y11, y12, y21, y22 = y
h1, h2 = h
dN1dt = (
N1 * r11 * R1 / (K11 + R1) +
N1 * r21 * R2 / (K21 + R2)
)
dN2dt = (
N2 * r12 * R1 / (K12 + R1) +
N2 * r22 * R2 / (K22 + R2)
)
dR1dt = (
-1/y11 * N1 * r11 * R1 / (K11 + R1)
-1/y12 * N2 * r12 * R1 / (K12 + R1)
)
dR2dt = (
h1 * N1 * r11 * R1 / (K11 + R1) +
h2 * N2 * r12 * R1 / (K12 + R1)
-1/y21 * N1 * r21 * R2 / (K21 + R2)
- 1/y22 * N2 * r22 * R2 / (K22 + R2)
)
return [dN1dt, dN2dt, dR1dt, dR2dt]
r11 = 1
r12 = 2.68
r21 = 7
r22 = 1
r = r11, r12, r21, r22
K11 = K12 = 50
K21 = K22 = 50
K = K11, K12, K21, K22
y11 = y12 = 1
y21 = y22 = 1
y = y11, y12, y21, y22
h1 = h2 = 1
h = h1, h2
R1_0 = 1
x0 = 0.4, 0.6, R1_0, 0
t = np.linspace(0, 50, 100)
x = odeint(ribeck_ode, x0, t, args=(K, r, y, h))
N1, N2, R1, R2 = x.T
red, blue, green, purple, orange = sns.color_palette('Set1', 5)
fig, axes = plt.subplots(2, 1, figsize=(5, 8))
ax = axes[0]
linekws = dict(lw=2)
ax.plot(t, N1, color=green, label='$N_1$', **linekws)
ax.plot(t, N2, color=red, label='$N_2$', **linekws)
# ax.plot(t, N1+N2, color=blue, label='$N_1+N_2$', **linekws)
ax.plot(t, R1, color=purple, label='$R_1$', alpha=0.75, **linekws)
ax.plot(t, R2, color=orange, label='$R_2$', alpha=0.75, **linekws)
ax.set(xlabel=('Time'), ylabel=('Density'))
legend_out(ax)
ax = axes[1]
ax.plot(t, N1/(N1+N2), color=green, label='$N_1$', **linekws)
ax.plot(t, N2/(N1+N2), color=red, label='$N_2$', **linekws)
ax.set(xlabel=('Time'), ylabel=('Frequency'), ylim=(0,1))
fig.tight_layout()
sns.despine()
w = np.linspace(0, 10)
plt.plot(w, 2*w/(1+w))
plt.plot([0,1,2],[0,1,2], 'ok')
plt.ylim(0,2)
c = 0.1
plt.axhline(2*(1-c), ls='--', color='k')
sns.despine()