%% info % this is a simulation for public goods game on % complete graphs function compg50(par1,par2,par3) % par1 is alpha1 % par2 is alpha2 % par3 is lambda %% parameters N = 50; % population size = graph size alpha1 = par1; % staying propensity of A individual alpha2 = par2; % staying propensity of B individual beta1 = 1; % benefit of interacting with A individual beta2 = 1; % benefit of interacting with B individual S = 0.03; % sensitivity to group members c = 0.04; % cost of cooperation v = 0.4; % benefits of cooperation lambda = par3; % cost of movement T = 10; % number of movement steps before returning home trials = 100000; % how many simulations will be averaged %% setup alpha = [alpha1 alpha2]; % staying propensities beta = [beta1 beta2]; % benefits A = 1; % type A individual B = 2; % type B individual cooperator = 1; % notation for cooperator strategy defector = 2; % notation for defector strategy istrat = [cooperator cooperator]; % interactive strategies of A % and B individuals % initialize the matrix containing information about each individual % column 1 = type (A or B) % column 2 = current position % column 3 = accumulated fitness % row number corresponds to the home location of the individual IS = zeros(N,3); % create placeholders for column indices type = 1; position = 2; fitness = 3; E = ones(N) - eye(N); % edge matrix for complete graph % construct a cell array containing indices of locations for possible % movement from any given node in column 1 % and number of neighboring locations in column 2 % the number of neighbors for each node could be different Neighbors = cell(N,2); for k=1:N Neighbors{k,1} = find(E(k,:)==1); % locations connected to node k Neighbors{k,2} = size(Neighbors{k,1},2); % number of such locations end % initialize the matrix containing demographics info for each graph node % column 1 = current number of As at the given position % column 2 = current number of Bs at the given position % the matrix will get actual initial values % and will be updated after every event Node_stats = zeros(N,2); RWM = zeros(N,N,N-1); % replacement weight matrix placeholder % last index corresponds to the number of A individuals in the population % we assume that A individuals occupy first nodes in the graph % followed by B individuals % compute the replacement weight matrix analytically p1 = (1-alpha1)/(N-1); % probability A goes to a given place p2 = (1-alpha2)/(N-1); % probability B goes to a given place for k=1:N-1 % k denotes the number of A individuals in the population % we start by computing all probabilities of different types of groups alone1 = alpha1*(1-p1)^(k-1)*(1-p2)^(N-k)... + (k-1)/(N-1)*(1-alpha1)^2*(1-p1)^(k-2)*(1-p2)^(N-k)... + (N-k)/(N-1)*(1-alpha1)*(1-alpha2)*(1-p1)^(k-1)*(1-p2)^(N-k-1); % probability individual A stays alone if there are k individuals A alone2 = alpha2*(1-p2)^(N-k-1)*(1-p1)^k... + (N-k-1)/(N-1)*(1-alpha2)^2*(1-p2)^(N-k-2)*(1-p1)^k... + k/(N-1)*(1-alpha2)*(1-alpha1)*(1-p2)^(N-k-1)*(1-p1)^(k-1); % probability individual B stays alone if there are k individuals A AmeetsA = 2*alpha1*p1*(1-p1)^(k-2)*(1-p2)^(N-k) ... + (k>=3)*(k-2)*p1^2*(1-alpha1)*(1-p1)^(k-3)*(1-p2)^(N-k) ... + (N-k)*p1^2*(1-p1)^(k-2)*(1-alpha2)*(1-p2)^(N-k-1); % probability a focal A meets another given A in a group of size 2 % larger groups will be accounted for in the following cycle for m=3:N-1 % here m stands for the size of the group where two A's meet l = min([k-2 m-2]); % max number of other As in the group a = max([0 m-2+k-N]); % need at least so many other As in the group groupa = 0; % they meet at focal A home or given A home for j=a:l % j is the number of other A's in the group groupa = groupa + nchoosek(k-2,j)*p1^j*(1-p1)^(k-2-j) ... *nchoosek(N-k,m-2-j)*p2^(m-2-j)*(1-p2)^(N-k-m+2+j); end groupa = groupa*alpha1*p1*2; % focal A and given A l1 = min([k-3 m-3]); % max number of other As in the group a1 = max([0 m-3+k-N]); % need at least so many other As in the group groupc1 = 0; % they meet at another A's home and host stays for j=a1:l1 groupc1 = groupc1 + nchoosek(k-3,j)*p1^j*(1-p1)^(k-3-j) ... *nchoosek(N-k,m-3-j)*p2^(m-3-j)*(1-p2)^(N-k-m+3+j); end groupc1 = groupc1*alpha1; % the host stays l2 = min([k-3 m-2]); % max number of other As in the group a2 = max([0 m-2+k-N]); % need at least so many other As in the group groupc2 = 0; % they meet at another A's home and host leaves for j=a2:l2 groupc2 = groupc2 + nchoosek(k-3,j)*p1^j*(1-p1)^(k-3-j) ... *nchoosek(N-k,m-2-j)*p2^(m-2-j)*(1-p2)^(N-k-m+2+j); end groupc2 = groupc2*(1-alpha1); % the host leaves groupc = (k>=3)*(k-2)*p1^2*(groupc1+groupc2); % all possible groups of type c n1 = min([k-2 m-3]); % max number of other As in the group b1 = max([0 m-2+k-N]); % need at least so many other As in the group groupd1 = 0; % they meet at some B's home and host stays for j=b1:n1 groupd1 = groupd1 + nchoosek(k-2,j)*p1^j*(1-p1)^(k-2-j) ... *nchoosek(N-k-1,m-3-j)*p2^(m-3-j)*(1-p2)^(N-k-m+2+j); end groupd1 = groupd1*alpha2; % the host stays n2 = min([k-2 m-2]); % max number of other As in the group b2 = max([0 m-1+k-N]); % need at least so many other As in the group groupd2 = 0; % they meet at some B's home and host leaves for j=b2:n2 groupd2 = groupd2 + nchoosek(k-2,j)*p1^j*(1-p1)^(k-2-j) ... *nchoosek(N-k-1,m-2-j)*p2^(m-2-j)*(1-p2)^(N-k-m+1+j); end groupd2 = groupd2*(1-alpha2); % the host leaves groupd = (N-k)*p1^2*(groupd1+groupd2); % all possible groups of type d AmeetsA = AmeetsA + (groupa+groupc+groupd)/(m-1); % adjust time end % we went through all the possible group sizes from 3 to N-1 % it now remains to account the group of size N AmeetsA = AmeetsA + (k*alpha1*p1^(k-1)*p2^(N-k) ... + (N-k)*p1^k*alpha2*p2^(N-k-1))/(N-1); % final value of AmeetsA AmeetsB = (alpha1*p2+p1*alpha2)*(1-p1)^(k-1)*(1-p2)^(N-k-1) ... + (k>=2)*(k-1)*p1*(1-alpha1)*(1-p1)^(k-2)*p2*(1-p2)^(N-k-1) ... + (N-k>=2)*(N-k-1)*p1*(1-p1)^(k-1)*p2*(1-alpha2)*(1-p2)^(N-k-2); % probability a focal A meets a given B in a group of size 2 % larger groups will be accounted for in the following cycle for m=3:N-1 % here m stands for the size of the group where A and B meet l = min([k-1 m-2]); % max number of other As in the group a = max([0 m-1+k-N]); % need at least so many other As in the group groupa = 0; % they meet at focal A home or given B home for j=a:l % j is the number of other A's in the group groupa = groupa + nchoosek(k-1,j)*p1^j*(1-p1)^(k-1-j) ... *nchoosek(N-k-1,m-2-j)*p2^(m-2-j)*(1-p2)^(N-k-m+1+j); end groupa = groupa*(alpha1*p2+p1*alpha2); % focal A and given B l1 = min([k-2 m-3]); % max number of other As in the group a1 = max([0 m-2+k-N]); % need at least so many other As in the group groupc1 = 0; % they meet at another A's home and host stays for j=a1:l1 groupc1 = groupc1 + nchoosek(k-2,j)*p1^j*(1-p1)^(k-2-j) ... *nchoosek(N-k-1,m-3-j)*p2^(m-3-j)*(1-p2)^(N-k-m+2+j); end groupc1 = groupc1*alpha1; % the host stays l2 = min([k-2 m-2]); % max number of other As in the group a2 = max([0 m-1+k-N]); % need at least so many other As in the group groupc2 = 0; % they meet at another A's home and host leaves for j=a2:l2 groupc2 = groupc2 + nchoosek(k-2,j)*p1^j*(1-p1)^(k-2-j) ... *nchoosek(N-k-1,m-2-j)*p2^(m-2-j)*(1-p2)^(N-k-m+1+j); end groupc2 = groupc2*(1-alpha1); % the host leaves groupc = (k>=2)*(k-1)*p1*p2*(groupc1+groupc2); % need k>=2 for this % all possible groups of type c n1 = min([k-1 m-3]); % max number of other As in the group b1 = max([0 m-1+k-N]); % need at least so many other As in the group groupd1 = 0; % they meet at another B's home and host stays for j=b1:n1 groupd1 = groupd1 + nchoosek(k-1,j)*p1^j*(1-p1)^(k-1-j) ... *nchoosek(N-k-2,m-3-j)*p2^(m-3-j)*(1-p2)^(N-k-m+1+j); end groupd1 = groupd1*alpha2; % the host stays n2 = min([k-1 m-2]); % max number of other As in the group b2 = max([0 m+k-N]); % need at least so many other As in the group groupd2 = 0; % they meet at some B's home and host leaves for j=b2:n2 groupd2 = groupd2 + nchoosek(k-1,j)*p1^j*(1-p1)^(k-1-j) ... *nchoosek(N-k-2,m-2-j)*p2^(m-2-j)*(1-p2)^(N-k-m+j); end groupd2 = groupd2*(1-alpha2); % the host leaves groupd = (N-k>=2)*(N-k-1)*p1*p2*(groupd1+groupd2); % need N-k>=2 % all possible groups of type d AmeetsB = AmeetsB + (groupa+groupc+groupd)/(m-1); % adjust time end % we went through all the possible group sizes from 3 to N-1 % it now remains to account the group of size N AmeetsB = AmeetsB + (k*alpha1*p1^(k-1)*p2^(N-k) ... + (N-k)*p1^k*alpha2*p2^(N-k-1))/(N-1); % final value of AmeetsB BmeetsB = 2*alpha2*p2*(1-p1)^k*(1-p2)^(N-k-2) ... + k*p2^2*(1-alpha1)*(1-p1)^(k-1)*(1-p2)^(N-k-2) ... + (N-k>=3)*(N-k-2)*p2^2*(1-p1)^k*(1-alpha2)*(1-p2)^(N-k-3); % probability a focal B meets another given B in a group of size 2 % larger groups will be accounted for in the following cycle for m=3:N-1 % here m stands for the size of the group where two A's meet l = min([k m-2]); % max number of other As in the group a = max([0 m+k-N]); % need at least so many other As in the group groupa = 0; % they meet at focal A home or given A home for j=a:l % j is the number of other A's in the group groupa = groupa + nchoosek(k,j)*p1^j*(1-p1)^(k-j) ... *nchoosek(N-k-2,m-2-j)*p2^(m-2-j)*(1-p2)^(N-k-m+j); end groupa = groupa*alpha2*p2*2; % focal A and given A l1 = min([k-1 m-3]); % max number of other As in the group a1 = max([0 m-1+k-N]); % need at least so many other As in the group groupc1 = 0; % they meet at another A's home and host stays for j=a1:l1 groupc1 = groupc1 + nchoosek(k-1,j)*p1^j*(1-p1)^(k-1-j) ... *nchoosek(N-k-2,m-3-j)*p2^(m-3-j)*(1-p2)^(N-k-m+1+j); end groupc1 = groupc1*alpha1; % the host stays l2 = min([k-1 m-2]); % max number of other As in the group a2 = max([0 m+k-N]); % need at least so many other As in the group groupc2 = 0; % they meet at another A's home and host leaves for j=a2:l2 groupc2 = groupc2 + nchoosek(k-1,j)*p1^j*(1-p1)^(k-1-j) ... *nchoosek(N-k-2,m-2-j)*p2^(m-2-j)*(1-p2)^(N-k-m+j); end groupc2 = groupc2*(1-alpha1); % the host leaves groupc = k*p2^2*(groupc1+groupc2); % all possible groups of type c n1 = min([k m-3]); % max number of other As in the group b1 = max([0 m+k-N]); % need at least so many other As in the group groupd1 = 0; % they meet at some B's home and host stays for j=b1:n1 groupd1 = groupd1 + nchoosek(k,j)*p1^j*(1-p1)^(k-j) ... *nchoosek(N-k-3,m-3-j)*p2^(m-3-j)*(1-p2)^(N-k-m+j); end groupd1 = groupd1*alpha2; % the host stays n2 = min([k m-2]); % max number of other As in the group b2 = max([0 m+1+k-N]); % need at least so many other As in the group groupd2 = 0; % they meet at some B's home and host leaves for j=b2:n2 groupd2 = groupd2 + nchoosek(k,j)*p1^j*(1-p1)^(k-j) ... *nchoosek(N-k-3,m-2-j)*p2^(m-2-j)*(1-p2)^(N-k-m-1+j); end groupd2 = groupd2*(1-alpha2); % the host leaves groupd = (N-k>=3)*(N-k-2)*p2^2*(groupd1+groupd2); % all possible groups of type d BmeetsB = BmeetsB + (groupa+groupc+groupd)/(m-1); % adjust time end % we went through all the possible group sizes from 3 to N-1 % it now remains to account the group of size N BmeetsB = BmeetsB + (k*alpha1*p1^(k-1)*p2^(N-k) ... + (N-k)*p1^k*alpha2*p2^(N-k-1))/(N-1); % final value of BmeetsB for i=1:k % individual i is of type A for j=1:k % individual j is of type A RWM(i,j,k) = AmeetsA; % time i and j spend together end for j=k+1:N % individual j is of type B RWM(i,j,k) = AmeetsB; % time i and j spend together end RWM(i,i,k) = alone1; % correct probability individual i stays alone end for i=k+1:N % individual i is of type B for j=1:k % individual j is of type A RWM(i,j,k) = AmeetsB; % time i and j spend together end for j=k+1:N % individual j is of type B RWM(i,j,k) = BmeetsB; % time i and j spend together end RWM(i,i,k) = alone2; % correct probability individual i stays alone end end % initialize the matrix containing results of the simulation % column 1 = total number of generations the simulation took % column 2 = flag (0 = indecisive, 1 = A won, 2 = B won) Results = zeros(trials,2); % initialize the vector containing aggregate results of the simulation Statistics = zeros(1,13); Statistics(:,1) = N; Statistics(:,2) = alpha1; Statistics(:,3) = alpha2; Statistics(:,4) = beta1; Statistics(:,5) = beta2; Statistics(:,6) = S; Statistics(:,7) = c; Statistics(:,8) = v; Statistics(:,9) = lambda; Statistics(:,10) = T; Statistics(:,11) = A; % who the mutant is % column 12 = average number of generations to reach fixation % column 13 = fixation probability of the mutant %% simulation close all; rng('shuffle'); % reset the random generator seed for run=1:trials steps = 0; % counter for the length of each simulation % populate the graph initA = 1; % number of type A individuals initially (1 or N-1) IS(1:initA,type) = A*ones(initA,1); IS(initA+1:N,type) = B*ones(N-initA,1); % start the evolution process numberA = sum(IS(:,type)==A); % number of type A individuals % run the process until one type fixates while (numberA>0 && numberA stay_prob % individual moves to another random location new_loc = randi(Neighbors{k_pos,2}); % number of neighbor k_new_pos = Neighbors{k_pos,1}(new_loc); % new position % update individual statistics and node statistics IS(k,position) = k_new_pos; temp_Node_stats(k_pos,k_type) = ... temp_Node_stats(k_pos,k_type)-1; temp_Node_stats(k_new_pos,k_type) = ... temp_Node_stats(k_new_pos,k_type)+1; % pay movement cost IS(k,fitness) = IS(k,fitness) - lambda; end end % end of movement phase % record the new disposition of individuals Node_stats = temp_Node_stats; % play the public goods game for k=1:N k_type = IS(k,type); % type of individual k k_pos = IS(k,position); % position of individual k % add payoff from the new round of the game n_total = Node_stats(k_pos,A)+Node_stats(k_pos,B); % total number of individuals at the location if n_total==1 % if individual is alone IS(k,fitness) = IS(k,fitness) + 1 ... - c*(istrat(k_type)==cooperator); else % if invidual is not alone n_coop = Node_stats(k_pos,A)*(istrat(A)==cooperator) ... + Node_stats(k_pos,B)*(istrat(B)==cooperator); % number of cooperators at the location IS(k,fitness) = IS(k,fitness) + 1 + ... (n_coop-1*(istrat(k_type)==cooperator))/(n_total-1)*v ... - c*(istrat(k_type)==cooperator); end end % end of the public good game phase end % end of T movement steps % evolve the population via BDB process repr = [rand rand]; % start by determining who is going to reproduce parent = find((cumsum(IS(:,fitness))>=repr(1)*sum(IS(:,fitness))),1); % find out who is going to be replaced using RWM dead = find((cumsum(RWM(parent,:,numberA))>=repr(2)*... sum(RWM(parent,:,numberA))),1); % update the state of the environment if there is a change if IS(parent,type) ~= IS(dead,type) if IS(parent,type)==A IS(numberA+1,type)=A; % one more A individual numberA = numberA+1; else IS(numberA,type)=B; % one more B individual numberA = numberA-1; end end end % end of the simulation, population reached fixation Results(run,1) = steps; % time to reach fixation switch numberA case 0 Results(run,2) = B; % type B won case N Results(run,2) = A; % type A won end end % end of all trials % average and record the results of the simulation Statistics(1,12) = mean(Results(:,1)); Awon = sum(Results(:,2)==A); Bwon = sum(Results(:,2)==B); Statistics(1,13) = Awon/(Awon+Bwon); % fixation probability of A dlmwrite('comp_pg_coopE_fine_N50.csv', Statistics, '-append'); end %