% Code for Figure 3 data generation % First run the models, then you put the figure together % June 2020 by Naven Narayanan % This code takes ~10 seconds to run clear all; close all; clc %%%%%%%%%% Set of all constants in code %%%%%%%%%%% %3 regimes for both cost_m and cost_i - lo(0.2), med(0.5), hi(0.8). Needs to be modified manually cost_m = 0.8; cost_i = 0.5; cost_r = 0.05; sig_r_s_base = 1.0; sig_r_s = 1.0*(1-cost_r); sig_r_i = sig_r_s_base*(1-cost_i)*(1-cost_r); sig_m_s = sig_r_s_base*(1-cost_m); sig_m_i = sig_r_s_base*(1-cost_i)*(1-cost_m); phi_s = 2.0; phi_i = 1.6; t1 = 0.5; t2 = 0.5; S_init = 1000;%arbitrary intial population sizes I_init = 500; %%%%%%% Initializing population array %%%%%%%%%%%% pop = [S_init;I_init]; % Introducing all possible ESS strategies theta_s = [0,1,0,1]; %zeroes imply no movement, ones imply complete movement theta_i = [0,0,1,1]; Eig_set = []; % Array that contains the leading eigenvalue of populations growth for each of the four strategies Strat = []; % 2d array defining which strategy is best for a (gamma, beta) combination. This array is used for plotting the graph d = 1; %counter for values of gamma used while plotting heat map for ga = 1.00:-0.05:0.05 c = 1; %counter for values of gamma used while plotting heat map for be = 0.05:0.05:1.0 %%%% The following terms are parts of alpha, beta, gamma, and delta %%%% which are in turn used to construct trans_mat (matrix %%%% connecting pop at 't' to 't+1' A = sig_r_s*exp(-be*(t1+t2)); B = sig_m_s*exp(-be*t1); C = sig_m_s*(1-(exp(-be*t1)))*(1-(exp(-ga*t2))); D = sig_m_s*(1-(exp(-ga*t2))); E = sig_r_i*(1-(exp(-be*t1))); F = sig_r_i*(1-(exp(-be*t2)))*(exp(-be*t1)); G = sig_m_i*(1-(exp(-be*t1)))*(exp(-ga*t2)); H = sig_r_i; J = sig_m_i*exp(-ga*t2); %a = 1; %defines the starting resident strategy -- can be either (1,2,3, or 4). Needs to be modified manually for a = 1:4 % the loop number of k -- number of times where we introduce mutant strategies to check if invasion occurs. Changing k does not change anything for j = 1:4 % 'j' defines the mutant strategy invading a population with some resident strategy, loop length of j = 3 for i = 1:200 % This loop simulated the population growth for resident strategy -- here 'a' %i S_r = (1-theta_s(a))*A*pop(1,i); I_r = (1-theta_i(a))*H*(pop(2,i)) + (1-theta_i(a))*E*pop(1,i) + (1-theta_s(a))*F*(pop(1,i)); S_m = theta_s(a)*B*pop(1,i) + theta_i(a)*D*(pop(2,i)) + pop(1,i)*C*theta_i(a); I_m = (pop(2,i))*J*theta_i(a) + pop(1,i)*G*theta_i(a); b_max = phi_s*(S_r+S_m)+phi_i*(I_r+I_m); DD = 1*b_max*exp(0.1*(-b_max)); alpha = (A*(1-theta_s(a)) + (B*theta_s(a)) + (C*theta_i(a))); beta = (E*(1-theta_i(a)) + (F*(1-theta_s(a))) + (G*theta_i(a))); gamma = D*theta_i(a); delta = (H*(1-theta_i(a))) + (J*theta_i(a)); trans_mat = [alpha+(DD*((phi_s*alpha)+(phi_i*beta))),gamma+(DD*((phi_s*gamma)+(phi_i*delta)));beta,delta]; %matrix linking pop at time 't' to time 't+1' pop(:,i+1) = trans_mat*pop(:,i); end %pop; Eig_set(a,j) = max(eig(trans_mat)); %Leading eigenvalue of resident strategy %%%% mutant strategy values m_alpha = (A*(1-theta_s(j)) + (B*theta_s(j)) + (C*theta_i(j))); m_beta = (E*(1-theta_i(j)) + (F*(1-theta_s(j))) + (G*theta_i(j))); m_gamma = D*theta_i(j); m_delta = (H*(1-theta_i(j))) + (J*theta_i(j)); trans_mat_ms = [m_alpha+(DD*((phi_s*m_alpha)+(phi_i*m_beta))),m_gamma+(DD*((phi_s*m_gamma)+(phi_i*m_delta)));m_beta,m_delta]; Eig_set(a,j) = max(eig(trans_mat_ms)); %Leading eigenvalues for all mutant strategies end end temp = []; for m = 1:4 if Eig_set(m,1:end~=m) < 1 temp = [temp m]; Strat(d,c) = m; if length(temp) > 1 disp('error1') end elseif Eig_set(m,1:end~=m) <= 1 %m for n = setdiff(1:4,m) if Eig_set(m,n) == Eig_set(m,m) if Eig_set(n,m) == Eig_set(n,n) Strat(d,c) = 10*n + 1*m; elseif Eig_set(n,n) > Eig_set(n,m) Strat(d,c) = n; %else disp('error2') end end end end end c = c+1; %counter used to go to next beta value Eig_set; end d = d+1; %counter used to go to next gamma value %Eig_set end save(strcat(['Fig3_cR=' num2str(cost_r) '_cM=' num2str(cost_m) '_cI=' num2str(cost_i) '.mat']))