Initial values used don't match those declared

[lhdp0110]

The graph produced by the program below shows the behavior of all the equations of my model once one of the equations is fitted to data points. However, the initial values on the graph do not match the initial values provided for some of the equations. For example, bP[0] appears to be 1 when it should be 0, same with bI[0], and NG1[0] appears to be 0 when it should be 1.

from symfit import variables, Parameter, Fit, D, ODEModel, exp
import numpy as np
import matplotlib.pyplot as plt

d_A = 0.5     
d_G = 0.02    
d_P = 0.01       
d_bP = 0.02     
d_bI = 0.01      
d_R = 0.01     
d_N = 0.01       
d_bN = 0.05   
d_NG1 = 0.05    
d_NG2 = 0.05    

s_A = d_A    
s_P = d_P    
s_B = d_R     
s_NG1 = d_NG1     
s_NG2 = d_NG1    

c_A = 0.8    
c_P = 0.5    
c_I = 1   
c_B = 4       
c_N = 3   
c_bN = 0.05      
c_NG1 = 0.1
c_NG2 = 0.2

K_I = 1
f_P =  9   
phi_NG1 = 0.07
phi_NG2 = 0.32
p_G = 2.5      
m_G = 2  

tdata = np.array([0, 2, 4, 8, 12, 24, 48])
concentration = np.array([1, 1.2, 1.8, 2.3, 2.1, 1.6, 1.5])

t, A, G, P, bP, bI, RB, RP, RN, bRN, NG1, NG2 = variables('t, A, G, P, bP, bI, RB, RP, RN, bRN, NG1, NG2')

c_A = Parameter('c_A', 0.8, min=0)

model = ODEModel({
    D(A, t): s_A + c_A * bRN - d_A * A,
    D(G, t): p_G * G - m_G * G * A - d_G * G,
    D(P, t): s_P - 2 * c_P * P**2 * G * exp(-phi_NG1 * NG1) - d_P * P,
    D(bP, t): c_P * P**2 * G * exp(-phi_NG1 * NG1) - d_bP * bP,
    D(bI, t): c_I * (1-bI/K_I) * bP * exp(-phi_NG2 * NG2) - d_bI * bI,  
    D(RB, t): s_B - c_B * RB * bI - d_R * RB,
    D(RP, t): c_B * RB * bI - f_P * RP - d_R * RP, 
    D(RN, t): f_P * RP - c_N * RN + c_bN * bRN - d_N * RN,
    D(bRN, t): c_N * RN - c_bN * bRN - d_bN * bRN,   
    D(NG1, t): s_NG1 + c_NG1 * bRN - d_NG1 * NG1,
    D(NG2, t): s_NG2 + c_NG2 * bRN - d_NG2 * NG2
    },
    initial={t: tdata[0], A: concentration[0], G: 1, P: 1, bP: 0, bI: 0, RB: 1, RP: 0, RN: 0,
             bRN: 0, NG1: 1, NG2: 1}
)

fit = Fit(model, t=tdata, A=concentration, G=None, P=None, bP=None, bI=None, RB=None, RP=None, RN=None,
          bRN=None, NG1=None, NG2=None)
fit_result = fit.execute()

print(fit_result)

taxis = np.linspace(0, 50)
A_fit, G_fit, P_fit, bP_fit, bI_fit, RB_fit, RP_fit, RN_fit, bRN_fit, NG1_fit, NG2_fit, = model(t=taxis, **fit_result.params)
plt.scatter(tdata, concentration)
plt.plot(taxis, A_fit, label='[A]')
plt.plot(taxis, G_fit, label='[G]')
plt.plot(taxis, P_fit, label='[P]')
plt.plot(taxis, bP_fit, label='[bP]')
plt.plot(taxis, bI_fit, label='[bI]')
plt.plot(taxis, RB_fit, label='[RB]')
plt.plot(taxis, RP_fit, label='[RP]')
plt.plot(taxis, RN_fit, label='[RN]')
plt.plot(taxis, bRN_fit, label='[bRN]')
plt.plot(taxis, NG1_fit, label='[NG1]')
plt.plot(taxis, NG2_fit, label='[NG2]')
plt.xlabel('Time in hours')
plt.ylabel('[X]')
plt.ylim(0, 3)
plt.xlim(0, 50)
plt.legend()
plt.show()

[pckroon]

It's working as intended :)

The problem is in this line: A_fit, G_fit, P_fit, bP_fit, bI_fit, RB_fit, RP_fit, RN_fit, bRN_fit, NG1_fit, NG2_fit, = model(t=taxis, **fit_result.params) The order in which you list the fit components is wrong.

result = model(t=taxis, **fit_result.params)
for compo in result.variables:
    print(f'{compo} = {result._asdict()[compo][0]}')

gives

A = 1.0
G = 1.0
NG1 = 1.0
NG2 = 1.0
P = 1.0
RB = 1.0
RN = 0.0
RP = 0.0
bI = 0.0
bP = 0.0
bRN = 0.0

as expected.

The components in a model's output are always sorted alphabetically. See also the last bit of https://symfit.readthedocs.io/en/stable/tutorial.html#named-models

[lhdp0110]

I did not realize that, thank you!

[pckroon]

No problem :) Happy fitting!