ODE Fitting with Symfit for Python: How to get estimations for intial values?

Question

I would like to know how to extract the intial values for the state variables.

Basically, the intial values are also considered as parameters.

I'm providing the initial values because they are needed for integration. However, these values are sometimes considered as additional parameters (next to the model parameters), so that the initial values provided are also considered are initial guesses. In fact, when doing growth experiments, we have some initial conditions, but we may not know all of them (depending on the specific experimental conditions and system under investigation).

Consider a simple microbial growth model with the growth rate mu governed by the well-known Monod kinetic (parameters mumax and ks) and a constant biomass-substrate yield (parameter yxs) and a growth-associated product formation (parameter ypx). See the code below.

For further implementations I would like to calculate first-order sensitivities over time, parameter correlations, etc.

from symfit import variables, parameters, ODEModel, Fit, D
import numpy as np
from matplotlib import pyplot as plt
# define ODE model
ks, mumax, ypx, yxs = parameters('ks, mumax, ypx, yxs')
p, s, x, t = variables('p, s, x, t')
model_dict = {
    D(x, t): mumax * s / (ks + s) * x,
    D(s, t): -1/yxs * mumax * s / (ks + s) * x,
    D(p, t): ypx * mumax * s / (ks + s) * x,
}
ode_model = ODEModel(model_dict, initial={t:0.0, x:0.1, s:20.0, p:0.0})

# Generate noisy measurement data
tSample = np.array([1, 2, 4, 6, 8, 12])
pSample, sSample, xSample = ode_model(t=tSample, ks=0.01, mumax=0.4, ypx=0.2, yxs=0.5)
xRelErr = 0.05
sRelErr = 0.10
pRelErr = 0.10
xNoise = xSample + xSample * xRelErr * np.random.randn(xSample.size)
sNoise = sSample + sSample * sRelErr * np.random.randn(sSample.size)
pNoise = pSample + pSample * pRelErr * np.random.randn(pSample.size)

# constraints for parameters
ks.min = 0.0
mumax.min = 0.0
ypx.min = 0.0
yxs.min = 0.0

# initial guesses for parameters
ks.value = 0.01
mumax.value = 0.4
ypx.value = 0.2
yxs.value = 0.5

# perform fitting
fit = Fit(ode_model, t=tSample, x=xNoise, s=sNoise, p=pNoise)
fit_result = fit.execute()

# show fit
print(fit_result)
print(fit_result.covariance_matrix)

Answer 1

I can confirm that this suggestion works in Symfit 0.5.1, with the caveat that the "initial" value needs to be specified as, e.g., x0*.value*. This creates an odd issue when you call "print(fit_result)", it complains that "'NoneType' object has no attribute 'sqrt'", but everything else seems fine.

Here's a (more or less) minimal example to demonstrate this:

from symfit.core.minimizers import DifferentialEvolution, BFGS, BasinHopping
import symfit as sf
import numpy as np
import matplotlib.pyplot as plt

signal = np.array([ 3.69798582e-04, -4.42968073e-04, -1.62901574e-04, -2.66074317e-03,
       -2.69579455e-03,  7.74354388e-04,  2.19358448e-03,  6.38607419e-03,
        3.78642692e-03, -3.27400548e-03, -4.34699571e-03,  2.41753615e-04,
        9.96158926e-03,  1.89990480e-02,  3.60202254e-02,  7.63199623e-02,
        1.43007912e-01,  2.48988030e-01,  3.95977066e-01,  5.65958061e-01,
        7.12988678e-01,  8.09083671e-01,  8.50628154e-01,  8.72826358e-01,
        8.83509851e-01,  8.81799211e-01,  8.81649458e-01,  8.88619439e-01,
        8.82045565e-01,  8.91226071e-01,  9.07541872e-01,  9.12674298e-01,
        9.06886981e-01,  8.98514353e-01,  9.02275457e-01,  9.07204574e-01,
        9.01857725e-01,  8.92016771e-01,  8.99608399e-01,  8.98652789e-01,
        8.98421894e-01,  8.94907751e-01,  9.09122059e-01,  9.13065656e-01,
        9.09254937e-01,  9.10644851e-01,  8.99605570e-01,  8.99581621e-01,
        9.12054065e-01,  9.06361431e-01,  9.03838466e-01,  9.09130042e-01,
        9.13443979e-01,  9.18176457e-01,  9.08160892e-01,  9.03154574e-01,
        9.03069088e-01,  9.02597396e-01,  8.98854582e-01,  9.07801262e-01])
cycles = np.arange(len(signal))

s,c = sf.variables('s,c')
r,K,s0 = sf.parameters('r,K,s0')
lg2_e = np.log2(np.exp(1))

s0.min = 1e-7
s0.value = 1e-6
s0.max = 1e-5
K.min = 0.1
K.value = 1
K.max = 2
r.min = 0.1/lg2_e
r.value = 1/lg2_e
r.max = 2/lg2_e


logistic_eqn = {
    sf.D(s,c): r*s*(1-s/K) * (sf.cos(s0)**2 + sf.sin(s0)**2)
}

logistic_model = sf.ODEModel(logistic_eqn, initial = {c:0.0,s:s0.value})

logistic_fit = sf.Fit(logistic_model,c=cycles,s=signal)#, minimizer=DifferentialEvolution)
fit_result = logistic_fit.execute()
#print(fit_result)
sig_fit, = logistic_model(c=cycles, **fit_result.params)
plt.plot(cycles, signal, 'o')
plt.plot(cycles, sig_fit)

Data points with fitted curve

Answer 2

Great question. The short answer is no, because this has not been implemented yet. (Though it has been on the list for a while, see Issue #53.)

However, in some cases there might already be a workaround. Lets say you want

x0 = Parameter('x0')
ode_model = ODEModel(model_dict, initial={t: 0.0, x: x0, s: 20.0, p: 0.0})

in your example above. symfit will currently only check model_dict for Parameter objects, so it will not see x0. So if you model explicitly depends on x0 it will work. Perhaps you can find some way to make x0 appear in your model. Worst case scenario, you could multiply one of the components with some trivial identity like cos(x0)**2 + sin(x0)**2:

model_dict = {
    D(x, t): mumax * s / (ks + s) * x * (cos(x0)**2 + sin(x0)**2),
    D(s, t): -1/yxs * mumax * s / (ks + s) * x,
    D(p, t): ypx * mumax * s / (ks + s) * x,
}

However, this is provided without warranty ;). Something I noticed in my own research when I had a model who's model_dict did explicitly depend on the initial parameter was that symfit had difficulty estimating the Jacobian for such parameters, so fitting was unstable. However, this was a couple of versions ago and many things have changed since then so give it a try!

You can also use a non-gradient based method such as Nelder-Mead or hopefully very soon Differential Evolution.

I will look into this, I opened an issue for this here. If you'd like to get involved in this any help would definitely be welcome, for example when it comes to providing a minimal working example that I could use as a test.