function mimicry_model(filename,link, nsteps); %%%%%%%%%%%%%%%%%%% % larval trait evolution model % % link options: % % 'nolink' - no linked traits % 'llink' - larval traits 1 & 2 linked (1 must be beneficial before 2 is) % 'plink-' - negative parent effects on positive larval traits,no linked larval traits % 'plink+' - positive parent effects on positive larval traits,no linked larval traits % 'lplink+' - positive parent effects on positive larval traits, linked larval traits % 'lplink+' - positive parent effects on positive larval traits, linked larval traits %%%%%%%%%%%%%%%%%%% tic n=1e6; % number of larvae m=1e3; % number of adults survival_factor=.0000005; mortality_rate=0.999995; mean_trait=zeros(nsteps,2); % first assume normal distribution of larval fish traits, say long spines % and bulbs on spines %define genes of initial potential parents for trait 1 & trait 2 % 2nd column is trait, 3rd column is genes for that trait parent_genes=randn(m,2,2); parent_age=ones(m,1); init_parent_genes=parent_genes; for i=1:nsteps %assign parent mortality parent_mort=ones(size(parent_age)).*0.025; %mate parents and assign genes to larvae larvae_genes=zeros(n,2,2); if i==1 parents=[ceil(m.*rand(n,1)), ceil(m.*rand(n,1))]; else parents=[ceil(length(parent_age).*rand(n,1)), ceil(length(parent_age).*rand(n,1))]; end %set larvae genes for trait 1 larvae_pr=rand(n,2); %random number to determine which of parent genes larvae gets larvae_genes(larvae_pr(:,1)<0.5,1,1)=parent_genes(parents(larvae_pr(:,1)<0.5,1),1,1); larvae_genes(larvae_pr(:,1)>=0.5,1,1)=parent_genes(parents(larvae_pr(:,1)>=0.5,1),1,2); larvae_genes(larvae_pr(:,2)<0.5,1,2)=parent_genes(parents(larvae_pr(:,2)<0.5,2),1,1); larvae_genes(larvae_pr(:,2)>=0.5,1,2)=parent_genes(parents(larvae_pr(:,2)>=0.5,2),1,2); %set larvae genes for trait 2 larvae_pr=rand(n,2); %random number to determine which of parent genes larvae gets larvae_genes(larvae_pr(:,1)<0.5,2,1)=parent_genes(parents(larvae_pr(:,1)<0.5,1),2,1); larvae_genes(larvae_pr(:,1)>=0.5,2,1)=parent_genes(parents(larvae_pr(:,1)>=0.5,1),2,2); larvae_genes(larvae_pr(:,2)<0.5,2,2)=parent_genes(parents(larvae_pr(:,2)<0.5,2),2,1); larvae_genes(larvae_pr(:,2)>=0.5,2,2)=parent_genes(parents(larvae_pr(:,2)>=0.5,2),2,2); %determine survival of larvae switch link case 'llink' %for larval linked trait model uncomment following two lines larv_trait=larvae_genes(:,1,1)+larvae_genes(:,1,2); larv_trait(larv_trait>4)=larv_trait(larv_trait>4)+larvae_genes(larv_trait>4,2,1)+larvae_genes(larv_trait>4,2,2); % benefit of trait 2 is dependent on a certain level of trait 1. case 'nolink' larv_trait=(larvae_genes(:,1,1)+larvae_genes(:,1,2)+larvae_genes(:,2,1)+larvae_genes(:,2,2)); %larv_trait=(larvae_genes(:,1,1)>1.65)+(larvae_genes(:,1,2)>1.65)+(larvae_genes(:,2,1)>1.65)+(larvae_genes(:,2,2)>1.65); case 'plink-' larv_trait=(larvae_genes(:,1,1)+larvae_genes(:,1,2)+larvae_genes(:,2,1)+larvae_genes(:,2,2)); %larv_trait=(larvae_genes(:,1,1)>1.65)+(larvae_genes(:,1,2)>1.65)+(larvae_genes(:,2,1)>1.65)+(larvae_genes(:,2,2)>1.65); parent_trait=10.*(parent_genes(:,1,1)+parent_genes(:,1,2)+parent_genes(:,2,1)+parent_genes(:,2,2)); parent_mort=parent_mort+parent_trait.*survival_factor; %harmful to parent case 'plink+' larv_trait=(larvae_genes(:,1,1)+larvae_genes(:,1,2)+larvae_genes(:,2,1)+larvae_genes(:,2,2)); %larv_trait=(larvae_genes(:,1,1)>1.65)+(larvae_genes(:,1,2)>1.65)+(larvae_genes(:,2,1)>1.65)+(larvae_genes(:,2,2)>1.65); parent_trait=10.*(parent_genes(:,1,1)+parent_genes(:,1,2)+parent_genes(:,2,1)+parent_genes(:,2,2)); parent_mort=parent_mort-parent_trait.*survival_factor; %beneficial to parent case 'lplink-' larv_trait=larvae_genes(:,1,1)+larvae_genes(:,1,2); larv_trait(larv_trait>4)=larv_trait(larv_trait>4)+larvae_genes(larv_trait>4,2,1)+larvae_genes(larv_trait>4,2,2); % benefit of trait 2 is dependent on a certain level of trait 1. parent_trait=10.*(parent_genes(:,1,1)+parent_genes(:,1,2)+parent_genes(:,2,1)+parent_genes(:,2,2)); parent_mort=parent_mort+parent_trait.*survival_factor; %harmful to parent case 'lplink+' larv_trait=larvae_genes(:,1,1)+larvae_genes(:,1,2); larv_trait(larv_trait>4)=larv_trait(larv_trait>4)+larvae_genes(larv_trait>4,2,1)+larvae_genes(larv_trait>4,2,2); % benefit of trait 2 is dependent on a certain level of trait 1. parent_trait=10.*(parent_genes(:,1,1)+parent_genes(:,1,2)+parent_genes(:,2,1)+parent_genes(:,2,2)); parent_mort=parent_mort-parent_trait.*survival_factor; %beneficial to parent end larv_mort=mortality_rate-survival_factor.*larv_trait; lsurvive_test=rand(n,1); larvae_genes=larvae_genes((lsurvive_test>=larv_mort)==1,:,:); %determine survival of parents if i==1 psurvive_test=rand(m,1); else psurvive_test=rand(length(parent_age),1); end %for constant age-dependent mortality uncomment following two lines. Mortality is 0.5 and increasese by 0.1 lineaerly after 8 years of age. parent_mort(parent_age>=8)=(parent_mort(parent_age>=8)+0.1.*(parent_age(parent_age>=8)-8)); %for parental linked trait model uncomment following line and one of the two after that lines %parent_trait=1000.*(parent_genes(:,1,1)+parent_genes(:,1,2)+parent_genes(:,2,1)+parent_genes(:,2,2)); %parent_mort=parent_mort+parent_trait.*survival_factor; %harmful to parent %parent_mort=parent_mort-parent_trait.*survival_factor; %beneficial to parent parent_genes=parent_genes(psurvive_test>=parent_mort,:,:); parent_age=parent_age(psurvive_test>=parent_mort); %determine new parent pool assuming all larvae survive to reproduce parent_genes=[parent_genes; larvae_genes]; parent_age=parent_age+1; parent_age=[parent_age; ones(size(larvae_genes,1),1)]; if(rem(i,floor(nsteps./100))==0) fprintf('%d%s%d%s\n',100.*i./nsteps,'% done in '); toc; end mean_trait(i,:)=mean(squeeze(mean(parent_genes))'); age_dist(i,:)=hist(parent_age,[1:25]); mean_pmort(i)=mean(parent_mort); mean_lmort(i)=mean(larv_mort); clear larvae_genes; end eval(['save -mat ',filename,'.mat parent_genes init_parent_genes mean_trait mean_pmort mean_lmort age_dist']); end