function [fitness_gradients,stand_times, P_data_out, Z_data_out, F_data_out, W_data_out, ... xP_data_out, xP_var_data_out, xF_data_out, xF_var_data_out, xZ_data_out, xZ_var_data_out, xW_data_out, xW_var_data_out,... t1,y1,final_gem_no_TC,final_gem_yes_TC,final_ode_no_TC,final_ode_yes_TC] = GEM_4_trophic_level_extinction_final(starting_mass,y0_4,a_scalings,h_scalings,a_prefactors,h_prefactors,h_targets,d0,dexp,r0,rexp,k,GGE,cv,h_2,num_replicates,t_max,run,EXT,m) % reconstitute parameters from passed vectors a0Z = a_prefactors(1); a0F = a_prefactors(2); a0W = a_prefactors(3); h0Z = h_prefactors(1); h0F = h_prefactors(2); h0W = h_prefactors(3); a1PM = a_scalings(1); a1ZM = a_scalings(2); a2ZM = a_scalings(3); a2FM = a_scalings(4); a3FM = a_scalings(5); a3WM = a_scalings(6); h1PM = h_scalings(1); h1ZM = h_scalings(2); h2ZM = h_scalings(3); h2FM = h_scalings(4); h3FM = h_scalings(5); h3WM = h_scalings(6); targethZonP = h_targets(1); targethFonZ = h_targets(2); targethWonF = h_targets(3); % determine rep number and standard times stand_times = 0:t_max/20:t_max; % standardized time steps for storing time series num_time_steps = length(stand_times); % preallocate fixed time matrices for abundances P_stand = nan(num_replicates,num_time_steps); % initial phytoplankton population size Z_stand = nan(num_replicates,num_time_steps); % initial zooplankton population size F_stand = nan(num_replicates,num_time_steps); % initial forage fish population size W_stand = nan(num_replicates,num_time_steps); % initial pisc fish population size % define initial population sizes P_stand(1:num_replicates,1) = y0_4(1); % initial phytoplankton population size Z_stand(1:num_replicates,1) = y0_4(2); % initial zooplankton population size F_stand(1:num_replicates,1) = y0_4(3); % initial forage fish population size W_stand(1:num_replicates,1) = y0_4(4); % initial pisc fish population size % preallocate fixed time matrices for traits and trait variances xP_stand = nan(num_replicates,num_time_steps); xZ_stand = nan(num_replicates,num_time_steps); xF_stand = nan(num_replicates,num_time_steps); xW_stand = nan(num_replicates,num_time_steps); xP_var_stand = nan(num_replicates,num_time_steps); xZ_var_stand = nan(num_replicates,num_time_steps); xF_var_stand = nan(num_replicates,num_time_steps); xW_var_stand = nan(num_replicates,num_time_steps); % put initial traits and variances into first elements xP_stand(1:num_replicates,1) = starting_mass(1); % initial mean trait xZ_stand(1:num_replicates,1) = starting_mass(2); xF_stand(1:num_replicates,1) = starting_mass(3); xW_stand(1:num_replicates,1) = starting_mass(4); xP_var_stand(1:num_replicates,1) = (cv*starting_mass(1))^2; xZ_var_stand(1:num_replicates,1) = (cv*starting_mass(2))^2; xF_var_stand(1:num_replicates,1) = (cv*starting_mass(3))^2; xW_var_stand(1:num_replicates,1) = (cv*starting_mass(4))^2; for i = 1:num_replicates % start Gillespie algorithm t = 0; % initial time i rng('shuffle'); % change random number seed % pre-allocate for the sliced standardized variables P_slice = nan(1,num_time_steps); Z_slice = nan(1,num_time_steps); F_slice = nan(1,num_time_steps); W_slice = nan(1,num_time_steps); xP_slice = nan(1,num_time_steps); xZ_slice = nan(1,num_time_steps); xF_slice = nan(1,num_time_steps); xW_slice = nan(1,num_time_steps); xP_var_slice = nan(1,num_time_steps); xZ_var_slice = nan(1,num_time_steps); xF_var_slice = nan(1,num_time_steps); xW_var_slice = nan(1,num_time_steps); % assign initial population sizes to sliced abundance variables P_slice(1) = y0_4(1); Z_slice(1) = y0_4(2); F_slice(1) = y0_4(3); W_slice(1) = y0_4(4); % also assign initial pop sizes to reduction variables (updated abundances) P = y0_4(1); Z = y0_4(2); F = y0_4(3); W = y0_4(4); % determine initial trait distributions MU_P = log(starting_mass(1)^2 / sqrt((cv*starting_mass(1))^2+starting_mass(1)^2)); % mean for lognormal SIGMA_P = sqrt(log((cv*starting_mass(1))^2/starting_mass(1)^2 + 1)); % std for lognormal xP_dist_init = lognrnd(MU_P,SIGMA_P,y0_4(1),1); % specify initial distribution of traits xP_dist = xP_dist_init; % reset trait distribution at the start of each simulation xP_slice(1) = mean(xP_dist); % pre-allocate for the ith slice of P_stand xP_var_slice(1) = var(xP_dist); % pre-allocate for the ith slice of P_stand MU_Z = log(starting_mass(2)^2 / sqrt((cv*starting_mass(2))^2+starting_mass(2)^2)); % mean for lognormal SIGMA_Z = sqrt(log((cv*starting_mass(2))^2/starting_mass(2)^2 + 1)); % std for lognormal xZ_dist_init = lognrnd(MU_Z,SIGMA_Z,y0_4(2),1); % specify initial distribution of traits xZ_dist = xZ_dist_init; % reset trait distribution at the start of each simulation xZ_slice(1) = mean(xZ_dist); % pre-allocate for the ith slice of P_stand xZ_var_slice(1) = var(xZ_dist); % pre-allocate for the ith slice of P_stand MU_F = log(starting_mass(3)^2 / sqrt((cv*starting_mass(3))^2+starting_mass(3)^2)); % mean for lognormal SIGMA_F = sqrt(log((cv*starting_mass(3))^2/starting_mass(3)^2 + 1)); % std for lognormal xF_dist_init = lognrnd(MU_F,SIGMA_F,y0_4(3),1); % specify initial distribution of traits xF_dist = xF_dist_init; % reset trait distribution at the start of each simulation xF_slice(1) = mean(xF_dist); % pre-allocate for the ith slice of P_stand xF_var_slice(1) = var(xF_dist); % pre-allocate for the ith slice of P_stand MU_W = log(starting_mass(4)^2 / sqrt((cv*starting_mass(4))^2+starting_mass(4)^2)); % mean for lognormal SIGMA_W = sqrt(log((cv*starting_mass(4))^2/starting_mass(4)^2 + 1)); % std for lognormal xW_dist_init = lognrnd(MU_W,SIGMA_W,y0_4(4),1); % specify initial distribution of traits xW_dist = xW_dist_init; % reset trait distribution at the start of each simulation xW_slice(1) = mean(xW_dist); % pre-allocate for the ith slice of P_stand xW_var_slice(1) = var(xW_dist); % pre-allocate for the ith slice of P_stand count = 1; % start counter to index steps while inside loop time_step_index = 2; % start a counter for standard times time_step = stand_times(time_step_index); % assign first standard time step while t < t_max if t >= EXT W = 0; xW_dist = []; end if P > 0 % as long as population size is > 0, pick another individual RP = randi(length(xP_dist),1); % randomly choose individual from the vector xP_next = xP_dist(RP); % pick the trait for that individual end if Z > 0 % as long as population size is > 0, pick another individual RZ = randi(length(xZ_dist),1); % randomly choose individual from the vector xZ_next = xZ_dist(RZ); % pick the trait for that individual end if F > 0 % as long as population size is > 0, pick another individual RF = randi(length(xF_dist),1); % randomly choose individual from the vector xF_next = xF_dist(RF); % pick the trait for that individual end if W > 0 % as long as population size is > 0, pick another individual RW = randi(length(xW_dist),1); % randomly choose individual from the vector xW_next = xW_dist(RW); % pick the trait for that individual end h_Z_next = h_prefactors(1).*mean(xP_dist).^h_scalings(1).*xZ_next.^h_scalings(2); % use pop average h h_F_next = h_prefactors(2).*mean(xZ_dist).^h_scalings(3).*xF_next.^h_scalings(4); % use pop average h h_W_next = h_prefactors(3).*mean(xF_dist).^h_scalings(5).*xW_next.^h_scalings(6); % use pop average h r_next = r0*xP_next^rexp; % estimate r from phyto mass a_Z_next = a0Z*xP_next^a1PM*xZ_next^a1ZM; % calculate a from phyto and zoop mass a_F_next = a0F*xZ_next^a2ZM*xF_next^a2FM; % turn current Zoo mass into a2 a_W_next = a0W*xF_next^a3FM*xW_next^a3WM; % turn current Fish mass into a3 d_W_next = d0*xW_next^dexp; d_W_next = (d_W_next - d_W_next*F/(F+1000)); d_F_next = d0*xF_next^dexp; d_F_next = (d_F_next - d_F_next*Z/(Z+3)); d_Z_next = d0*xZ_next^dexp; d_Z_next = (d_Z_next - d_Z_next*P/(P+10)); e_Z_next = GGE*xP_next/xZ_next; % e for Z e_F_next = GGE*xZ_next/xF_next; % e for F e_W_next = GGE*xF_next/xW_next; % e for W % set up rates of each possible event %%%%% Phytoplankton % Birth rate of Phytoplankton b_P = r_next*P; % Natural death rate of Phytoplankton d_P_1 = r_next*P^2/k; % otherwise r and k are fixed % Phytoplankton mortality rate from predation by zooplankton d_P_2 = a_Z_next*P*(Z/1e3)^(1+m)/(1+a_Z_next*h_Z_next*P*Z^m); %%%%% Zooplankton % Zooplankton birth rate b_Z = 1e3*e_Z_next*d_P_2; % Zooplankton death rate d_Z_1 = d_Z_next*Z; % Zooplankton mortality rate from predation by Fish d_Z_2 = a_F_next*Z*(F/1e6)^(1+m)/(1+a_F_next*h_F_next*(Z/1e3)*F^m); %%%%% Fish % Fish birth rate b_F = 1e3*e_F_next*d_Z_2; % Natural death rate of Fish d_F_1 = d_F_next*F; % Fish mortality rate from predation by Pisc Fish d_F_2 = a_W_next*F*(W/1e9)^(1+m)/(1+a_W_next*h_W_next*(F/1e6)*W^m); %%%%% Piscivorous Fish % Pisc Fish birth rate (convert predated Fish into Pisc Fish) b_W = 1e3*e_W_next*d_F_2; % Natural death rate of Pisc Fish d_W_1 = d_W_next*W; % sum the events to make wheel of fortune CS_vector = cumsum([b_P d_P_1 d_P_2 b_Z d_Z_1 d_Z_2 b_F d_F_1 d_F_2 b_W d_W_1]); Slice_widths = CS_vector./CS_vector(end); LI = rand < Slice_widths; Event_index = find(LI,1,'first'); % now choose actual events % Phytoplankton Birth if Event_index == 1 xP_parent = (1-h_2)*mean(log(xP_dist)) + h_2*log(xP_next); % weight the expected value of the offspring between the initial mean offP_std = sqrt(1-h_2^2)*((1-h_2)*std(log(xP_dist_init))+h_2*std(log(xP_dist))); xP_dist(length(xP_dist)+1) = exp(pearsrnd(xP_parent,offP_std,0,3)); % Phytoplankton natural death elseif Event_index == 2 xP_dist(RP) = []; % reduce dist by lost individual % Phytoplankton predation death from Zooplankton elseif Event_index == 3 xP_dist(RP) = []; % reduce dist by lost individual % Zooplankton birth elseif Event_index == 4 xZ_parent = (1-h_2)*mean(xZ_dist) + h_2*xZ_next; % weight the expected value of the offspring between the initial mean % and the parent mean, with h_2 as the weighting offZ_std = sqrt(1-h_2^2)*((1-h_2)*std(xZ_dist_init)+h_2*std(xZ_dist)); MU_Z = log(xZ_parent^2 / sqrt(offZ_std^2+xZ_parent^2)); SIGMA_Z = sqrt(log(offZ_std^2/xZ_parent^2 + 1)); xZ_dist(length(xZ_dist)+1) = lognrnd(MU_Z,SIGMA_Z,1,1); % pick offspring from lognormal dist and add to trait distribution % Zooplankton natural death elseif Event_index == 5 xZ_dist(RZ) = []; % reduce dist by lost individual % Zooplankton predation death from Fish elseif Event_index == 6 xZ_dist(RZ) = []; % reduce dist by lost individual % Fish birth elseif Event_index == 7 xF_parent = (1-h_2)*mean(xF_dist) + h_2*xF_next; % weight the expected value of the offspring between the initial mean % and the parent mean, with h_2 as the weighting offF_std = sqrt(1-h_2^2)*((1-h_2)*std(xF_dist_init)+h_2*std(xF_dist)); MU_F = log(xF_parent^2 / sqrt(offF_std^2+xF_parent^2)); SIGMA_F = sqrt(log(offF_std^2/xF_parent^2 + 1)); xF_dist(length(xF_dist)+1) = lognrnd(MU_F,SIGMA_F,1,1); % pick offspring from lognormal dist and add to trait distribution % Fish natural death elseif Event_index == 8 xF_dist(RF) = []; % reduce dist by lost individual % Fish predation death from Piscivore Fish elseif Event_index == 9 xF_dist(RF) = []; % reduce dist by lost individual % Piscivorous Fish birth elseif Event_index == 10 xW_parent = (1-h_2)*mean(xW_dist) + h_2*xW_next; % weight the expected value of the offspring between the initial mean % and the parent mean, with h_2 as the weighting offW_std = sqrt(1-h_2^2)*((1-h_2)*std(xW_dist_init)+h_2*std(xW_dist)); MU_W = log(xW_parent^2 / sqrt(offW_std^2+xW_parent^2)); SIGMA_W = sqrt(log(offW_std^2/xW_parent^2 + 1)); xW_dist(length(xW_dist)+1) = lognrnd(MU_W,SIGMA_W,1,1); % pick offspring from lognormal dist and add to trait distribution % Piscivorous Fish natural death elseif Event_index == 11 xW_dist(RW) = []; % reduce dist by lost individual end P = length(xP_dist); % add up phytos Z = length(xZ_dist); % took away a Zooplankton F = length(xF_dist); % removed a Fish W = length(xW_dist); % removed a Pisc Fish if t > time_step P_slice(time_step_index) = P; % assign current values to sliced standard times Z_slice(time_step_index) = Z; F_slice(time_step_index) = F; W_slice(time_step_index) = W; xP_slice(time_step_index) = mean(xP_dist); xZ_slice(time_step_index) = mean(xZ_dist); xF_slice(time_step_index) = mean(xF_dist); xW_slice(time_step_index) = mean(xW_dist); xP_var_slice(time_step_index) = var(xP_dist); xZ_var_slice(time_step_index) = var(xZ_dist); xF_var_slice(time_step_index) = var(xF_dist); xW_var_slice(time_step_index) = var(xW_dist); time_step_index = time_step_index+1; time_step = stand_times(time_step_index); end t = t + exp(-1/CS_vector(end))/(CS_vector(end)); count = count+1; end % pass in current values to the final standardized time P_slice(end) = P; Z_slice(end) = Z; F_slice(end) = F; W_slice(end) = W; xP_slice(end) = mean(xP_dist); xZ_slice(end) = mean(xZ_dist); xF_slice(end) = mean(xF_dist); xW_slice(end) = mean(xW_dist); xP_var_slice(end) = var(xP_dist); xZ_var_slice(end) = var(xZ_dist); xF_var_slice(time_step_index) = var(xF_dist); xW_var_slice(end) = var(xW_dist); % pass sliced variable to standard time matrix P_stand(i,:) = P_slice; Z_stand(i,:) = Z_slice; F_stand(i,:) = F_slice; W_stand(i,:) = W_slice; xP_stand(i,:) = xP_slice; xZ_stand(i,:) = xZ_slice; xF_stand(i,:) = xF_slice; xW_stand(i,:) = xW_slice; xP_var_stand(i,:) = xP_var_slice; xZ_var_stand(i,:) = xZ_var_slice; xF_var_stand(i,:) = xF_var_slice; xW_var_stand(i,:) = xW_var_slice; end % calculate ci's for time series upper_ci_level = 75; % choose ci levels lower_ci_level = 25; % choose ci levels % Phyto abundance test(:,:) = P_stand(:,:); ci_P_up = prctile(test,lower_ci_level); ci_P_down = prctile(test,upper_ci_level); median_P = prctile(test,50); P_data_out = [median_P; ci_P_up; ci_P_down]; % Zoo abundance test(:,:) = Z_stand(:,:); ci_Z_up = prctile(test,lower_ci_level); ci_Z_down = prctile(test,upper_ci_level); median_Z = prctile(test,50); Z_data_out = [median_Z; ci_Z_up; ci_Z_down]; % Fish abundance test(:,:) = F_stand(:,:); ci_F_up = prctile(test,lower_ci_level); ci_F_down = prctile(test,upper_ci_level); median_F = prctile(test,50); F_data_out = [median_F; ci_F_up; ci_F_down]; % Pisc Fish abundance test(:,:) = W_stand(:,:); ci_W_up = prctile(test,lower_ci_level); ci_W_down = prctile(test,upper_ci_level); median_W = prctile(test,50); W_data_out = [median_W; ci_W_up; ci_W_down]; % Phytoplankton median mass through time test(:,:) = xP_stand(:,:); ci_xP_up = prctile(test,lower_ci_level); ci_xP_down = prctile(test,upper_ci_level); median_xP = prctile(test,50); xP_data_out = [median_xP; ci_xP_up; ci_xP_down]; % Phytoplankton variance in mass through time test(:,:) = xP_var_stand(:,:); ci_xP_var_up = prctile(test,lower_ci_level); ci_xP_var_down = prctile(test,upper_ci_level); median_xP_var = prctile(test,50); xP_var_data_out = [median_xP_var; ci_xP_var_up; ci_xP_var_down]; % Zooplankton median mass through time test(:,:) = xZ_stand(:,:); ci_xZ_up = prctile(test,lower_ci_level); ci_xZ_down = prctile(test,upper_ci_level); median_xZ = prctile(test,50); xZ_data_out = [median_xZ; ci_xZ_up; ci_xZ_down]; % Zooplankton variance in mass through time test(:,:) = xZ_var_stand(:,:); ci_xZ_var_up = prctile(test,lower_ci_level); ci_xZ_var_down = prctile(test,upper_ci_level); median_xZ_var = prctile(test,50); xZ_var_data_out = [median_xZ_var; ci_xZ_var_up; ci_xZ_var_down]; % Fish median mass through time test(:,:) = xF_stand(:,:); ci_xF_up = prctile(test,lower_ci_level); ci_xF_down = prctile(test,upper_ci_level); median_xF = prctile(test,50); xF_data_out = [median_xF; ci_xF_up; ci_xF_down]; % Fish variance in mass through time test(:,:) = xF_var_stand(:,:); ci_xF_var_up = prctile(test,lower_ci_level); ci_xF_var_down = prctile(test,upper_ci_level); median_xF_var = prctile(test,50); xF_var_data_out = [median_xF_var; ci_xF_var_up; ci_xF_var_down]; % Piscivorous fish median mass through time test(:,:) = xW_stand(:,:); ci_xW_up = prctile(test,lower_ci_level); ci_xW_down = prctile(test,upper_ci_level); median_xW = prctile(test,50); xW_data_out = [median_xW; ci_xW_up; ci_xW_down]; % Piscivorous fish variance in mass through time test(:,:) = xW_var_stand(:,:); ci_xW_var_up = prctile(test,lower_ci_level); ci_xW_var_down = prctile(test,upper_ci_level); median_xW_var = prctile(test,50); xW_var_data_out = [median_xW_var; ci_xW_var_up; ci_xW_var_down]; % call function for fitness gradient through time P_fit_grad = fitness_gradient_phytoplankton(r0,median_xP,rexp,median_P,k,a0Z,a1PM,median_xZ,a1ZM,median_Z,h0Z); Z_fit_grad = fitness_gradient_zooplankton(median_xP,a0Z,a1PM,median_xZ,a1ZM,median_P,median_Z,median_F,GGE,a0F,a2ZM,median_xF,a2FM,d0,dexp,h2ZM,h2FM,h0F,h1PM,h1ZM,h0Z); % adjust pisc fish states to prevent NaN after extinction median_xW_alt = median_xW; median_xW_alt(isnan(median_xW)==1) = 1; median_W_alt = median_W; median_W_alt(median_W_alt == 0) = 1; F_fit_grad = fitness_gradient_fish(median_F,a0F,a2ZM,GGE,median_xF,a2FM,median_xZ,median_Z,h0F,median_W_alt,a3FM,a0W,median_xW_alt,a3WM,h0W,d0,dexp,h3WM,h3FM,h2FM,h2ZM); W_fit_grad = fitness_gradient_walleye(median_F,GGE,median_W,a0W,a3WM,median_xF,a3FM,median_xW,h0W,d0,dexp,h3FM,h3WM); fitness_gradients = [P_fit_grad; Z_fit_grad; F_fit_grad; W_fit_grad]; % realized initial parameters r = r0*median_xP(1)^rexp; eZ = GGE*median_xP(1)/median_xZ(1); % e for Z eF = GGE*median_xZ(1)/median_xF(1); % e for F eW = GGE*median_xF(1)/median_xW(1); % e for W aZ = a0Z*(median_xZ(1)^a1ZM)*(median_xP(1)^a1PM); aF = a0F*(median_xF(1)^a2FM)*(median_xZ(1)^a2ZM); aW = a0W*(median_xW(1)^a3WM)*(median_xF(1)^a3FM); hZ = h0Z*(median_xZ(1)^h1ZM)*(median_xP(1)^h1PM); hF = h0F*(median_xF(1)^h2FM)*(median_xZ(1)^h2ZM); hW = h0W*(median_xW(1)^h3WM)*(median_xF(1)^h3FM); dZ = d0*median_xZ(1)^dexp; dF = d0*median_xF(1)^dexp; dW = d0*median_xW(1)^dexp; ode_4 = @(t,y) MR_model_four_vol_V5(t,y,r,k,aZ,aF,aW,hZ,hF,hW,eZ,eF,eW,dZ,dF,dW); % compile function and call [t1,y1] = ode45(ode_4, [0 t_max], y0_4); % return time and population density vectors % rerun solution starting at time = EXT tmp = abs(t1-EXT); minima = min(tmp); idx = find(tmp==minima); if run == 2 % removed top predator scenario [t1TC,y1TC] = ode45(ode_4, [EXT t_max], [y1(idx,1:3) 0]); % return time and population density vectors % combine ODE solutions < and > loss of top predator t1 = [t1(t1