function [I_res,N_res,m_res_realized,eigJ,IavgNM,IavgM] = do_run_simulation(eflag,lambda,beta_min,gamma,alpha_s,mu,psi,phi_max,alpha_f,I_res,N_res,m_res_strategy,m_mut_strategy)
% written by Allison Shaw, Contact ashaw@umn.edu for assistance.
% started: 14 September 2021
% last updated: 18 October 2021
%
% This simulates a population under one of two scenarios:
%   'eco': ecological dynamics only, simulates a single population over
%       many years for a given set of parameters and outputs the final
%       population details. 
%   'evo': evolutionary dynamics, simulates a population with both an
%       abundant resident type and a rare mutant type over Tevo (=1) year for
%       a given set of parameters and outputs the final population details.
%
% INPUTS:
%   eflag = whether to run 'eco' or 'evo' simulations
%   lambda = time to half change in infection
%   beta_min = scalar of infection dynamics at x=1
%   gamma = midpoint of survival-v-intensity
%   alpha_s = slope of survival-v-intensity
%   mu = slope of survival-v-migration
%   I_res = starting intensity of each intensity stage (column vector) (use
%     for residents and migrants in evolutionary simulations)
%   psi = degree of migratory culling (interaction infection/migration)
%   phi_max = maximum per capita fecundity
%   alpha_f = slope of fecundity-v-intensity sims)
%   N_res = for residents, starting number of individuals in each intensity stage (column vector)
%   m_res_strategy = for residents, migration strategy (column vector) for all possible intensity stages
%   m_mut_strategy = for mutant, migration strategy (column vector) for all possible intensity stages (if applicable)
%
% OUTPUTS:
%   I_res = final intensity of each intensity stage (column vector) (use
%     for residents and migrants in evolutionary simulations)
%   N_res = for residents, final number of individuals in each intensity stage (column vector)
%   m_res_realized = for residents, vector of the migration distances that each intensity stage
%      present in the population goes
%   eigJ = for mutants, growth rate from matrix
%   IavgNM = average infection intensity of Non-Migrants
%   IavgM = average infection intensity of Migrants

%%---MODEL SET UP---------------------------------------------------------%
    tau_1 = 0.2;     % time of outward migration
    tau_2 = 1-tau_1; % time of return migration
    delta = 0.001;   % strength of density-dependence
    beta_max = (tau_2 + lambda)./(tau_1+tau_2+lambda); % scalar of infection dynamics at x=0

    Nvia = 0.1;  % viable stage size - remove stages where size falls below this
    Teco = 200;  % number of years to run for
    Teco_max = 1000;  % absolute max number of years to run for
    Tevo = 1; % number of years to run for (not set up to run eco for more than 1 year)
     % as is I_mut doesn't change over time b/c assuming only one year
    eps1 = 0.1;   % error for population size (must be within this to be same size)
%%---MODEL SET UP---------------------------------------------------------%

% check parameter values
if alpha_f>phi_max;   error('negative fecundity: decrease alpha_f or increase phi_max'); end
if beta_min>beta_max; error('migration increases infection: decrease beta_min'); end

if strcmp(eflag,'eco') % if ecological simulations
    T = Teco;
    I_mut = []; % ignore mutants
elseif strcmp(eflag,'evo') % if evolutionary simulations
    T = Tevo;
    I_mut = I_res; % for mutant, starting intensity of each intensity stage (column vector)
else
    error('unrecognizable eflag')
end
    
[Nsum_res,Iavg_res] = deal(NaN(T,1));
[IavgNM,IavgM] = deal(NaN(4,1));

t = 0;

% run simulations for T time or until population crashes
while t < T && isempty(N_res)==0
    t=t+1;
    
    % determine the migration strategy, based on I
    m_res_realized = do_get_migration_distance(I_res,m_res_strategy);
    if strcmp(eflag,'eco')
        m_mut_realized = []; % ignore mutants
    elseif strcmp(eflag,'evo')
        m_mut_realized = do_get_migration_distance(I_mut,m_mut_strategy); % get mutant strategy
    end

    % find index corresponding to Non-Migrants than Migrants
    indNM = find(m_res_realized==0);
    indM = find(m_res_realized>0);

    % calculate average infection intensity of Non-Migrants and Migrants
    IavgNM(1) = sum(I_res(indNM).*N_res(indNM))./sum(N_res(indNM));
    IavgM(1)  = sum(I_res(indM).*N_res(indM))./sum(N_res(indM));
    
    % determine survival, based on I and m
    %   (don't let survival go negative)
    s1_res = max(0,1 - 1./(1+exp(-alpha_s.*(I_res-gamma))) - mu.*m_res_realized - psi.*I_res.*m_res_realized);
    s1_mut = max(0,1 - 1./(1+exp(-alpha_s.*(I_mut-gamma))) - mu.*m_mut_realized - psi.*I_mut.*m_mut_realized);
    
    % calculate population size
    N_tau1_res = s1_res.*N_res;
    
    % calculate average infection intensity of Non-Migrants and Migrants
    IavgNM(2) = sum(I_res(indNM).*N_tau1_res(indNM))./sum(N_tau1_res(indNM));
    IavgM(2)  = sum(I_res(indM).*N_tau1_res(indM))./sum(N_tau1_res(indM));

    %---TIME TAU 2 PASSES HERE ---%

    % get intensity scaling parameter
    beta_tmp_res = (beta_min-beta_max)*m_res_realized + beta_max;
    beta_tmp_mut = (beta_min-beta_max)*m_mut_realized + beta_max;
    
    % calculate infection intensity
    I_tau2_res = NaN(size(I_res));
    for j = 1:length(beta_tmp_res)
        beta = beta_tmp_res(j);
        if beta>=0 % for positive beta
            I_tau2_res(j) = I_res(j) + beta*(1-I_res(j)).*tau_2./(tau_2+lambda);
        else % for negative beta
            I_tau2_res(j) = I_res(j) + beta*I_res(j).*tau_2./(tau_2+lambda);
        end
    end
    I_tau2_mut = NaN(size(I_mut));
    for j = 1:length(beta_tmp_mut)
        beta = beta_tmp_mut(j);
        if beta>=0 % for positive beta
            I_tau2_mut(j) = I_mut(j) + beta*(1-I_mut(j)).*tau_2./(tau_2+lambda);
        else % for negative beta
            I_tau2_mut(j) = I_mut(j) + beta*I_mut(j).*tau_2./(tau_2+lambda);
        end
    end
    if sum(isnan(I_tau2_res)) > 0; error('NaN infection'); end
    
    % calculate average infection intensity of Non-Migrants and Migrants
    IavgNM(3) = sum(I_tau2_res(indNM).*N_tau1_res(indNM))./sum(N_tau1_res(indNM));
    IavgM(3)  = sum(I_tau2_res(indM).*N_tau1_res(indM))./sum(N_tau1_res(indM));

    % determine survival, based on I_tau2 and m
    %   (don't let survival go negative)
    s2_res = max(0,1 - 1./(1+exp(-alpha_s.*(I_tau2_res-gamma))) - mu.*m_res_realized - psi.*I_tau2_res.*m_res_realized);
    s2_mut = max(0,1 - 1./(1+exp(-alpha_s.*(I_tau2_mut-gamma))) - mu.*m_mut_realized - psi.*I_tau2_mut.*m_mut_realized);
   
    % calculate population size
    N_tau2_res = s2_res.*N_tau1_res;

    % calculate average infection intensity of Non-Migrants and Migrants
    IavgNM(4) = sum(I_tau2_res(indNM).*N_tau2_res(indNM))./sum(N_tau2_res(indNM));
    IavgM(4)  = sum(I_tau2_res(indM).*N_tau2_res(indM))./sum(N_tau2_res(indM));

    % reproduction
    f_res = phi_max - alpha_f*I_tau2_res; % density-independent fecundity
    f_mut = phi_max - alpha_f*I_tau2_mut; % density-independent fecundity
    
    d = exp(-delta*(sum(N_tau2_res)));  % strength of density-dependence
    newborns_res = sum(f_res.*N_tau2_res*d);  % total number of newborns

    % save overall populations, adding newborns at the start with no infection
    I_new = [0; I_tau2_res];
    N_new = [newborns_res; N_tau2_res];
    
    % find intensity stages that have too few individuals and remove
    ind = find(N_new<Nvia);
    N_new(ind) = [];
    I_new(ind) = [];
    clear ind
    
    %---TIME TAU 1 PASSES HERE ---%
    
    % calculate infection intensity
    I_tau1_res = I_new + (1-I_new).*tau_1./(tau_1+lambda);
    I_tau1_mut = I_new + (1-I_new).*tau_1./(tau_1+lambda);
    if sum(isnan(I_tau1_res)) > 0; error('NaN infection'); end

    % census population and record metrics
    I_res = I_tau1_res;
    N_res = N_new;
    Nsum_res(t) = sum(N_res);
    if ~isempty(I_res)
        Iavg_res(t) = sum(I_res.*N_res)./sum(N_res);
    end
    
    % check for convergence and extend time if possible
    if t==Teco
        if t>=Teco_max
            disp('max simulation time reached')
        elseif abs(Nsum_res(end-1)-Nsum_res(end))>eps1
            Teco = Teco + 100;
            T = Teco;
        end
    end

end

if strcmp(eflag,'eco')
    eigJ = []; % ignore mutants and growth rate
    
    % check for convergence
    if abs(Nsum_res(end-1)-Nsum_res(end))>eps1; error('population size not converged'); end
    if abs(Iavg_res(end-1)-Iavg_res(end))>eps1; error('average intensity not converged'); end

elseif strcmp(eflag,'evo')
    % construct Jacobian matrix for ESS analysis
    Svec = s2_mut.*s1_mut;
    Fvec = f_mut*d';
    J = diag(Svec,-1);
    J(1,1:end-1) = Fvec;
    eigJ = max(eig(J));
end

m_res_realized = do_get_migration_distance(I_res,m_res_strategy);
