% Experiemental data on 12 months
time = 1:12 ; % experimental time 12 months
Yc = [20 18 19 24 28 26 25 24 20 19 21 19]' ; % Bu cases
W = [0.14 0.08 0.18 0.08 0 0.37 0.24 0.18 0 0.1 0.08 0.05]'; % MU prevalence in vectors
U = [0.05 0.02 0.05 0.07 0 0.09 0.1 0.16 0.05 0.13 0.07 0.06]'; %MU positivity
beta1 = 0.1; beta2 = 0.1; sigma = 0.2; eps = 0.16; gamma = 0.17;
eta = 0.7; alpha = 0.8; Delta = 0.7; theta = 0.9;
y0 = [6500; 0; 20; 0; 0; 0.63; 0.37; 1] ; a=1; b=1;
mu = 0.00216; % Mortality/fertility rate 2,6% by year
Lambda = @(t) ??? % Your time-dependent function for Lambda
p0 = [beta1; beta2; sigma; eps; gamma; eta; alpha; Delta; theta; y0; a; b] ; % Initial guess for parameters
%Parameter upper, lower, and Using lsqcurvefit to adjust model parameters
lb = [0, 0, 1/5, 1/6, 1/6, 0, 0, 0, 0, 0,0,0,0,0,0,0,0, -inf,-inf]; % lower bound
ub = [1, 1, 1/3, 1.0, 1/2, 1, 1, 1, 1, 10^5,50,50,50,50,2,2,3, inf,inf]; % upper bound
opts = optimoptions('lsqcurvefit');
opts.MaxFunEvals = 10000;
[p_fit,resnorm,residual,exitflag,output,jacobian] = lsqcurvefit(@(p,t) fun_model(p,t,mu,Lambda),p0,time,[Yc,W,U],lb,ub,opts);
function Z = fun_model(p,t,mu,Lambda)
% Unknown parameters
beta1 = p(1);
beta2 = p(2);
sigma = p(3);
eps = p(4);
gamma = p(5);
eta = p(6);
alpha = p(7);
Delta = p(8);
theta = p(9);
y0 = p(10:17); % initial conditions
a = p(18);
b= p(19);
% systems
f = @(t,y) [ mu*(y(1)+y(2)+y(3)+y(4)+y(5)) - mu*y(1) - beta1*y(1)*y(7) - beta2*y(1)*y(8); ...
beta1*y(1)*y(7) + beta2*y(1)*y(8) - sigma*y(2) - mu*y(2); ...
sigma*y(2) - eps*y(3) - mu*y(3); ...
eps*y(3) - gamma*y(4) - mu*y(4); ...
gamma*y(4) - mu*y(5); ...
eta*(y(6)+y(7)) - alpha*y(6)*y(8) - eta*y(6); ...
alpha*y(6)*y(8) - eta*y(7); ...
(a+b*Lambda(t))*y(8) - Delta*y(8) - theta*eta*y(7) ];
% Resolution with ode45
options = odeset('RelTol',1e-6,'AbsTol',1e-6); % Adjust tolerances for performance
[~,ypred] = ode45(f,t,y0,options);
% Return the desired output
Z = [ypred(:,3), ypred(:,7), ypred(:,8)];
end
