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.