function [prop_cooperators_traj, modifiability_traj, correlation_coeff_traj]=multilocus_model(timesteps,popsize,mu,m,b,d,c,allowmodification) % Grafen ESS vs standard mixed ESS % there are 1+2N loci, stored in 2+2N columns: all but the first locus are % 'null' in the beginning. So, in the beginning it works like in Grafen's % discrete case % the first locus has two types of alleles: 0 for defect and 1 for % cooperate these alleles are initially fully modifiable, but can mutate to % become modifiable to various extents - a property encoded in the second % column as a number between 0 and 2 (with a mean of 1 under random drift). % At values >1, modifiers have their full weight. But at <1, modifiers are % weighted by the value of the 'modifiability' variable. % then there are N modifiers that can be 0 or 1, and whose sum/N is the % probability of reversing the Grafenian default 'defect' to 'cooperate' % there are also N modifiers for the opposite switch % input %timesteps=2000; trajectories=1; % record data for every time step %popsize=5000; % number of individuals (should be even) %mu=0.001; % mutation rate that applies to all loci mu1=mu; % mutation rate at first locus mu2=mu; % mutation rate at modifier loci %m=3; % number of modifier loci per phenotype penetrance=1; % probability that the modifier loci, if all "on", can collectively reverse the phenotype specified by locus 1 baseline=(-1).*min([b+d-c,-c,b,0]); % to avoid negative payoffs %b=5; % additive payoff component from cooperating %d=-2.5; % non-additive payoff component from mutual cooperation %c=1; % cost of cooperating fecundity=4; % should be even %allowmodification=1; % shall the first locus be modifiable? 0= no; 1= yes, 100%; 2= yes, gradually, according to a continuous, evolvable property of the first locus (in the range 0..2, making 1 the mean under drift) inds=zeros(popsize,1); % initialise female matrix lastgene=2+2*m; % just a shorthand for the column holding the last gene phenotype=lastgene+1; % the column storing an individual's phenotype; i.e. its probability to cooperate payoff=lastgene+2; % the column storing an individual's payoff from interacting with a sibling totaloffspring=popsize/2*fecundity; % initialise inds(:,1)=randi([0 1],popsize,1); % the first locus, with "0" for defect and "1" for cooperate inds(:,2)=(allowmodification==0).*zeros(popsize,1) + (allowmodification==1) .* ones(popsize,1) + (allowmodification==2).*rand(popsize,1)*2; % inds(:,3:lastgene)=zeros(length(inds),2*m); % initialise 2*loci modifier genes per individual; the first half are for switching the defector gene to 'cooperation; the second half are for the reverse switch if trajectories==1 % store means over time prop_cooperators_traj=zeros(1,timesteps); modifiability_traj=zeros(1,timesteps); correlation_coeff_traj=zeros(1,timesteps); else prop_cooperators_traj=0; modifiability_traj=0; correlation_coeff_traj=0; end for t=1:timesteps % if t/1000==ceil(t/1000), t %#ok % end sires=inds(1:popsize/2,:); % the first half are assigned as "males" - an arbitrary label, as we assume no sex differences dams=inds(popsize/2+1:popsize,:); % the second half are assigned as "females" % inheritance maternalattributes=kron(dams,ones(fecundity,1)); % a list of offspring, with colums showing their mothers' attributes paternalattributes=kron(sires,ones(fecundity,1)); throwdice=rand(length(maternalattributes(:,1)),lastgene); % for each locus, decide randomly between maternal and paternal inheritance offspring=(throwdice<0.5).*maternalattributes(:,1:lastgene)+~(throwdice<0.5).*paternalattributes(:,1:lastgene); % assign to offspring what genes they inherited from which parent offspring(:,2)=(throwdice(:,1)<0.5).*maternalattributes(:,2)+~(throwdice(:,1)<0.5).*paternalattributes(:,2); % in the second column, assign modifiability corresponding to the allele inherited at the first locus defection_modification_score =sum(offspring(:,3:2+m),2)*penetrance/m; % total score of first half of modifier loci cooperation_modification_score=sum(offspring(:,2+m+1:lastgene),2)*penetrance/m; % total score of second half of modifier loci modifiability= min(offspring(:,2),1 ); % At values >1, modifiers have their full weight. But at <1, modifiers are weighted by the value of the 'modifiability' variable. cooperation_propensity= (offspring(:,1)==0) .* (modifiability.*defection_modification_score) + (offspring(:,1)==1) .* (offspring(:,1)-modifiability.*cooperation_modification_score); % calculate propensity to cooperate offspring(:,phenotype)=cooperation_propensity>rand(size(cooperation_propensity)); % decide phenotype (1=cooperate) based on cooperation propensity odd_offspring = offspring(1:2:end,:); % odd-numbered offspring even_offspring= offspring(2:2:end,:); % even-numbered offspring [r,P]=corrtest(odd_offspring(:,phenotype),even_offspring(:,phenotype)); % test for phenotypic correlation correlation_coeff_traj(t)=r; % coefficient of phenotypic correlation % let the games begin odd_offspring(:,payoff) = baseline + ... odd_offspring(:,phenotype).* even_offspring(:,phenotype) *(b+d-c)... % personal payoff from mutual cooperation + odd_offspring(:,phenotype).* ~even_offspring(:,phenotype) *(-c)... % personal payoff from unilateral cooperation + ~odd_offspring(:,phenotype).* even_offspring(:,phenotype) *(b)... % personal payoff from unilateral defection + ~odd_offspring(:,phenotype).* ~even_offspring(:,phenotype) *(0); % personal payoff from mutual defection even_offspring(:,payoff) = baseline + ... even_offspring(:,phenotype).* odd_offspring(:,phenotype) *(b+d-c)... % personal payoff from mutual cooperation + even_offspring(:,phenotype).* ~odd_offspring(:,phenotype) *(-c)... % personal payoff from unilateral cooperation + ~even_offspring(:,phenotype).* odd_offspring(:,phenotype) *(b)... % personal payoff from unilateral defection + ~even_offspring(:,phenotype).* ~odd_offspring(:,phenotype) *(0); % personal payoff from mutual defection offspring=[odd_offspring; even_offspring]; % stack odd and even offspring to combine them in single matrix % generate candidate mutations in phenotype encoded by 1st locus mutations=randi([0 1],totaloffspring,1); % corresponding_modifiabilities=(allowmodification==0).*zeros(totaloffspring,1) + (allowmodification==1).*ones(totaloffspring,1) + (allowmodification==2).* rand(totaloffspring,1)*2; % these correspond to the potential mutations of the previous line % implement mutations mutate=rand(totaloffspring,1)