scipy.integrate odeint

[Jose Guzman]

Hello everybody!

I am trying to solve a basic 2-state kinetic scheme of the form:

# Two-state kinetic model
#
#        alpha*[cc]
# (1-n) ----------> n
#       <----------
#          beta


where [cc] is a function step function. To solve it numerically, I used scipy.integrate.odeint. Everything works fine if I define a step function that takes the value 1 between 2 and 3 .


from numpy import linspace
from scipy.integrate import odeint

# define a square pulse with the concentration of Glu

def pulse(t):
    """ t in ms """
    if t<=2 or t>=3:
        return 0.0
    else:
        return 1.0 # in mM
   
# independent variable
time = linspace(0,10, 1000) # 10 ms


# define differential equation
def diff(n, t, alpha, beta):
    """
    alpha in 1/(mM*ms)
    beta in 1/ms
    dn/dt = alpha*cc*(1-n) - beta*n
    """
    cc = pulse(t) # square pulse
    return alpha*cc*(1-n)-beta*n   

# and solve it for alpha = 1 and beta = 1
y = odeint(func = diff, y0=0, t= time, args=(1,1))


However, if i change the pulse function to return 1 between, say 5 and 6 the solution to the differential equation returns 0 for all the interval. I cannot figure out why.

# new function
def pulse(t):
    """ t in ms """
    if t<=5 or t>=6:
        return 0.0
    else:
        return 1.0 # in mM
y = odeint(func = diff, y0=0, t= time, args=(1,1)) # returns all zeros

I would greatly appreciate some help here

Best

Jose