Fitting a system of ODEs to data using lsqcurvefit to estimate system's parameters

Hi

I have a problem in using lsqcurvefit to fit data. I am not very good in matlab. My problem is about paramater estimation of a system of ODEs. I tried to use the lsqcurvefit in a way similar to fminsearchbnd. I managed to get a good fit and acceptable estimation of the parameters when using fminsearchbnd. However, when using lsqcurvefit I ended up with errors! One of the error tha I couldn't fix is: "Error using lsqcurvefit (line 229) FUN must have two input arguments"

I have attached the code and the excel sheet for data.

Thanks in advance for your help.

[Alex Sha]

In your xlsx file, the 2nd column data named "new cases", does this data corresponded to " I = y(5);" in your code? if so, the I0 = 4411, is this value correct? since this value is too bigger than all data values in "new cases" column, looks strange.

[Maryam Alyahyai]

Yes, you are right. I forgot to update the initial conditions! I0 = 86, R0= 233, S0=5215663 and the remaining are the same since they are only assumed. (These initials are based on total population N0=5223376).

Thanks for you reply @Alex Sha

[Alex Sha]

if taking "t=1, S = 5215663,E = 3779,A = 2103,P = 1512,I = 86,R = 233,N = 5215663+3779+2103+1512+86+233" as initial conditions, the result obtained by using 1stOpt is:

Sum Squared Error (SSE): 4596805.3510258
Root of Mean Square Error (RMSE): 211.256190833071
Correlation Coef. (R): 0.907255641665208
R-Square: 0.823112799333349

Parameter          	Best Estimate
--------------------	-------------
mu	0.00762819733960001
beta_1	0.9
beta_2	0.115133226217187
beta_3	1.00000000124386E-5
lambda_1	1.00000000752235E-5
lambda_2	0.218867856024406
lambda_3	0.0170358412000396
gamma_1	1.00000000000006E-5
gamma_2	0.89999999999981
sigma	1.00000000000118E-5
delta	0.899999999999996

if you want to use Matlab to solve this problem, refer to Star Strider's replies like:

1: Fitting and parameter estimation of kinetic model ODE system

2: how to model kinetic reaction system

He has done many excellcent job on this topic.

[Maryam Alyahyai]

Thanks a lot, @Alex Sha. Would you mind please sending me the code so that I know what exactly error I made? Or did you use another solver other than MatLab?

[Alex Sha]

Hi, Maryam Alyahyai, I have told in previous answer that the result had been obtained by using another package named "1stOpt", the codes are as below, very simple. For Matlab solution, you may refer to Star Strider's great replies.

Parameter mu, beta_1, beta_2, beta_3, lambda_1, lambda_2, lambda_3, gamma_1, gamma_2, sigma, delta;
ParameterDomain = [0.00001,0.9];
InitialODEValue t=1, S = 5215663,E = 3779,A = 2103,P = 1512,I = 86,R = 233,N = 5215663+3779+2103+1512+86+233;
Variable t,I;
ODEFunction
S'=mu-beta_1*A*S/N - beta_2*P*S/N  - beta_3*I*S/N-mu*S ;
E'=beta_1*A*S/N  + beta_2*P*S/N  + beta_3*I*S/N - (mu + lambda_1 + lambda_2 + lambda_3)* E ;
A'= lambda_1*E - (mu+gamma_1)*A;
P'= lambda_2*E - (mu+sigma)*P;
I'= lambda_3*E + sigma*P- (mu + gamma_2 + delta)* I;
R'= gamma_1*A + gamma_2*I - mu*R;
N'= -delta*I;
Data;
t=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103;
I=86,144,98,106,102,74,115,93,51,82,143,74,99,36,85,69,98,168,55,154,112,175,174,148,298,322,284,404,157,193,292,372,327,424,463,513,400,348,255,636,811,603,1014,786,576,738,778,770,930,866,604,712,689,1067,1117,1006,1404,1043,745,810,739,852,896,905,1605,1318,1142,1366,1132,919,1197,910,1010,1124,1361,1374,1177,1072,1557,1262,1210,1518,1889,1083,1318,2164,1839,1679,1327,1619,1311,1157,1739,1487,1660,1099,1145,1067,1147,1053,846,665,590;

[Maryam Alyahyai]

Oh, I first thought that "1stOpt" is a tool within MatLab :) Many thanks for your help, replies and patience.