Curve fitting of Monod growth/degradation equations to the experimental data

Question

So the problem that is being faced here is the curve fitting of the Monod equations to the experimental data. The model of bacteria growth and degradation of the organic carbon looks like this:

dX/dt = (u * S * X )/(K + S)

dS/dt = ((-1/Y) * u * S * X )/(K + S)

These equations are solved using the scipy odeint function. Results after integration are stored into two vectors, one for growth, and the another one for degradation. The next step is to curve fit this model to the experimentally observed data and estimate the model parameters: u, K and Y. Once the code is run, the following error is produced:

File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\minpack.py", line 392, in leastsq
    raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))

TypeError: Improper input: N=3 must not exceed M=2"

For the convenience, curve fitting part is commented out, so the plot of the expected result can be generated. Bellow is the code sample:

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

"""Experimental data!"""
t_exp = np.array([0, 8, 24, 32, 48, 96, 168])
S_exp = np.array([5.5, 4.7, 3.7, 2.5, 1.5, 0.7, 0.5])
X_exp = np.array([10000, 17000, 30000, 40000, 60000, 76000, 80000])

"Model of the microbial growth and the TOC degradation"
# SETTING UP THE MODEL
def f(t, u, K, Y):
     'Function that returns mutually dependent variables X and S'
     def growth(x, t):
         X = x[0]
         S = x[1]
         "Now differential equations are defined!"
         dXdt = (u * S * X )/(K + S)
         dSdt = ((-1/Y) * u * S * X )/(K + S)
         return [dXdt, dSdt]
     # INTEGRATING THE DIFFERENTIAL EQUATIONS
     "initial Conditions"
     init = [10000, 5]
     results = odeint(growth, init, t)
     "Taking out desired column vectors from results array"
     return results[:,0], results[:,1]

# CURVE FITTING AND PARAMETER ESTIMATION
"""k, kcov = curve_fit(f, t_exp, [X_exp, S_exp], p0=(1, 2, 2))
u = k[0]
K = k[1]
Y = k[2]"""

# RESULTS OF THE MODEL WITH THE ESTIMATED MODEL PARAMETERS
t_mod = np.linspace(0, 168, 100)
compute = f(t_mod, 0.8, 75, 13700)# these fit quite well, but estimated manually
X_mod = compute[0]
S_mod = compute[1]

# PLOT OF THE MODEL AND THE OBSERVED DATA
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.plot(t_exp, X_exp, "yo")
ax1.plot(t_mod, X_mod, "g--", linewidth=3)
ax1.set_ylabel("X")

ax2 = ax1.twinx()
ax2.plot(t_exp, S_exp, "mo", )
ax2.plot(t_mod, S_mod, "r--", linewidth=3)
ax2.set_ylabel("S", color="r")
for tl in ax2.get_yticklabels():
    tl.set_color("r")
plt.show()

Any advice of how to cope with this problem and proceed further would be highly appreciated. Thanks in advance.

Answer 1

The result of f() needs to have the same shape as the experimental data you feed into curve_fit as third parameter. In the last line of f() you just take the t = 0s value of the solution for both ODEs and return that, but you should return the complete solution. When fitting several sets of data at once using curve_fit, just concat them (stack horizontally), i.e.

def f(t, u, K, Y):
   .....
   return np.hstack((results[:,0], results[:,1]))

and call curve_fit like

k, kcov = curve_fit(f, t_exp, np.hstack([X_exp, S_exp]), p0=(1, 2, 2))

You will have to adapt the plotting part of your script, too:

compute = f(t_mod, u, K, Y)
compute = compute.reshape((2,-1))