I am using first-order differential equations to model viral spread and am using a likelihood approach to fit the ODE to experimental data. However, whenever the code is run, the fitted initial value of the ODE changes each time the code is run, while all of the other parameters in the model converge to the same value every time. My code is shown below.
import numpy as np
from scipy.integrate import odeint
from lmfit import minimize, Parameters, Parameter, report_fit
import matplotlib.pyplot as plt
#====================================================================
''' Data on which the model will be fitted '''
TROMAS_DATA = [[3,3,1,7219,33],[3,3,2,7378,102],[3,3,3,8978,66],[3,3,4,8105,116],[3,3,5,11941,70],
[5,3,1,12349,214],[5,3,2,8958,154,],[5,3,3,10661,156],[5,3,4,11924,305],[5,3,5,9848,417],
[7,3,1,8854,429],[7,3,2,11098,886],[7,3,3,12757,470],[7,3,4,10233,594],[7,3,5,8347,721],
[10,3,1,8347,721],[10,3,2,8895,586],[10,3,3,9554,999],[10,3,4,11021,940],[10,3,5,9416,917]]
#====================================================================
''' Parameter list '''
likeParams = Parameters()
likeParams.add('I0', value = .00372, min = .0000001, max = 1.0000)
likeParams.add('b', value = .871, min = .0001, max = 1.0000)
likeParams.add('psi3', value = .083, min = .0001, max = 1.0000)
#====================================================================
''' Initial conditions '''
I0 = [likeParams['I0'].value]
#====================================================================
''' Time axis '''
ZERO_DAYS_AXIS = [0, 3, 5, 7, 10]
DAYS_AXIS = [3, 5, 7, 10]
t = np.linspace(0, 10, 10000)
#====================================================================
''' Define the model '''
def model(M3, t, parameters):
try: # Get parameters
b = parameters['b'].value
psi3 = parameters['psi3'].value
except KeyError:
b, psi3 = parameters
if (M3 < psi3):
S3 = (1 - (M3 / psi3))
else:
S3 = 0
dM3dt = b * M3 * S3
return dM3dt
#====================================================================
''' Compute negative log likelihood of Tromas' data given the model '''
def negLogLike(parameters):
# Solve ODE system to get model values; parameters are not yet fitted
M3 = odeint(model, I0, ZERO_DAYS_AXIS, args=(parameters,))
nll = 0
epsilon = 10**-10
for t in range(4): # Iterate through days
for p in range(5): # Iterate through replicates
Vktp = TROMAS_DATA[5 * t + p][4] # Number of infected cells
Aktp = TROMAS_DATA[5 * t + p][3] + Vktp # Total number of cells observed
Iktp = M3[t + 1] # Frequency of cellular infection
if (Iktp <= 0):
Iktp = epsilon
elif (Iktp >= 1):
Iktp = 1 - epsilon
nll += Vktp * np.log(Iktp) + (Aktp - Vktp) * np.log(1 - Iktp)
return -nll
#====================================================================
''' Miminize negative log likelihood with differential evolution algorithm '''
result_likelihood = minimize(negLogLike, likeParams, method = 'differential_evolution')
#====================================================================
''' Compare previous and fitted values '''
print("Original I0 = .00372, My I0 = ", result_likelihood.params['I0'].value)
print("Original beta = .871, My beta = ", result_likelihood.params['b'].value)
print("Original psi3 = .083, My psi3 = ", result_likelihood.params['psi3'].value)
In order to help illustrate the problem, the output of three runs of the code is shown below. I also know that the "correct" value for I0 should be roughly around 0.003.
1st time:
Original I0 = .00372, My I0 = 0.9711066426171252
Original beta = .871, My beta = 0.44374676021094367
Original psi3 = .083, My psi3 = 0.11330842894454747
2nd time:
Original I0 = .00372, My I0 = 0.09183577426999114
Original beta = .871, My beta = 0.4437477735310379
Original psi3 = .083, My psi3 = 0.11330747932246728
3rd time:
Original I0 = .00372, My I0 = 0.47857244584676956
Original beta = .871, My beta = 0.44374801498547956
Original psi3 = .083, My psi3 = 0.11330824660617532
If you could help me learn how to correctly fit the initial condition I would greatly appreciate it. Thank you!
Your model is slightly hard for me to follow. I wonder why define
I0 = [likeParams['I0'].value]
and then use that within your negLogLike function and do not actually use the varying value for I0 in your fitting model. Perhaps that should be
M3 = odeint(model, parameters['I0'].value, ZERO_DAYS_AXIS, args=(parameters,))
or some variation of that? That seems to give a more satisfying result. I tried it with method='Nelder' to get:
Original I0 = .00372, My I0 = 0.0010018942980410726
Original beta = .871, My beta = 0.7074140684730617
Original psi3 = .083, My psi3 = 0.08807904093119709