Example: Multiple species Reaction Kinetics using ODEModel

In this example we shall fit to a complex system of ODEs, based on that published by Polgar et al. However, we shall be generating some mock data instead of using the real deal.

from symfit import (
    variables, parameters, ODEModel, D, Fit
)
from symfit.core.support import key2str
import numpy as np
import matplotlib.pyplot as plt

First we build a model representing the system of equations.

t, F, MM, FMM, FMMF = variables('t, F, MM, FMM, FMMF')
k1_f, k1_r, k2_f, k2_r = parameters('k1_f, k1_r, k2_f, k2_r')

MM_0 = 10  # Some made up initial amount of [FF]

model_dict = {
    D(F, t): - k1_f * MM * F + k1_r * FMM - k2_f * FMM * F + k2_r * FMMF,
    D(FMM, t): k1_f * MM * F - k1_r * FMM - k2_r * FMM * F + k2_f * FMMF,
    D(FMMF, t): k2_f * FMM * F - k2_r * FMMF,
    D(MM, t): - k1_f * MM * F + k1_r * FMM,
}
model = ODEModel(
    model_dict,
    initial={t: 0.0, MM: MM_0, FMM: 0.0, FMMF: 0.0, F: 2 * MM_0}
)
print(model)

Derivative(F, t; k1_f, k1_r, k2_f, k2_r) = -k1_f*F*MM + k1_r*FMM - k2_f*F*FMM + k2_r*FMMF
Derivative(FMM, t; k1_f, k1_r, k2_f, k2_r) = k1_f*F*MM - k1_r*FMM + k2_f*FMMF - k2_r*F*FMM
Derivative(FMMF, t; k1_f, k1_r, k2_f, k2_r) = k2_f*F*FMM - k2_r*FMMF
Derivative(MM, t; k1_f, k1_r, k2_f, k2_r) = -k1_f*F*MM + k1_r*FMM

Generate mock data.

tdata = np.linspace(0, 3, 20)
data = model(t=tdata, k1_f=0.1, k1_r=0.2, k2_f=0.3, k2_r=0.3)._asdict()
sigma_data = 0.3
np.random.seed(42)
for var in data:
    data[var] += np.random.normal(0, sigma_data, size=len(tdata))

plt.scatter(tdata, data[MM], label='[MM]')
plt.scatter(tdata, data[FMM], label='[FMM]')
plt.scatter(tdata, data[FMMF], label='[FMMF]')
plt.scatter(tdata, data[F], label='[F]')
plt.xlabel(r'$t$')
plt.ylabel(r'$C$')
plt.legend()
plt.show()

Perform the fit. Let’s pretend that for experimental reasons, we can only measure the concentration for MM and F, but not for the intermediate FMM nor the product FMMF. This is no problem, as we can tell symfit to ignore those components by setting the data for them to None.

k1_f.min, k1_f.max = 0, 1
k1_r.min, k1_r.max = 0, 1
k2_f.min, k2_f.max = 0, 1
k2_r.min, k2_r.max = 0, 1

fit = Fit(model, t=tdata, MM=data[MM], F=data[F],
          FMMF=None, FMM=None,
          sigma_F=sigma_data, sigma_MM=sigma_data)
fit_result = fit.execute()
print(fit_result)

Parameter Value        Standard Deviation
k1_f      9.540422e-02 4.440701e-03
k1_r      1.065140e-01 7.165730e-02
k2_f      2.706137e-01 5.305089e-02
k2_r      2.633632e-01 5.647271e-02
Status message         CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
Number of iterations   30
Objective              <symfit.core.objectives.LeastSquares object at 0xffff62aa3ef0>
Minimizer              <symfit.core.minimizers.LBFGSB object at 0xffff8d659280>

Goodness of fit qualifiers:
chi_squared            33.985494538721184
objective_value        16.992747269360592
r_squared              0.9936568366393397

taxis = np.linspace(tdata.min(), tdata.max(), 1000)
model_fit = model(t=taxis, **fit_result.params)._asdict()
for var in data:
    plt.scatter(tdata, data[var], label='[{}]'.format(var.name))
    plt.plot(taxis, model_fit[var], label='[{}]'.format(var.name))
plt.legend()
plt.show()

We see that the lack of data for some components is not a problem, they are predicted quite nicely.