I'm trying to replicate the results of "Modeling Plasma Virus Concentration during Primary HIV Infection", by fitting the data sets to my ODE. Essencially I need to fit 4 parameters (d, k, delta, p) to a dif. equation given a data set.
Up until now I got the following code done (I'm a begginer):
c = 3 #global constant
v0= 10e-6
#dados
timeDataPatient1=np.array([0,22,43,78,106,146,183,230,268,358,435,489,519,534,584,610,687,778])
VirDataPatient1=10e3*np.array([v0*10e-3,27.2,210,85.9,81.1,46.2,60.1,82.8,103,72.1,79.4,70.4,207,42.6,10.8,54.2,22.3,40.8])
xdata = timeDataPatient1
ydata = VirDataPatient1
#inicial condition and max time of data set
IC = [10,0,10e-6]
t_max = max(xdata)
def diffeqs(t,y,d,k,delta,p):
'''Model of the differencial equation'''
global lam, c
T=y[0]
I=y[1]
V=y[2]
dT = 10*d - d*T-k*T*V
dI = k*T*V - delta*I
dV = p*I-c*V
return [dT,dI,dV]
def Vtime(theta):
'''solves the differencial equation, and returns the values of the model on the time points corresponding to the data set'''
global t_max, IC, xdata
sol = solve_ivp(diffeqs, (0,t_max), IC, args = (theta),t_eval=xdata)
return sol.y[2]
#now I define the objetive function to minimize. For a parameter theta, that corresponds to d,k,delta,p in the ODE model,it calculates the difference between the log of the data given and the predicts by the ODE model.
def calcErrorHIV(theta):
'''objetive function given in the paper'''
global ydata
dif = np.log(ydata)-np.log(Vtime(theta))
return sum(dif**2)
#should return theta parameters that best fits the data, but It doesn't compute.
sp.optimize.least_squares(calcErrorHIV,[0.013,0.46e-2,0.40,980])
but when I try to run the otimize or least_squares function I get no result.
Does my code have a problem or the 4 parameters problem takes a long time to solve?
Thanks for the help.
Edit: I tried the inicial guess [0.013,0.46 10*(-3),0.40,980]
Edit2: Added extra info about the problem.