%9/16/2014
%EbolaNetwork models SEIR
%with one transmission rate within families
%another transmission rate within an extended network of N families
%the extended network occurs within a radius:
%interaction between R nearest neighbors on either side
%the larger population effectively infinite
%...error message if the larger population becomes limiting
%three scales
%household (H members)
%community (comprised of C=2*r+1 families)
%population (effectively infinite)
%rates of infection
%RzeroH: average number of people in household infected by 1 sick person
%RzeroC: average number of people in community infected by 1 sick person
%LATTICE STATES
%0=healthy
%1=exposed
%2=infectious
%3=refractory
function [Re,NumFamiliesExposed,symptomatic_cases,success]=EbolaNetwork0918(H,C,RzeroH,RzeroC,MaxCases);
%Fixed model parameters
avg_days_incubation=5.3;%or 7 %this is 1/sigma in the model
avg_days_infectious=5.61;%or 10 %this is 1/gamma in the model
global Community
global Parentage
global symptomatic_cases
global new_cases_per_day
global NumFamiliesExposed
RofHinf=RzeroH/(avg_days_infectious*(H-1)); %daily rate per person of infection within family
RofCinf=RzeroC/(avg_days_infectious*(C*H-1));%daily rate of infecting someone in near family
community_radius=((C-1)/2);
if fix(community_radius)~=community_radius
'error-community radius is not an integer, C must be an odd number'
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Initialize Lattices And Data Arrays %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TotalNumFamilies=50*max(2*community_radius,1);%this is the parameter 'L' in the manuscript
%['total population size is ', num2str(H*TotalNumFamilies),' people']
%an error message may mean you need to increase this
Community=zeros(H,TotalNumFamilies);
Parentage=zeros(size(Community));
new_cases_per_day=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Initial Conditions: One exposed person %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Community(1,fix(TotalNumFamilies/2))=1;
symptomatic_cases=sum(sum(Community>=2));
day=1;
%this is what happens each day
keep_going=1;
success=0;
while keep_going==1
%infections occur within a household
for i=1:TotalNumFamilies
%infections occur within families
ThisFamily=Community(:,i); %single out the ith household
%list of infected persons
infectious=(ThisFamily==2);
num_infectious=sum(sum(infectious));
if num_infectious>0
[Ix,Iy]=find(infectious>0);
susceptible=(ThisFamily==0);
for j=1:num_infectious
newly_exposed=susceptible.*(rand(H,1)0);
num_infectious=sum(sum(infectious));
for j=1:num_infectious
%community_radius=((C-1)/2);
RadiusSusceptible=zeros(size(Community));
RadiusSusceptible(:,(Iy(j)-community_radius):(Iy(j)+community_radius))=1;
susceptible=(Community==0).*RadiusSusceptible;
newly_exposed=susceptible.*(rand(size(Community))=2)));
%infected individuals may become refractory
rateItoR=1/avg_days_infectious;
Community=Community+(Community==2).*(rand(size(Community))1))>MaxCases
keep_going=0;
%[num2str(MaxCases),' Cases Reached']
success=1;
end
end
NumFamiliesExposed=sum(sum(Community)>0);
%calculating R0
Re=mean(Parentage(Community==3));