Improved SEIR virus model dynamic simulation and theoretical curve

SEIR model

\( \frac{\mathrm{d}}{\mathrm{d}t}E=r\beta I\frac{S}{N}-\alpha E+r_2 \beta_2 E\frac{S}{N} \\ \frac{\mathrm{d}}{\mathrm{d}t}S=-r\beta I\frac{S}{N}-r_2 \beta_2 E\frac{S}{N} \\ \frac{\mathrm{d}}{\mathrm{d}t}R=\gamma I \\ \frac{\mathrm{d}}{\mathrm{d}t}I=\alpha E-\gamma I\)

View all models and descriptions: Infectious disease model simulation and theoretical trends

clear
figure
N = 100;    % total population
E = 0;      % lurker
I = 1;      % infector
S = N-I;    % susceptible
R = 0;      % recovered
 
r = 2;      % number of infected persons contacting susceptible persons
B = 0.005;  % probability of infection
a = 0.1;    % Probability of conversion from latent to infected
r2 = 15;    % of the number of lurkers contacting the susceptible
B2 = 0.015; % probability of latent persons infecting normal people
y = 0.025;  % recovery probability
D = 20;     % Days begin isolation
days = 400;
Dynamic = false;  %Set to true to show the dynamic process

state=zeros(1,N);
index=randperm(N,I);
state(index)=1;
if Dynamic
    axis([0 N 0 N]);
    hold on
end
for i=1:days
    xpos=randperm(N);
    ypos=randperm(N);
    if i>=D
        r=5;
        r2=5;
    end
    for j=1:N
        if state(j)==1
            if rand()<=y
                state(j)=2;
                continue
            end
            dis=sqrt((xpos-xpos(j)).^2+(ypos-ypos(j)).^2);
            l=0;
            for k=1:N
                if state(k)==1 || state(k)==2 || state(k)==3
                    continue
                end
                l=l+1;
                peo(l).dis=dis(k);
                peo(l).num=k;
            end
            T = struct2table(peo);
            sortedT = sortrows(T,'dis');
            sortedS = table2struct(sortedT);
            tp=sortedT{1:r,2};
            for k=1:length(tp)
                if rand()<=B
                    state(tp(k))=3;
                end
            end
        elseif state(j)==3
            if rand()<=a
                state(j)=1;
            end
            dis=sqrt((xpos-xpos(j)).^2+(ypos-ypos(j)).^2);
            l=0;
            for k=1:N
                if state(k)==1 || state(k)==2 || state(k)==3
                    continue
                end
                l=l+1;
                peo(l).dis=dis(k);
                peo(l).num=k;
            end
            T = struct2table(peo);
            sortedT = sortrows(T,'dis');
            sortedS = table2struct(sortedT);
            tp=sortedT{1:r2,2};
            for k=1:length(tp)
                if rand()<=B2
                    state(tp(k))=3;
                end
            end
        end
        
    end
    Idata(i)=length(find(state==1));
    Sdata(i)=length(find(state==0));
    Rdata(i)=length(find(state==2));
    Edata(i)=length(find(state==3));
    %a=-8;
    %b=8;
    %xpos=xpos+ a + (b-a).*rand(1,N);
    %ypos=ypos+ a + (b-a).*rand(1,N);
    if Dynamic && rem(i,10)==0
        phd=scatter(xpos,ypos,[],state,"filled");
        title(['Day ' num2str(i)])
        drawnow
        delete(phd)
    end
end
figure
set(gcf,'visible',true)
days=length(Sdata);
plot(1:days,Sdata,1:days,Idata,1:days,Rdata,1:days,Edata)
hold on
T = 1:days+100;
for idx = 1:length(T)-1
    if S(idx)<0
        S(idx)=0;
    elseif S(idx)>N
        S(idx)=N;
    end
    if E(idx)<0
        E(idx)=0;
    elseif E(idx)>N
        E(idx)=N;
    end
    if I(idx)<0
        I(idx)=0;
    elseif I(idx)>N
        I(idx)=N;
    end
    if R(idx)<0
        R(idx)=0;
    elseif R(idx)>N
        R(idx)=N;
    end
    if idx>=D
        r=5;
        r2=5;
    end
    S(idx+1) = S(idx) - r*B*S(idx)*I(idx)/N(1) - r2*B2*S(idx)*E(idx)/N;
    E(idx+1) = E(idx) + r*B*S(idx)*I(idx)/N(1)-a*E(idx) + r2*B2*S(idx)*E(idx)/N;
    I(idx+1) = I(idx) + a*E(idx) - y*I(idx);
    R(idx+1) = R(idx) + y*I(idx);
    
end

plot(T,S,T,E,T,I,T,R);
plot([D D],[0 N])
legend ('Simulation of susceptible person', 'Simulation of infectious person', 'Simulation of rehabilitation person', 'Simulation of lurking person', 'Suspension person', 'Latency person', 'Infector', 'Rehabilitation person', 'Execution Martial law measures')
xlabel ('day'); ylabel ('number of people')
title ('SEIR Model 2-Martial Law Measures')Code language: Matlab (matlab)

The final screenshot of the simulation is as follows (the dynamic process is not given):

Not implementing martial law measures:

Starting martial law measures from 20 days: