Reproducing results from Goldbeter, Albert. "A model for circadian oscillations in the Drosophila period protein (PER)." Proc. R. Soc. Lond. B 261.1362 (1995): 319-324.

Goldbeter95Model.py

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


def PERModel(y, t, vd = 0.95):
    #here y0 is Pn , y1 is M, y2 is P0, y3 is P1, y4 is P2
    vs, vm = 0.76, 0.65
    Km, Ki, Kd, = 0.5, 1, 0.2
    ks = 0.38
    K1 = K2 = K3 = K4 = 2
    k1, k2 = 1.9, 1.3
    V1, V2, V3, V4 = 3.2, 1.58, 5, 2.5
    n = 4

    dMdt = vs * (Ki**n / (Ki**n + y[0]**n)) - vm*(y[1] / (Km + y[1]))
    dP0dt = ks * y[1] - V1 * (y[2] / (K1 + y[2])) + V2 * (y[3] / (K2 + y[3]))
    dP1dt = V1 * (y[2] / (K1 + y[2])) - V2 * (y[3] / (K2 + y[3])) - V3 * (y[3] / (K3 + y[3])) + V4 * (y[4] / (K4 + y[4]))
    dP2dt = V3 * (y[3] / (K3 + y[3])) - V4 * (y[4] / (K4 + y[4])) - k1 * y[4] + k2 * y[0] - vd * (y[4] / ([Kd + y[4]]))
    dPndt = k1 * y[4] - k2 * y[0]
    return [dPndt, dMdt, dP0dt, dP1dt, dP2dt]

y0 = [0.25, 1, 0.25, 0.25, 0.25]
t = np.linspace(0, 100, num=1000)
y = odeint(PERModel, y0, t, (0.95,))
plt.figure()
plt.plot(t, y[:,0])
plt.plot(t, y[:,1])
plt.plot(t, y[:,2])
plt.plot(t, y[:,3])
plt.show()

y0_central = [0.8, 1.5, 0.8, 0.8, 0.8]
z = odeint(PERModel, y0_central, t)

plt.figure()
plt.plot((y[:,0] + y[:,2] + y[:,3]),y[:, 1])
plt.plot((z[:,0] + z[:,2] + z[:,3]),z[:, 1])
plt.show()

vdRange = np.linspace(0.6, 1.5, num=15)
for i in range(len(vdRange)):
    y = odeint(PERModel, y0, t, (vdRange[i],))
    plt.figure()
    plt.plot(t, y[:, 0])
    plt.plot(t, y[:, 1])
    plt.plot(t, y[:, 2])
    plt.plot(t, y[:, 3])
    plt.show()
    plt.figure()
    plt.plot((y[:, 0] + y[:, 2] + y[:, 3]), y[:, 1])
    plt.show()

# Conclusion: as we increase the vd parameter - degradation rate - the period of the waves and thus
# the circadian cycle gets longer