0

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.

Daniel Fonseca
  • 1
  • 2
  • Can you try t modify the code to make it a working example, please. – mikuszefski Dec 15 '22 at 15:15
  • @mikuszefski I will update it but basically I want to minimize calcErrorHIV, given the inicial guess [0.013,0.46 10*(-3),0.40,980]. – Daniel Fonseca Dec 15 '22 at 16:28
  • @mikuszefski Added extra comments to the code – Daniel Fonseca Dec 15 '22 at 17:23
  • I think the solver works, although there are some parameter sets that take very long to solve. The error function should only return `diff`. The sum and square is done by `least_squares`. I would play with some manual solutions to check the behavior first. Plot the result vs experimental data. Here I also would remove the`t_eval` and check later if the fixed time is actually sufficient to get a decent ode solution. Initial conditions are fixed? Later, when running the optimization, you can add a `print` to get a progress update. – mikuszefski Dec 16 '22 at 08:31
  • Concerning the `t_eval` I am not 100% sure what it is doing. If it is calculating additional points in-between, but returning the values only for the given `t`, it would be ok. – mikuszefski Dec 16 '22 at 08:34
  • Parameters can also be related to the saturation behavior. If the equations converge to a steady state you can relate this to the parameters by setting `dT=0`, `dI=0`, `dV=0` and solve manually. – mikuszefski Dec 16 '22 at 09:09
  • BTW, the `t_eval` is ok. – mikuszefski Dec 16 '22 at 09:10
  • ...and please put the correct imports and remove unnecessary variables...hence, minimal working example – mikuszefski Dec 16 '22 at 09:13

0 Answers0