%################################################################################################################################################## % Description: Main raccoon demography and rabies dynamics code to support analyses presented in: % Variation in host home range size decreases rabies vaccination effectiveness in host populations by increasing the spatial spread of rabies virus % Katherine M. McClure, Amy T. Gilbert, Rich Chipman, Erin E. Rees, Kim M. Pepin % % Author: Kim M. Pepin % Contact: Kim.M.Pepin@usda.gov % Code version date: 28 April 2019 % Matlab R2016b Version 9.1.0 %################################################################################################################################################## %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % OUTPUTS: %store(time,16): for each time step % 1:8: 0-1, 1-2, 2-3, 3+ females, 0-1, 1-2, 2-3, 3+ males, % 9:10: number of groups, average size of sounders, % 11:13: number of solo males, number new conceptions, number of new births, % 14:15: number of new solo males, new groups, % 16: number of groups of young males; % 17:18: number and proportion of natural deaths % 19:20: number and proportion of culled individuals % 21:22: number and proportion of newly sterilized individuals % 23:24: group members that move due to overcrowding, solomales that move due to overcrowding % SEIR: for each time step % 1:4: S, E, I, R % 5: new cases % 6: P function mat = main_noSterCull(time,initialk,K,rt,D_age,D_d,male_group_age,dmeth,L_size,L_num,L_age,min_P,gest,long,c_time,... vary,tinfect,within,between,dis_ind_mort,space,net_list,tune,alpha,gamma,rho,... int_vac,prop_vac,space_vac,num_vac,maxlong,Imm0) %I0,int_ster,prop_ster,space_ster,time_ster,min_ster,num_ster,int_cull,prop_cull,space_cull,num_cull, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % INITIALIZE THE POPULATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % INITIAL AGE-STRUCTURE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% coln = {'raccoon id', 'age', 'sex', 'group id', 'D_p', 'D_age', 'D_d', 'age at death', 'age at first conception', 'litter size', ... 'gestation period', 'reproductive status','pregnancy clock', 'post natal clock', 'NHR X', 'NHR Y','DHR X','DHR Y','Dispersal angle',... 'Litters to date','Offspring to date','Sterility clock','S','E','I','R','Habitat grid cell','Ptiter','Pclock','Ppeak','Pdecay','Cell K','Cell class'}; N = sum(initialk); % Initialize population matrix X = []; % Assign age distribution xmax = maxlong; ymax = maxlong; rep = N; rate = [52 long]; temp = zeros(rep,3); % select age & lifespan for i = 1:rep x = xmax+1; y = ymax+1; yy = 0; while x > xmax x = exprnd(rate(1)); end while y > ymax || yy < x y = exprnd(rate(2)); yy = y; end temp(i,:) = [x randsample([0 1],1) y]; % age, sex, lifespan end temp = round(temp); X(:,[2 3 8]) = temp; X(:,1) = 1:size(X,1); %add unique ID's X(:,23) = 1; % Make all susceptible to start X(X(:,2) <= male_group_age & X(:,3) == 1,6) = D_age(2); % enter dispersal age X(X(:,5) == 1 & X(:,2) >= X(:,6),5) = 0; % for the ones that are older than dispersal age already, assume they already dispersed X(X(:,2) <= male_group_age & X(:,3) == 1,5) = 1; %for ones that haven't dispersed yet, set dispersal status to 1 (so they will disperse when it is time) X(X(:,2) > male_group_age & X(:,3) == 1,4) = nan; % assign group id to adult males %%% ASSIGN GROUP MEMBERSHIP %%% % subadult males id = find(X(:,2) > X(:,6) & X(:,2) <= male_group_age & X(:,5) == 1 & X(:,3) == 1); % get young males for grouping temp = randsample(K(2),length(id),'true'); temp2 = []; for i = 1:length(temp) if i == 1 temp2 = [temp2; -1.*ones(temp(i),1)]; else temp2 = [temp2; (min(temp2)-1).*ones(temp(i),1)]; end end X(id,4) = temp2(1:length(id)); % group assignment for subadult males %family groups % reproductive females id = find(X(:,2) > 52 & X(:,3) == 0); % get reproductive-age females temp = randsample(K(2),length(id),'true'); temp2 = []; for i = 1:length(temp) if i == 1 temp2 = [temp2; 1.*ones(temp(i),1)]; else temp2 = [temp2; (max(temp2)+1).*ones(temp(i),1)]; end end X(id,4) = temp2(1:length(id)); % group assignment for reproductive-age females %assign juveniles to groups ave = length(X(:,4) == 0)/length(X(:,4) > 0); % average number that should be in each group id = find(X(:,4) == 0); % get the rest temp2 = []; for i = 1:length(id) grpsize = poissrnd(ave)+1; %choose group size to add if i == 1 temp2 = [temp2; 1.*ones(grpsize,1)]; else temp2 = [temp2; (max(temp2)+1).*ones(grpsize,1)]; end end X(id,4) = temp2(1:length(id)); %%%%%%%%%%%%%%%%%%%%%%%%%% ASSIGN GRID CELLS AND COORDINATES %%%%%%%%%%%%%%%%%%%% % adult males id = find(isnan(X(:,4)) == 1); for i = 1:length(id) hab = randsample(net_list(:,1),1); % pick a habitat id at random HRx = rand*(net_list(hab,6)-net_list(hab,4)) + net_list(hab,4); % get x coordinate for home range centroid of pigs HRy = rand*(net_list(hab,7)-net_list(hab,5)) + net_list(hab,5); X(id(i),[17 18 27 32 33]) = [HRx HRy net_list(hab,[1 8 9])]; end % subadult males id = find(X(:,4) < 0); temp = X(id,4); gid = unique(temp); for i = 1:length(gid) hab = randsample(net_list(:,1),1); % pick a habitat id at random HRx = rand*(net_list(hab,6)-net_list(hab,4)) + net_list(hab,4); % get x coordinate for home range centroid of pigs HRy = rand*(net_list(hab,7)-net_list(hab,5)) + net_list(hab,5); idd = find(X(:,4) == gid(i)); X(idd,[17 18 27 32 33]) = repmat([HRx HRy net_list(hab,[1 8 9])],length(idd),1); end % family groups id = find(X(:,4) > 0); temp = X(id,4); gid = unique(temp); for i = 1:length(gid) hab = randsample(net_list(:,1),1); % pick a habitat id at random HRx = rand*(net_list(hab,6)-net_list(hab,4)) + net_list(hab,4); % get x coordinate for home range centroid of pigs HRy = rand*(net_list(hab,7)-net_list(hab,5)) + net_list(hab,5); idd = find(X(:,4) == gid(i)); X(idd,[17 18 27 32 33]) = repmat([HRx HRy net_list(hab,[1 8 9])],length(idd),1); end den0 = [mean(hist(X(:,27),net_list(:,1))) std(hist(X(:,27),net_list(:,1)))]; % SPATIAL PLOTS TO CHECK THAT INITIAL CONDITIONS ARE WORKING % figure(100); clf; % xid = 0:max(max(net_list(:,[4 6]))); % yid = 0:max(max(net_list(:,[5 7]))); % for i = 1:length(yid) % plot(xid,yid(i).*ones(length(xid)),'color',[0.8 0.8 0.8],'linewidth',0.5); hold on; % end % for i = 1:length(xid) % plot(xid(i).*ones(length(yid)),yid,'color',[0.8 0.8 0.8],'linewidth',0.5);hold on; % end % set(gca,'box','off','xtick',xid,'ytick',yid); [~, ~, xc] = intersect(unique(X(X(:,4)>0,4)),X(:,4)); group0 = [hist(X(X(:,4)>0,4),unique(X(X(:,4)>0,4)))' X(xc,17) X(xc,18)]; [~, ~, xc] = intersect(unique(X(X(:,4)<0,4)),X(:,4)); males0 = [hist(X(X(:,4)<0,4),unique(X(X(:,4)<0,4)))' X(xc,17) X(xc,18)]; ids = find(isnan(X(:,4)) == 1); solo0 = [ones(length(ids),1) X(ids,17) X(ids,18)]; % figure(100); % sc = [5 5]; % scatter(group0(:,2),group0(:,3),'markeredgecolor','r','sizedata',group0(:,1).*sc(1)); hold on; % sounders % scatter(solo0(:,2),solo0(:,3),'markeredgecolor','b','markerfacecolor','b','sizedata',solo0(:,1).*sc(2)); %solo males % scatter(males0(:,2),males0(:,3),'markeredgecolor','g','markerfacecolor','g','sizedata',males0(:,1).*sc(2));%young males % xlabel 'X'; ylabel 'Y'; % group males %%%%%%% REPRODUCTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % age at first conception id = find(X(:,3) == 0); if vary(5) == 0 X(id,9) = L_age(1); else X(id,9) = ones(length(id),1).*normrnd(L_age(1),L_age(2),length(id),1); end X(id,14) = min_P+1; % set initial post-natal clock to above P_min so they can get preganant the first time once they reach their age at first conception %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% START THE MAIN SIMULATION LOOP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %store = zeros(time,25); %0-1, 1-2, 2-3, 3+ females, 0-1, 1-2, 2-3, 3+ males, number of groups, average size of sounders, number of solo males, %new conceptions, new births, number of males going out on their own, new groups %F = []; %Keep track of mother ID, age and litter size of all births %surv = []; POP = zeros(length(time),1); SEIR = zeros(length(time),22); %IMMUNITY = zeros(length(time),4); indexcases = []; %R0 = zeros(1,I0); maxuniqueid = max(X(:,1)); %breach = []; barrier = [];%initialize spread = []; %disp('We started the main loop!') for i = 1:time % For each day if size(X,1) > 0 tempid = find(X(:,24)>0 | X(:,25)>0); temp = unique(X(tempid,27)); % habitats with exposed of infectious individuals if isempty(temp) == 0 SEIR(i,16) = max(net_list(temp,7))-1; %furthest distance from origin end sounders = length(unique(X(find(X(:,4)>0),4))); % current number of family groups X(:,2) = X(:,2)+1; % update ages id = find(X(:,13) > 0); X(id,13) = X(id,13)+1; %update pregnancies id = find(X(:,3) == 0 & X(:,9) < X(:,2) & X(:,12) == 0 & X(:,22) == 0); X(id,12) = 1; % update reproductive status of non-sterilized females id = find(X(:,14) > 0); X(id,14) = X(id,14)+1; %update postnatal clock for non-pregnants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% DISEASE TRANSMISSION if sum(X(:,23)) < size(X,1) %if there are any disease clocks to deal with %[X, newE] = dis(X,space,within,between,alpha,gamma,rho,tinfect,indexcases,tune,dis_ind_mort,i,vary);%dis(X,within,between,alpha,gamma,rho,space,fixed,net_list,i,tinfect,indexcases,tune,dis_ind_mort); [X, newE] = dis2(X,space,within,between,alpha,gamma,rho,dis_ind_mort,vary); SEIR(i,1:15) = [sum(X(:,23)) sum(X(:,24) > 0) sum(X(:,25) > 0) sum(X(:,26) > 0) newE]; % get age-structure of R too % if i > find(tinfect > 0,1,'first') && i < find(tinfect > 0,1,'first')+alpha+gamma+10 % only bother if we are within the infectious period of the first case + some buffer % R0 = R0+R0tally; % end else SEIR(i,1:15) = [sum(X(:,23)) sum(X(:,24) > 0) sum(X(:,25) > 0) sum(X(:,26) > 0) zeros(1,11)]; end % %%%%% CHECK FOR BREACH % id = find(X(:,33) == 4); % get id's for cells beyond barrier % if sum(X(id,24)) + sum(X(id,25)) > 0 % breach = [breach; i]; % record time of breach % %X(id,:) = []; % kill everyone on the other side of the barrier to reduce computation time % end % %%%%% CHECK FOR CASES IN BARRIER ZONE % id = find(X(:,33) == 3); % get id's for cells beyond barrier % if sum(X(id,24)) + sum(X(id,25)) > 0 % barrier = [barrier; i]; % record time to reach barrier % %i = time; % end simulation % end % %%%%% GET IMMUNITY by zone and age class % ageclass = 0:52:52*10; POP(i,1) = size(X,1); %sum(store(:,1:8),2); % for aa = 1:length(ageclass) % if aa == length(ageclass) % id = find(X(:,2) > ageclass(aa)); % else % id = find(X(:,2) > ageclass(aa) & X(:,2) <= ageclass(aa+1)); % end % POP(i,aa+1) = length(id)/size(X,1); % end % SEIR(i,17) = length(find(X(:,33) == 2 & X(:,26) > 0 & X(:,2) <= 52))/length(find(X(:,33) == 2)); SEIR(i,18) = length(find(X(:,33) == 2 & X(:,26) > 0 & X(:,2) > 52 & X(:,2) <= 52*2))/length(find(X(:,33) == 2)); SEIR(i,19) = length(find(X(:,33) == 2 & X(:,26) > 0 & X(:,2) > 52*2 & X(:,2) <= 52*3))/length(find(X(:,33) == 2)); SEIR(i,20) = length(find(X(:,33) == 2 & X(:,26) > 0 & X(:,2) > 52*3 & X(:,2) <= 52*4))/length(find(X(:,33) == 2)); SEIR(i,21) = length(find(X(:,33) == 2 & X(:,26) > 0 & X(:,2) > 52*4 & X(:,2) <= 52*5))/length(find(X(:,33) == 2)); SEIR(i,22) = length(find(X(:,33) == 2 & X(:,26) > 0 & X(:,2) > 52*5))/length(find(X(:,33) == 2)); %%%%% INITIATE DISEASE IN BUFFER ZONE if tinfect(i) > 0 seed = find(X(:,33) == 1); % get X ids in seed zone (input correct seed zone code) if isempty(seed) == 0 temp = round(mean(X(seed,27))); % habitat id in the middle of seeding zone id = find(X(:,27) == temp); % id in X to seed infection if sum(tinfect(1:i)) == 1 % if it is the index case, get R0 indexcases = X(id,1); % remember index cases to get R0 %R0 = zeros(1,length(indexcases)); if Imm0 > 0 % if we want to start with some immune, pick them randomly (not the index case) immunity0 = find(X(:,33) < 3); id2 = randsample(immunity0,round(length(immunity0)*Imm0)); X(id2,26) = poissrnd(rho,length(id2),1)+10; %start with some proportion immune; make sure they have at least 10 days of immunity left; choose time left at random so they don't recover in a pulse X(id2,28) = 1; %make infection history = 1 X(id2,23) = 0; % make sure they are not susceptible SEIR(i,4) = sum(X(:,26) > 0); end end % update info for index case X(id,[23 25 26]) = zeros(length(id),3); X(id,24) = alpha; % set alpha SEIR(i,5) = 1; %record index case as a new infection end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% VACCINATION if int_vac(i) > 0 X = vaccinate(X,prop_vac,space_vac,num_vac,rho); end %%% NATURAL MORTALITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% id = find(X(:,2) > X(:,8)); % get the ones that reach their age of natural death %store(i,17) = length(id); % record the number of deaths from natural causes %store(i,18) = length(id)/size(X,1); % record the proportion of the population dying from natural causes %surv = [surv; X(id,2)]; % store the ages at death to plot survival curve X(id,:) = []; % erase them id = find(X(:,2) > maxlong); X(id,:) = []; % erase ones that are older than maximum lifespan %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% DETERMINE DISPERSAL AND UPDATE DHR %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ADULT MALES ids = find(X(:,2) > male_group_age & X(:,3) == 1 & X(:,4) <= 0); % get males that should go solo now if isempty(ids) == 0 % if there are dispersing adult males offgrid = []; for ww = 1:length(ids) [X, del] = density_dep_disp_v4(X,ids(ww),'solomale',D_d,dmeth,vary,net_list); offgrid = [offgrid del']; X(ids(ww),4) = nan; end %store(i,14) = length(ids); %record the number of males going out on their own X(offgrid,:) = []; % delete the ones going off the grid end % if isempty(ids) % YOUNG MALES id = find(X(:,2) > X(:,6) & X(:,2) <= male_group_age & X(:,5) == 1 & X(:,3) == 1 & X(:,4) > 0); % get the YOUNG males that should disperse now dym = X(id,:); % extract dispersing young males if isempty(id) == 0 % if there are dispersing young males sid = unique(dym(:,4)); %get unique groups of males that will disperse (number of new nodes, etc we need) X(id,5) = 0; % turn off dispersal ability once it happens %store(i,11) = length(id); %add the number of young males that disperse to the number of adult males that disperse offgrid = []; for ww = 1:length(sid) % for each group id_in_all = find(X(:,4) == sid(ww) & ismember(X(:,1),dym(:,1)) == 1); [X, del] = density_dep_disp_v4(X,id_in_all,'male',D_d,dmeth,vary,net_list); %make the ones from the same sounder disperse together offgrid = [offgrid del']; X(id_in_all,4) = min(min(-1,X(:,4)))-1; % disassociate from a sounder and give them a group male id (a negative number) end X(offgrid,:) = []; % delete the ones going off the grid end % if isempty(id) % FEMALES:OVERSIZED FAMILY GROUPS ufem = histc(X(:,4),unique(X(X(:,4)>0,4))); %# of individuals in each family group s_id = unique(X(X(:,4)>0,4)); %ids of unique family group difid = find(ufem-K(1) > 0); g_id = s_id(difid); %get family group id's that have grown > K difK = ufem(difid)-K(1); % get surplus that needs to leave if isempty(g_id) == 0 % if there are some offgrid = []; for k = 1:length(g_id) % for each family group that got too big, disperse a proportion of the group to form a new group % by choosing a random individuals including males idf = find(X(:,4) == g_id(k) & X(:,2) > D_age(1)*30 & X(:,3) == 0); % get all FEMALES from AN oversized group that are > the youngest age for leaving if length(idf) == difK(k) id_move = idf; elseif length(idf) < difK(k) temp_id = find(X(:,4) == g_id(k)); %all ids from family group [~,rem_id] = intersect(temp_id,idf); %ids that are the same as idf temp_id(rem_id) = []; % get rid of the idf ones id_move = [idf; randsample(temp_id,difK(k))]; %choose randomly from the others in the group elseif length(idf) > difK(k) id_move = idf(1:difK(k)); end [X, del] = density_dep_disp_v4(X,id_move,'female',D_d,dmeth,vary,net_list); offgrid = [offgrid del']; X(id_move,4) = max(max(1,X(:,4)))+1; % assign them a new group number (positive value) %store(i,26) = store(i,26)+length(id_move); % store number of females that disperse end X(offgrid,:) = []; % delete the ones going off the grid end % if isempty(g_id) %%%%%%%%% HABITAT-DEPENDENT MORTALITY (FROM OVERCROWDING, choose youngest to die) pops = hist(X(:,27),net_list(:,1)); % get population size in each cell celldif = pops'-net_list(:,8); cellid = find(celldif > 0); % get the cells that are overcrowded if isempty(cellid) == 0 % if the cell is overpopulated, figure out which ones should move vec = []; for k = 1:length(cellid) id_all = find(X(:,27) == cellid(k)); % get row ids of all the pigs in that cell sortage = sortrows([id_all X(id_all,2)],2); % make a matrix with id's and ages and sort die = sortage(1:min(size(sortage,1),abs(celldif(cellid(k))))); % choose ones that will die (picking the youngest ones) vec = [vec die]; end X(unique(vec),:) = []; %store(i,25) = length(unique(vec)); %add up ones that will die this time step end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% NEW CONCEPTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Get all reproductive-age, non-pregnant females, that have been non-pregnant longer than the min time between giving birth and conceiving again id = find(X(:,12) == 1 & X(:,13) == 0 & X(:,14) > min_P & X(:,22) == 0 & X(:,2) <= 78); % 0-1 age-class con = c_time(i)*L_num(1) > rand(length(id),1); %determine who conceives on this day for 0-1 con_id = id(con); %get id's for those that conceive id = find(X(:,12) == 1 & X(:,13) == 0 & X(:,14) > min_P & X(:,22) == 0 & X(:,2) > 78); % 0-1 age-class con = c_time(i)*L_num(2) > rand(length(id),1); %determine who conceives on this day for >1.5 con_id = [con_id; id(con)]; %add id's for those that conceive if isempty(con_id) == 0 % if there are new conceptions... store(i,12) = length(con_id); % store the number of conceptions that happen today for j = 1:length(con_id) % march through and determine litter size based on age class if X(con_id(j),2) <= 78 % get 0-1.5 class if vary(3) == 0 % this will be the litter size once they are born X(con_id(j),10) = L_size(1,1); else X(con_id(j),10) = max(1,round(normrnd(L_size(1,1),L_size(2,1),1,1))); %(con_id(j),10) = max(1,poissrnd(L_size(1,1),1,1)); % make sure there are no -ve's end elseif (X(con_id(j),2) > 78) % >1.5 year class if vary(3) == 0 X(con_id(j),10) = L_size(1,2); else X(con_id(j),10) = max(1,round(normrnd(L_size(1,2),L_size(2,2),1,1))); %X(con_id(j),10) = max(1,poissrnd(L_size(1,2),1,1)); % make sure there are no -ve's end end % set gestation time (does not depend on age-class) X(con_id(j),11) = gest; end % set their preganancy clocks to 1 X(con_id,13) = 1; X(con_id,14) = 0; % turnoff their post-natal clocks end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% NEW BIRTHS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% id = find(X(:,13) > X(:,11)); % get the pregnant females that reached their target gestation time if isempty(id) == 0 % if there are some births, get each litter and initialize each individual in it %newb = 0; for j = 1:length(id) % for each female's litter Y = zeros(X(id(j),10),size(X,2)); % make a matrix for the new individuals from mother id(j) based on the litter size Y(:,23) = 1; %make everyone susceptible Y(:,[27 32 33]) = repmat(X(id(j),[27 32 33]),size(Y,1),1); % assign the mother's habitat info ff = round(X(id(j),10)*rt(1)); mm = X(id(j),10)-ff; % determine the number of each sex in the litter using the rt values Y(ff+1:ff+mm,3) = 1; % label the males % determine dispersal info for the males Y(ff+1:ff+mm,5) = 1; %binornd(1,D_p(2),mm,1); % enter 1 for males that disperse according to prob of dispersal (set to 1 if all males will disperse) if vary(1) == 0 Y(ff+1:ff+mm,6) = D_age(2); % enter dispersal age else Y(ff+1:ff+mm,6) = repmat(normrnd(D_age(2),D_age(4),1),mm,1); %gives a random age for dispersal - same for all males in the litter end % determine age of first births for females if vary(5) == 0 Y(1:ff,9) = L_age(1); else Y(1:ff,9) = ones(ff,1).*normrnd(L_age(1),L_age(2),ff,1); end Y(1:ff,14) = min_P+1; % set initial post-natal clock to above P_min so they can get preganant the first time once they reach their age at first conception Y(:,2) = 1; % add everybody's age Y(:,1) = maxuniqueid+(1:size(Y,1))'; % add everybody's unique id maxuniqueid = max(Y(:,1)); %update max unique identifier before culling so we don't reassign any numbers Y(:,[4 15 16 17 18]) = repmat(X(id(j),[4 17 18 17 18]),size(Y,1),1); % put everybody in their mother's group and current HR Y(:,8) = round(exprnd(long,size(Y,1),1)); % give them each a death age X = [X; Y]; % add the newborns to the main population %newb = newb+size(Y,1); % store the new births % Keep track of some reproductive stats %X(id,20) = X(id,20)+1; %add 1 to litters %X(id,21) = X(id,21)+size(Y,1); %add offspring to litters end %store(i,13) = newb; % store the new births % Reset the pregnancy clocks, gestation time and litter size for females that gave birth X(id,[10 11 13]) = 0; X(id,14) = 1; % start their post-natal clock end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %convert outputs to monthly scale % dist (1), all cases (1), cases in each zone (2), new cases by sex/age (0-1, 1, 2, 3+) (8), immunity % by age class (6), pop (1) % divide year into quarters: 13 week periods scale = 13; mat = zeros(length(0:scale:time),18); mat(:,1) = Retime(SEIR(:,16),scale,'point'); %distance mat(:,2) = Retime(SEIR(:,5),scale,'sum'); %all cases mat(:,3) = Retime(SEIR(:,14),scale,'sum'); %all cases zone 2 mat(:,4) = Retime(SEIR(:,15),scale,'sum'); % all cases zone 3 mat(:,5) = Retime(SEIR(:,6),scale,'sum'); % cases: F 0-1 mat(:,6) = Retime(SEIR(:,7),scale,'sum'); % cases: M 0-1 mat(:,7) = Retime(SEIR(:,8),scale,'sum'); % cases: F 1 mat(:,8) = Retime(SEIR(:,9),scale,'sum'); % cases: M 1 mat(:,9) = Retime(SEIR(:,10),scale,'sum'); % cases: F 2 mat(:,10) = Retime(SEIR(:,11),scale,'sum'); % cases: M 2 mat(:,11) = Retime(SEIR(:,12),scale,'sum'); % cases: F 3 mat(:,12) = Retime(SEIR(:,13),scale,'sum'); % cases: M 3 mat(:,13) = Retime(SEIR(:,17),scale,'med'); %zone 2 immunity 0-1 mat(:,14) = Retime(SEIR(:,18),scale,'med'); %zone 2 immunity 1 mat(:,15) = Retime(SEIR(:,19),scale,'med'); %zone 2 immunity 2 mat(:,16) = Retime(SEIR(:,20),scale,'med'); %zone 2 immunity 3 mat(:,17) = Retime(SEIR(:,21),scale,'med'); %zone 2 immunity 4 mat(:,18) = Retime(SEIR(:,22),scale,'med'); %zone 2 immunity 5 mat(:,19) = Retime(POP,scale,'point'); %abundance %%% RECORD OUTPUTS % Probability of a rabies case breaching the ORV barrier % Time until breach % Case reduction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % left = find(X(:,4) == 0); % end