Inference of Parameters of an Ordinary Differential Equation using SciPy and nlopt

gistfile1.py

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import leastsq
import nlopt

def michelis_menten(y, t, *args):
    Vmax = args[0][0]
    km = args[0][1]
    St = args[0][2]
    P = y[0]
    S = St - P
    
    dP = Vmax * (S / (S+km))
    return dP
 
def mm_residuals(params, *args):

	t_exp = args[0]
	y_exp = args[1]
	# Experimental noisy data

	P_0 = 0
	# Initial Condition is 0

	params = tuple(params)
   # SciPy wants the params in a tuple

	t_sim = np.linspace(0, 100, 100)
	y_sim = odeint(michelis_menten, [P_0], t_sim, args = (params,))
	y_sim = y_sim.flatten()
	mapped_t = np.searchsorted(t_sim, t_exp)
	# Maps the experimental timepoints to the closest simulated timepoint

	return (y_sim[mapped_t] - y_exp)
 
def lmopt_residuals(x, grad):
    
    res = mm_residuals(x, *(noisy_t, noisy_P))
    res = np.sum(res**2)
    return res
 

# Parameters MM
Vmax = 1
km = 3
St = 10
mm_params = (Vmax, km, St)

# Initial Conditions MM
P_0 = 0

# Timesteps
n_steps = 1000
t = np.linspace(0, 100, n_steps)

num_P = odeint(michelis_menten, [P_0], t, args = (mm_params,))

noisy_P = num_P + np.random.randn(n_steps, 1)
noisy_P[noisy_P < 0] = 0
noisy_P = noisy_P[::10].flatten()
noisy_t = t[::10]

# Michelis-Menten Least Squares Using SciPy
mm_init_guess = (0, 0, 12)
out = leastsq(mm_residuals, mm_init_guess, args = (noisy_t, noisy_P),
              full_output = 1)
mm_lsq_params = out[0].flatten()
mm_lsq_P = odeint(michelis_menten, P_0, t,
                  args = (tuple(mm_lsq_params),)).flatten()
min_ls = np.sum(mm_residuals(mm_lsq_params, *(noisy_t, noisy_P))**2)

# Michelis-Menten Least Squares Using nlopt
opt = nlopt.opt(nlopt.GN_DIRECT_L, 3)
opt.set_lower_bounds([0.0, 0.0, 0.0])
opt.set_upper_bounds([10.0, 5.0, 15.0])
opt.set_min_objective(lmopt_residuals)
opt.set_maxtime(10)
x = opt.optimize(np.array([0.1, 0.1, 12]))
minf = opt.last_optimum_value()

print "NL OPT:"
print "optimum at ", x
print "minimum value = ", minf
print "result code = ", opt.last_optimize_result()

print "\nLS:"
print "optimum at ", mm_lsq_params
print "minimum value = ", min_ls