function [thetaF,thetaH, d, P, Q, R, n] = idealized(thetaF,thetaH, d, P, Q, R, n) %ASSUME: % -NO STRESS: % MxF: M=0.18, F=0.82, H=0.0 % MxH: M=1 % HxH: M=0.087, F=0.913, H=0.00 % -STRESS: % MxF: M=0.18, H=0.82 % MxH: M=1 % HxH: M=0.0136, F=0.0, H=0.863 %thetaF = probability that a female is without a male mate %thetaH = probability that a hermaphrodite is w/o a male %d = fitness cost for selfed offspring %s = presence of stress %0=no stress %1=stress %P = frequency of females %Q = frequency of hermaphrodites %R = frequency of males P0 = P; Q0 = Q; R0 = R; p = P; q = Q; r = R; %% no stress A = [p q r]; P0=p; Q0=q; R0=r; for i=1:1:n s=0; d=0; w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; %add frequency for each generation A = [A; P Q R]; P0=P; Q0=Q; R0=R; end %% intermediate stress B = [p q r]; P0=p; Q0=q; R0=r; for i=1:1:n s=0.5; w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; %add frequency for each generation B = [B; P Q R]; P0=P; Q0=Q; R0=R; end %% constant stress X = [p q r]; P0=p; Q0=q; R0=r; for i=1:1:n s=1; w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; %add frequency for each generation X = [X; P Q R]; P0=P; Q0=Q; R0=R; end %% binary stress for generations Y = [p q r]; P0=p; Q0=q; R0=r; for i=1:1:n s = randi([0,1],1,1); %random 0 or 1 w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; %add frequency for each generation Y = [Y; P Q R]; P0=P; Q0=Q; R0=R; end %% scale stress for generations E = [p q r]; P0=p; Q0=q; R0=r; for i=1:1:n s = rand; w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; %add frequency for each generation E = [E; P Q R]; P0=P; Q0=Q; R0=R; end %% new location (hold stress for a few generations) F = [p q r]; P0=p; Q0=q; R0=r; ST=[]; for i=1:1:20 Time = randi([7 9],1,1); %how many generations it stays s=rand*rand; for j=1:Time w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; %add frequency for each generation F = [F; P Q R]; P0=P; Q0=Q; R0=R; ST=[ST; s]; s=s+abs(1-s)*rand*rand; if s> 0.7 break end end end %% linear relationship of stress %stress from 0 to 1 C = []; P0 =p; Q0 =q; R0 =r; for i=0:0.1:1 s = i; w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; %add frequency for each generation C = [C; P Q R]; end %% di to tri Z = []; P0=p; Q0=q; R0=r; for i=1:1:n/2 Q0=0; s=rand; w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; Q = 0; Z = [Z; P Q R]; P0=P; Q0=Q; R0=R; end P Q R for i=n/2+1:1:n s = rand; %random 0 to 1 w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; %add frequency for each generation Z = [Z; P Q R]; P0=P; Q0=Q; R0=R; end %% scale stress for generations P0=0; Q0=1; R0=0; A2 = [P0 Q0 R0]; for i=1:1:n s = rand; w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; %add frequency for each generation A2 = [A2; P Q R]; P0=P; Q0=Q; R0=R; end %% fluctuating stress -> no stress G = []; P0=p; Q0=q; R0=r; for i=1:1:n/2 Q0=0; s=rand; w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; G = [G; P Q R]; P0=P; Q0=Q; R0=R; end for i=n/2+1:1:n s = 0.5; w = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)+ s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)+ ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0); P = (1-s)*(0.82*(1-thetaF)*P0+0.913*(1-d)*thetaH*Q0)/w; Q = s*(0.82*(1-thetaF)*P0+0.863*(1-d)*thetaH*Q0)/w; R = ((1-thetaH)*Q0+0.18*(1-thetaF)*P0+0.087*(1-s)*(1-d)*thetaH*Q0+0.136*s*(1-d)*thetaH*Q0)/w; %add frequency for each generation G = [G; P Q R]; P0=P; Q0=Q; R0=R; end %% %no stress figure(1) plot(A, 'linewidth', 2); T=title('No Stress'); set(T, 'FontSize', 20) x=xlabel('Generations'); set(x, 'FontSize', 20); xticks([0 n/2 n]); y=ylabel('Frequency'); yticks([0 0.5 1]); set(y, 'FontSize', 20); l=legend('females', 'hermaphrodites', 'males'); set(l, 'FontSize', 18); set(gca, 'FontSize', 18); ylim([0 1]); xlim([0 n]); %constant stress figure(2) plot(X, 'linewidth', 2); T=title('Constant Stress'); set(T, 'FontSize', 20); x=xlabel('Generations'); set(x, 'FontSize', 20); xticks([0 n/2 n]); y=ylabel('Frequency'); yticks([0 0.5 1]); set(y, 'FontSize', 20); l=legend('females', 'hermaphrodites', 'males'); set(l, 'FontSize', 18); set(gca, 'FontSize', 18); ylim([0 1]); xlim([0 n]); %intermediate stress figure(3) plot(B, 'linewidth', 2); T=title('Intermediate Stress'); set(T, 'FontSize', 20); x=xlabel('Generations'); set(x, 'FontSize', 20); xticks([0 n/2 n]); y=ylabel('Frequency'); yticks([0 0.5 1]); set(y, 'FontSize', 20); l=legend('females', 'hermaphrodites', 'males'); set(l, 'FontSize', 18); set(gca, 'FontSize', 18); ylim([0 1]); xlim([0 n]); %scaled stress figure(4) plot(E, 'linewidth', 2); T=title('Fluctuating Stress'); set(T, 'FontSize', 20); x=xlabel('Generations'); set(x, 'FontSize', 20); xticks([0 n/2 n]); y=ylabel('Frequency'); yticks([0 0.5 1]); set(y, 'FontSize', 20); l=legend('females', 'hermaphrodites', 'males'); set(l, 'FontSize', 18); set(gca, 'FontSize', 18); ylim([0 1]); xlim([0 n]); %linear relationship of stress figure(5) plot(C, 'linewidth', 2); T=title('Linear Relationship of Stress'); set(T, 'FontSize', 20); x=xlabel('Stress'); set(x, 'FontSize', 20); y=ylabel('Frequency'); set(y, 'FontSize', 20); l=legend('females', 'hermaphrodites', 'males'); set(l, 'FontSize', 18); ticks = [0:0.1:1]; xticklabels(ticks); set(gca, 'FontSize', 18); ylim([0 1]); xlim([1 11]); %with moving figure(6) plot(F, 'linewidth', 1.5); T=title('With Movement to New Locations'); set(T, 'FontSize', 20); x=xlabel('Generations'); set(x, 'FontSize', 20); y=ylabel('Frequency'); set(y, 'FontSize', 20); l=legend('females', 'hermaphrodites', 'males'); set(l, 'FontSize', 17); set(gca, 'FontSize', 18); xlim([0 50]); ylim([0 1]); yticks([0 0.5 1]); xlim([0 n]); %stress w/ movement figure(7) plot(ST, 'linewidth', 1.5); T=title('Stress With Movement to New Locations'); set(T, 'FontSize', 20); x=xlabel('Generations'); set(x, 'FontSize', 20); y=ylabel('Stress'); set(y, 'FontSize', 20); set(gca, 'FontSize', 18); yticks([0 0.5 1]); ylim([0 1]); xlim([0 n]); %dioecy -> trioecy figure(8) plot(Z, 'linewidth', 2); T=title('dioecy -> trioecy'); set(T, 'FontSize', 20); x=xlabel('Generations'); set(x, 'FontSize', 20); xticks([0 n/2 n]); y=ylabel('Frequency'); yticks([0 0.5 1]); set(y, 'FontSize', 20); l=legend('females', 'hermaphrodites', 'males'); set(l, 'FontSize', 18); set(gca, 'FontSize', 18); ylim([0 1]); xlim([0 n]); %Q0=1 figure(9) plot(A2, 'linewidth', 2); T=title('Single Hermaphrodite Colonization'); set(T, 'FontSize', 20); x=xlabel('Generations'); set(x, 'FontSize', 20); xticks([0 n/2 n]); y=ylabel('Frequency'); yticks([0 0.5 1]); set(y, 'FontSize', 20); l=legend('females', 'hermaphrodites', 'males'); set(l, 'FontSize', 18); set(gca, 'FontSize', 18); ylim([0 1]); xlim([0 n]); %Fluctuating stress to no stress figure(10) plot(G, 'linewidth', 2); T=title('Fluctuating to Intermediate Stress'); set(T, 'FontSize', 20); x=xlabel('Generations'); set(x, 'FontSize', 20); xticks([0 n/2 n]); y=ylabel('Frequency'); yticks([0 0.5 1]); set(y, 'FontSize', 20); l=legend('females', 'hermaphrodites', 'males'); set(l, 'FontSize', 18); set(gca, 'FontSize', 18); ylim([0 1]); xlim([0 n]); end