Problem with simulating with SEIR variant

I am currently trying to get results for simulating a SEIR variant model but I am not getting any curves for 3 out of 7 of the variables (I, Q and D). I think the problem came in when I implemented xA and xI into the equations. What i was hoping for is that I wanted to let xA and xI to remain a constant (which will be subtracted from the (3) & (4) respectively. But at the same time, should Y(3) and Y(4) be smaller than the value of xA and xI, the constant will change accordingly to be equal to Y(3) and Y(4), as Y(3) and Y(4) should not be a negative number. Lowest value to remain at 0.

Been trying to troubleshoot it but to no avail. Not sure if anyone can provide me with some insights on this too? Thanks alot :)

%SEIAQRD_Main

to = 1;
tf = 365;
tspan = [to:1:tf];
N = 5000000;
y0 = [N-100 0 50 50 0 0 0]; %Initial conditions
[t,Y] = ode45(@SEIAQRD, tspan, y0);
plot(t,Y)



%SEIAQRD_function_file
function dYdt = SEIAQRD(t,Y)

tests_conducted = 2900;

%Input R0 value for simulation
R0 = 2; %Range of 2.0 to 3.5

%Calculating beta(B)
Duration = 11;
B = R0/Duration;
Trans = 0.13;
c = B/Trans;

%Parameters
N = 5703600;
ep = 1/5.2;
gamma = 1/Duration; 
gammaQ = 1/21; 
alpha = 0.5; 
Miu = (30/50000)*100; 
avg_positive = ((74+49+65+75+120+66+106)/7)/2900;
x = tests_conducted*avg_positive; 
xI = (2/3)*x; 
xA = (1/3)*x; 
tau = 1/1; 

%Define equations
%Y(1)=S, Y(2)=E, Y(3)=I, Y(4)=A
%Y(5)=Q, Y(6)=R, Y(7)=D, 
dYdt = zeros(7,1);
lambda = (c*(Y(3)/N)) + (c*(Y(4)/N));

%Susceptible (S)
dYdt(1) = -(lambda*Y(1));
%Exposed (E)
dYdt(2) = (lambda*Y(1))-(ep*Y(2));
%Symptomatic (I)
if xI>Y(3)
    xI1=Y(3);
elseif xI<=Y(3)
    xI1=xI;
end
dYdt(3) = (ep*alpha*Y(2))-(gamma*Y(3))-(xI1*tau);
%Asymptomatic (A)
if xA>Y(4)
    xA1=Y(4);
elseif xA<=Y(4)
        xA1=xA;
end
dYdt(4) = (ep*(1-alpha)*Y(2))-(gamma*Y(4))-(xA1*tau);
%Quarantine (Q)
dYdt(5) = ((xI1+xA1)*tau)-(gammaQ*Y(5))-(Miu*Y(5));
%Recovered (R)
dYdt(6) = (gammaQ*Y(5))+(gamma*Y(4));
%Death (D)
dYdt(7) = Miu*Y(5);
end

Answer 1

[Alan Stevens]

I get your results for I, Q and D. See them clearly by changing your Main section to;

%SEIAQRD_Main
to = 1;
tf = 365;
tspan = to:tf;
N = 5000000;
y0 = [N-100 0 50 50 0 0 0]; %Initial conditions
[t,Y] = ode45(@SEIAQRD, tspan, y0);
plot(t,Y)
legend('S','E','I','A','Q','R','D')
figure(2)
plot(t,Y(:,3))
legend('I')
figure(3)
plot(t,Y(:,5),t,Y(:,7))
legend('Q','D')

[Clement Koh]

I see, so i guess it's because that it is too small to be able to see on the normal scale?

[Alan Stevens]

Yes.

[Clement Koh]

Thank you!