传染病模型SIS及相应的matlab代码_sis传染病模型matlab代码

传染病模型：
常见的传染病模型按照具体的传染病的特点可分为 SI、SIS、SIR、SIRS、SEIR 模型。其中“S”“E”“I”“R”的现实含义如下：

S (Susceptible)，易感者，指缺乏免疫能力健康人，与感染者接触后容易受到感染；

E (Exposed)，暴露者 ，指接触过感染者但暂无传染性的人，可用于存在潜伏期的传染病；

I (Infectious)，患病者，指有传染性的病人，可以传播给 S，将其变为 E 或 I ；

R (Recovered)，康复者，指病愈后具有免疫力的人，如是终身免疫性传染病，则不可被重新变为 S 、E 或 I ，如果免疫期有限，就可以重新变为 S 类，进而被感染。

SIS传染病模型假设：
1.在疾病传播期内，研究区域内总人数N不变,既不考虑生死，也不考虑迁移。将该区域内的人群分为：易感染者、感病者和病愈免疫的移出者三类，以下简称易感染者，病人和移出者。在t时刻将这三类人群的数量分别记为S(t)、I(t)和R(t)，t时刻这三类人在总人数中所占的比例分别记为s(t)、i(t)和r(t)。
2.每个病人每天有效接触的平均人数是常数λ，称为日接触率。当病人与健康者有效接触时,会使健康者受感染变为病人。
3.每天被治愈的病人占病人总数的比例为常数μ，称为日治愈率。病人被治愈后仍有可能被感染为病人，故假设1/μ是该疾病的平均传染期。
4.在初始时刻，只有少数个体处于感染状态，其他都是易感染状态。
5.假设疾病的时间尺度远小于个体生命周期，从而不考虑个体的出生和自然死亡对数据统计造成的影响。
6.完全混合：每一个个体与其他个体接触的机会均等。
根据以上假设可得：每个病人每天可使λ·s(t)个健康者变为病人，每天有μ·Ni个病人被治愈。
σ=λ/μ
由该式中，接触数σ是病人在平均传染期内有效接触的人数。

SIS传染病模型构建：
由假设1可得，整个群体由易感染者，感染者，移出者构成，三者之间的关系如下：
根据假设1可得：

由上式可得，关于S(t),I(t),R(t)的微分方程：

进而得到关于s(t),i(t),r(t)的微分方程：

将数据带入上述公式，求解微分方程即可得到各个变量随时间变化的情况。

matlab代码：

function [t,S,I] = Program_7_7(N,n,tau,gamma,MaxTime,Type)

% It is an SIS disease spread through a network. Allowed
% network types are 'Random','Spatial','Lattice' and 'SmallWorld'
%
% We assume N individuals, each with an averge of n contacts.

% In this model we define an individual by their status flag:
% Status=1   =>  Susceptible
% Status=2   =>  Infectious
% Status=0   =>  Recovered


% Sets up default parameters if necessary.
if nargin == 0
    N=100;
    n=4;
    tau=1;
    gamma=0.1; 
    MaxTime=1000;
    Type='Random';
end

% Checks all the parameters are valid
CheckGreater(N,0,'Number of individuals N');
CheckGreater(n,0,'Number of neighbours n');
CheckGreater(tau,0,'tau');
CheckGreater(gamma,0,'gamma');
CheckGreater(MaxTime,0,'MaxTime');

%Initialise the Network
% (X,Y) is location, G is the network graph matrix
% this means we use S and I for the number of susceptibles and infecteds

[X,Y,G,N]=Create_Network(N,n,Type);
Status=1+0*X; Status(1)=2;
Rate=0*X; Rate(1)=gamma; Rate(find(G(:,1)))=tau;

t=0; i=1; S=N-1; I=1;
% The main iteration
subplot(2,1,1);
[j,k,s]=find(G);
plot([X(j) X(k)]',[Y(j) Y(k)]','-k');
hold on
Col=[0.7 0.7 0.7; 0 1 0; 1 0 0];
for k=1:N
    H(k)=plot(X(k),Y(k),'ok','MarkerFaceColor','g');
end
set(H(1),'MarkerFaceColor','r');
hold off;
drawnow;

while (t<MaxTime & I(end)>0)
    [step,Rate,Status,e]=Iterate(Rate,Status,G,tau,gamma);
    i=i+1;
    t(i)=t(i-1)+step;
    S(i)=length(find(Status==1));
    I(i)=length(find(Status==2));
    set(H(e),'MarkerFaceColor',Col(Status(e)+1,:));

    %     subplot(2,1,1);
    %     [j,k,s]=find(G);
    %     plot([X(j) X(k)]',[Y(j) Y(k)]','-k');
    %     hold on
    %     s=find(Status==0); plot(X(s),Y(s),'ok','MarkerFaceColor',[0.7 0.7 0.7],'MarkerSize',8);
    %     s=find(Status==1); plot(X(s),Y(s),'ok','MarkerFaceColor','g','MarkerSize',8);
    %     s=find(Status==2); plot(X(s),Y(s),'ok','MarkerFaceColor','r','MarkerSize',10);
    %     hold off;

    subplot(4,1,3);
    plot(t,S,'-g');
    ylabel 'Number of Susceptibles'
    subplot(4,1,4);
    plot(t,I,'-r');
    ylabel 'Number of Infecteds'
    xlabel 'Time'
    drawnow;
    [];
end

% Create the Network
function [X,Y,G,N]=Create_Network(N,n,Type);

if n > (N-1)
    error('Impossible to have an average of %d contacts with a population size of %d',n,N);
end

G=sparse(1,1,0,N,N);
X=rand(N,1); Y=rand(N,1);
switch Type

    case {'Random','random'}
        contacts=0;
        while(contacts<n*N)
            i=ceil(rand(1,1)*N); j=ceil(rand(1,1)*N);
            if i~=j & G(i,j)==0
                contacts=contacts+2;
                G(i,j)=1; G(j,1)=1;
            end
        end
                
    case {'Spatial','spatial'}
        D=(X*ones(1,N) - ones(N,1)*X').^2 + (Y*ones(1,N) - ones(N,1)*Y').^2;
        Prob=exp(-D*5)-rand(N,N); Prob=triu(Prob,1)-1e100*tril(ones(N,N),0);
        [y i]=sort(reshape(Prob,N*N,1));
        p=i(end+[(1-n*N/2):1:0]);
        i=1+mod(p-1,N); j=1+floor((p-1)/N);
        G=sparse([i; j],[j; i],1,N,N);

    case {'Lattice','lattice'}
        if mod(sqrt(N),1)~=0
            warning('N=%d is not a square number, rounding to %d',N,round(sqrt(N)).^2);
            N=round(sqrt(N)).^2;
        end
        [X,Y]=meshgrid([0:(sqrt(N)-1)]/(sqrt(N)-1));
        X=reshape(X,N,1); Y=reshape(Y,N,1);
        D=(X*ones(1,N) - ones(N,1)*X').^2 + (Y*ones(1,N) - ones(N,1)*Y').^2;
        Prob=triu(D,1)+1e100*tril(ones(N,N),0);
        [y i]=sort(reshape(Prob,N*N,1));
        p=i(1:(n*N/2 - 2*sqrt(N)));
        i=1+mod(p-1,N); j=1+floor((p-1)/N);
        G=sparse([i; j],[j; i],1,N,N);

    case {'SmallWorld','Smallworld','smallworld'}
        if mod(sqrt(N),1)~=0
            warning('N=%d is not a square number, rounding to %d',N,round(sqrt(N)).^2);
            N=round(sqrt(N)).^2;
        end
        [X,Y]=meshgrid([0:(sqrt(N)-1)]/(sqrt(N)-1));
        X=reshape(X,N,1); Y=reshape(Y,N,1);
        D=(X*ones(1,N) - ones(N,1)*X').^2 + (Y*ones(1,N) - ones(N,1)*Y').^2;
        Prob=triu(D,1)+1e100*tril(ones(N,N),0);
        [y i]=sort(reshape(Prob,N*N,1));
        p=i(1:(n*N/2 - 2*sqrt(N)));
        i=1+mod(p-1,N); j=1+floor((p-1)/N);
        G=sparse([i; j],[j; i],1,N,N);
        % Already made a lattice, now add R random connections
        R=10;
        for k=1:R
            i=1; j=1;
            while (i==j | G(i,j)==1)
                i=ceil(rand(1,1)*N); j=ceil(rand(1,1)*N);
            end
            G(i,j)=1; G(j,1)=1;
        end
        
    otherwise
        error('Do not recognise network type %s',Type);
end



%Do the Up-Dating.
function [step,Rate,Status,Event]=Iterate(Rate,Status,G,tau,gamma);

Sum=sum(Rate); Cum=cumsum(Rate);

Event=min(find(Cum>rand(1,1)*Sum));

Status(Event)=mod(Status(Event)+1,3);

contacts=find(G(:,Event) & Status==1);
switch Status(Event)
    case 1
        
    case 2
        Rate(Event)=gamma;
        Rate(contacts)=Rate(contacts)+tau;
        G(Event,:)=0;
    case 0
        % For SIR type dynamics we require the following 2 lines
        Rate(Event)=0;
        Rate(contacts)=Rate(contacts)-tau;

        % For SIS type dynamics we require the following 3 lines
    %Status(Event)=1;
        %Rate(Event)=tau*length(find(G(:,Event) & Status==2));
        %Rate(contacts)=Rate(contacts)-tau;

end
Rate=Rate.*sign(Status);

step=-log(rand(1,1))/Sum;

% Does a simple check on the value
function []=CheckGreaterOrEqual(Parameter, Value, str)

m=find(Parameter<Value);
if length(m)>0
    error('Parameter %s(%g) (=%g) is less than %g',str,m(1),Parameter(m(1)),Value);
end

function []=CheckGreater(Parameter, Value, str)

m=find(Parameter<=Value);
if length(m)>0
    error('Parameter %s(%g) (=%g) is less than %g',str,m(1),Parameter(m(1)),Value);
end
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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198

结果