%% info
% this is a simulation for public goods game on
% star graphs
function starpg_bc50(par1,par2,par3,par4)
% par1 is alpha1
% par2 is alpha2
% par3 is lambda
% par4 is v/c ratio
%% 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 = c*par4; % 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 = zeros(N); % edge matrix for star graph
E(:,1) = ones(N,1);
E(1,:) = ones(1,N);
E(1,1) = 0;
% 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,2); % replacement weight matrix placeholder
% third index corresponds to the number of A individuals in the population
% last index corresponds to the type of the individual in the center
% we assume that the center of the star is the first node
% followed by A individuals
% followed by B individuals
% first, we look at the case when A is in the center
alphac = alpha1; % A is in the center
pnc = (N-2+alphac)/(N-1); % probability center doesn't move to a given leaf
for k=1:N-1
% k denotes the number of A individuals in the population
center_alone = alpha1^k*alpha2^(N-k) ...
+ ((k-1)*(1-alphac)*(1-alpha1) ...
+ (N-k)*(1-alphac)*(1-alpha2))/(N-1);
leafA_alone = alpha1*pnc ...
+ (1-alpha1)*(1-alphac)*alpha1^(k-2)*alpha2^(N-k);
leafB_alone = alpha2*pnc ...
+ (1-alpha2)*(1-alphac)*alpha1^(k-1)*alpha2^(N-k-1);
% now we will compute center --> leaf A weight
center_meets_leafA = (1-alphac)/(N-1)*alpha1 ...
+ (1-alpha1)*alpha1^(k-1)*alpha2^(N-k); % group of size 2
for m=3:N-1
% group of size m
l = min([k-2 m-2]); % max number of As in the group
a = max([0 m-2+k-N]); % need at least this many As
tempCA = 0; % this will accumulate the formula for m-groups
for j=a:l
tempCA = tempCA + nchoosek(k-2,j)*(1-alpha1)^j*alpha1^(k-2-j) ...
*nchoosek(N-k,m-2-j)*(1-alpha2)^(m-2-j)*alpha2^(N-k-m+2+j);
end
center_meets_leafA = center_meets_leafA ...
+ tempCA*alphac*(1-alpha1)/(m-1); % accounted for m-group
end
center_meets_leafA = center_meets_leafA ...
+ alphac*(1-alpha1)^(k-1)*(1-alpha2)^(N-k)/(N-1); % N-group
% this concludes the computation of center --> leaf A weight
% now we will compute center --> leaf B weight
center_meets_leafB = (1-alphac)/(N-1)*alpha2 ...
+ (1-alpha2)*alpha1^k*alpha2^(N-k-1); % group of size 2
for m=3:N-1
% group of size m
l = min([k-1 m-2]); % max number of As in the group
a = max([0 m-1+k-N]); % need at least this many As
tempCB = 0; % this will accumulate the formula for m-groups
for j=a:l
tempCB = tempCB + nchoosek(k-1,j)*(1-alpha1)^j*alpha1^(k-1-j) ...
*nchoosek(N-k-1,m-2-j)*(1-alpha2)^(m-2-j)*alpha2^(N-k-m+1+j);
end
center_meets_leafB = center_meets_leafB ...
+ tempCB*alphac*(1-alpha2)/(m-1); % accounted for m-group
end
center_meets_leafB = center_meets_leafB ...
+ alphac*(1-alpha1)^(k-1)*(1-alpha2)^(N-k)/(N-1); % N-group
% this concludes the computation of center --> leaf B weight
% now we will compute leaf A --> leaf A weight
leafA_meets_leafA = (1-alphac)*(1-alpha1)^2*alpha1^(k-3)*alpha2^(N-k);
% group of size 2
for m=3:N-1
% group of size m
l1 = min([k-3 m-3]); % max number of As in the group
a1 = max([0 m-3+k-N]); % need at least this many As
tempAA1 = 0; % this will accumulate the formula for m-groups when
% the center stays home
for j=a1:l1
tempAA1 = tempAA1 + nchoosek(k-3,j)*(1-alpha1)^j*alpha1^(k-3-j) ...
*nchoosek(N-k,m-3-j)*(1-alpha2)^(m-3-j)*alpha2^(N-k-m+3+j);
end
l2 = min([k-3 m-2]); % max number of As in the group
a2 = max([0 m-2+k-N]); % need at least this many As
tempAA2 = 0; % this will accumulate the formula for m-groups when
% the center leaves
for j=a2:l2
tempAA2 = tempAA2 + nchoosek(k-3,j)*(1-alpha1)^j*alpha1^(k-3-j) ...
*nchoosek(N-k,m-2-j)*(1-alpha2)^(m-2-j)*alpha2^(N-k-m+2+j);
end
tempAA = (alphac*tempAA1+(1-alphac)*tempAA2)*(1-alpha1)^2;
leafA_meets_leafA = leafA_meets_leafA + tempAA/(m-1);
% accounted for m-group
end
leafA_meets_leafA = leafA_meets_leafA ...
+ alphac*(1-alpha1)^(k-1)*(1-alpha2)^(N-k)/(N-1); % N-group
% this concludes the computation of leaf A --> leaf A weight
% now we will compute leaf A --> leaf B weight
leafA_meets_leafB = (1-alphac)*(1-alpha1)*(1-alpha2) ...
*alpha1^(k-2)*alpha2^(N-k-1); % group of size 2
for m=3:N-1
% group of size m
l1 = min([k-2 m-3]); % max number of As in the group
a1 = max([0 m-2+k-N]); % need at least this many As
tempAB1 = 0; % this will accumulate the formula for m-groups when
% the center stays home
for j=a1:l1
tempAB1 = tempAB1 + nchoosek(k-2,j)*(1-alpha1)^j*alpha1^(k-2-j) ...
*nchoosek(N-k-1,m-3-j)*(1-alpha2)^(m-3-j)*alpha2^(N-k-m+2+j);
end
l2 = min([k-2 m-2]); % max number of As in the group
a2 = max([0 m-1+k-N]); % need at least this many As
tempAB2 = 0; % this will accumulate the formula for m-groups when
% the center leaves
for j=a2:l2
tempAB2 = tempAB2 + nchoosek(k-2,j)*(1-alpha1)^j*alpha1^(k-2-j) ...
*nchoosek(N-k-1,m-2-j)*(1-alpha2)^(m-2-j)*alpha2^(N-k-m+1+j);
end
tempAB = (alphac*tempAB1+(1-alphac)*tempAB2)*(1-alpha1)*(1-alpha2);
leafA_meets_leafB = leafA_meets_leafB + tempAB/(m-1);
% accounted for m-group
end
leafA_meets_leafB = leafA_meets_leafB ...
+ alphac*(1-alpha1)^(k-1)*(1-alpha2)^(N-k)/(N-1); % N-group
% this concludes the computation of leaf A --> leaf B weight
% now we will compute leaf B --> leaf B weight
leafB_meets_leafB = (1-alphac)*(1-alpha2)^2*alpha1^(k-1)*alpha2^(N-k-2);
% group of size 2
for m=3:N-1
% group of size m
l1 = min([k-1 m-3]); % max number of As in the group
a1 = max([0 m-1+k-N]); % need at least this many As
tempBB1 = 0; % this will accumulate the formula for m-groups when
% the center stays home
for j=a1:l1
tempBB1 = tempBB1 + nchoosek(k-1,j)*(1-alpha1)^j*alpha1^(k-1-j) ...
*nchoosek(N-k-2,m-3-j)*(1-alpha2)^(m-3-j)*alpha2^(N-k-m+1+j);
end
l2 = min([k-1 m-2]); % max number of As in the group
a2 = max([0 m+k-N]); % need at least this many As
tempBB2 = 0; % this will accumulate the formula for m-groups when
% the center leaves
for j=a2:l2
tempBB2 = tempBB2 + nchoosek(k-1,j)*(1-alpha1)^j*alpha1^(k-1-j) ...
*nchoosek(N-k-2,m-2-j)*(1-alpha2)^(m-2-j)*alpha2^(N-k-m+j);
end
tempBB = (alphac*tempBB1+(1-alphac)*tempBB2)*(1-alpha2)^2;
leafB_meets_leafB = leafB_meets_leafB + tempBB/(m-1);
% accounted for m-group
end
leafB_meets_leafB = leafB_meets_leafB ...
+ alphac*(1-alpha1)^(k-1)*(1-alpha2)^(N-k)/(N-1); % N-group
% this concludes the computation of leaf B --> leaf B weight
% special formulas apply to the center
RWM(1,1,k,1) = center_alone;
for j=2:k
% leaf j is occupied by A
RWM(1,j,k,1) = center_meets_leafA;
RWM(j,1,k,1) = center_meets_leafA;
end
for j=k+1:N
% leaf j is occupied by B
RWM(1,j,k,1) = center_meets_leafB;
RWM(j,1,k,1) = center_meets_leafB;
end
for i=2:k
% leaf i is occupied by A
for j=2:k
% leaf j is occupied by A
RWM(i,j,k,1) = leafA_meets_leafA;
end
for j=k+1:N
% leaf j is occupied by B
RWM(i,j,k,1) = leafA_meets_leafB;
end
RWM(i,i,k,1) = leafA_alone; % correct formula for itself
end
for i=k+1:N
% leaf i is occupied by B
for j=2:k
% leaf j is occupied by A
RWM(i,j,k,1) = leafA_meets_leafB;
end
for j=k+1:N
% leaf j is occupied by B
RWM(i,j,k,1) = leafB_meets_leafB;
end
RWM(i,i,k,1) = leafB_alone; % correct formula for itself
end
end % went over all possible numbers of A individuals given A center
% second, we look at the case when B is in the center
alphac = alpha2; % B is in the center
pnc = (N-2+alphac)/(N-1); % probability center doesn't move to a given leaf
for k=1:N-1
% k denotes the number of A individuals in the population
center_alone = alpha1^k*alpha2^(N-k) ...
+ (k*(1-alphac)*(1-alpha1) ...
+ (N-k-1)*(1-alphac)*(1-alpha2))/(N-1);
leafA_alone = alpha1*pnc ...
+ (1-alpha1)*(1-alphac)*alpha1^(k-1)*alpha2^(N-k-1);
leafB_alone = alpha2*pnc ...
+ (1-alpha2)*(1-alphac)*alpha1^k*alpha2^(N-k-2);
% now we will compute center --> leaf A weight
center_meets_leafA = (1-alphac)/(N-1)*alpha1 ...
+ (1-alpha1)*alpha1^(k-1)*alpha2^(N-k); % group of size 2
for m=3:N-1
% group of size m
l = min([k-1 m-2]); % max number of As in the group
a = max([0 m-1+k-N]); % need at least this many As
tempCA = 0; % this will accumulate the formula for m-group
for j=a:l
tempCA = tempCA + nchoosek(k-1,j)*(1-alpha1)^j*alpha1^(k-1-j) ...
*nchoosek(N-k-1,m-2-j)*(1-alpha2)^(m-2-j)*alpha2^(N-k-m+1+j);
end
center_meets_leafA = center_meets_leafA ...
+ tempCA*alphac*(1-alpha1)/(m-1); % accounted for m-group
end
center_meets_leafA = center_meets_leafA ...
+ alphac*(1-alpha1)^k*(1-alpha2)^(N-k-1)/(N-1); % N-group
% this concludes the computation of center --> leaf A weight
% now we will compute center --> leaf B weight
center_meets_leafB = (1-alphac)/(N-1)*alpha2 ...
+ (1-alpha2)*alpha1^k*alpha2^(N-k-1); % group of size 2
for m=3:N-1
% group of size m
l = min([k m-2]); % max number of As in the group
a = max([0 m+k-N]); % need at least this many As
tempCB = 0; % this will accumulate the formula for m-groups
for j=a:l
tempCB = tempCB + nchoosek(k,j)*(1-alpha1)^j*alpha1^(k-j) ...
*nchoosek(N-k-2,m-2-j)*(1-alpha2)^(m-2-j)*alpha2^(N-k-m+j);
end
center_meets_leafB = center_meets_leafB ...
+ tempCB*alphac*(1-alpha2)/(m-1); % accounted for m-group
end
center_meets_leafB = center_meets_leafB ...
+ alphac*(1-alpha1)^k*(1-alpha2)^(N-k-1)/(N-1); % N-group
% this concludes the computation of center --> leaf B weight
% now we will compute leaf A --> leaf A weight
leafA_meets_leafA = (1-alphac)*(1-alpha1)^2*alpha1^(k-2) ...
*alpha2^(N-k-1); % group of size 2
for m=3:N-1
% group of size m
l1 = min([k-2 m-3]); % max number of As in the group
a1 = max([0 m-2+k-N]); % need at least this many As
tempAA1 = 0; % this will accumulate the formula for m-groups when
% the center stays home
for j=a1:l1
tempAA1 = tempAA1 + nchoosek(k-2,j)*(1-alpha1)^j*alpha1^(k-2-j) ...
*nchoosek(N-k-1,m-3-j)*(1-alpha2)^(m-3-j)*alpha2^(N-k-m+2+j);
end
l2 = min([k-2 m-2]); % max number of As in the group
a2 = max([0 m-1+k-N]); % need at least this many As
tempAA2 = 0; % this will accumulate the formula for m-groups when
% the center leaves
for j=a2:l2
tempAA2 = tempAA2 + nchoosek(k-2,j)*(1-alpha1)^j*alpha1^(k-2-j) ...
*nchoosek(N-k-1,m-2-j)*(1-alpha2)^(m-2-j)*alpha2^(N-k-m+1+j);
end
tempAA = (alphac*tempAA1+(1-alphac)*tempAA2)*(1-alpha1)^2;
leafA_meets_leafA = leafA_meets_leafA + tempAA/(m-1);
% accounted for m-group
end
leafA_meets_leafA = leafA_meets_leafA ...
+ alphac*(1-alpha1)^k*(1-alpha2)^(N-k-1)/(N-1); % N-group
% this concludes the computation of leaf A --> leaf A weight
% now we will compute leaf A --> leaf B weight
leafA_meets_leafB = (1-alphac)*(1-alpha1)*(1-alpha2) ...
*alpha1^(k-1)*alpha2^(N-k-2); % group of size 2
for m=3:N-1
% group of size m
l1 = min([k-1 m-3]); % max number of As in the group
a1 = max([0 m-1+k-N]); % need at least this many As
tempAB1 = 0; % this will accumulate the formula for m-groups when
% the center stays home
for j=a1:l1
tempAB1 = tempAB1 + nchoosek(k-1,j)*(1-alpha1)^j*alpha1^(k-1-j) ...
*nchoosek(N-k-2,m-3-j)*(1-alpha2)^(m-3-j)*alpha2^(N-k-m+1+j);
end
l2 = min([k-1 m-2]); % max number of As in the group
a2 = max([0 m+k-N]); % need at least this many As
tempAB2 = 0; % this will accumulate the formula for m-groups when
% the center leaves
for j=a2:l2
tempAB2 = tempAB2 + nchoosek(k-1,j)*(1-alpha1)^j*alpha1^(k-1-j) ...
*nchoosek(N-k-2,m-2-j)*(1-alpha2)^(m-2-j)*alpha2^(N-k-m+j);
end
tempAB = (alphac*tempAB1+(1-alphac)*tempAB2)*(1-alpha1)*(1-alpha2);
leafA_meets_leafB = leafA_meets_leafB + tempAB/(m-1);
% accounted for m-group
end
leafA_meets_leafB = leafA_meets_leafB ...
+ alphac*(1-alpha1)^k*(1-alpha2)^(N-k-1)/(N-1); % N-group
% this concludes the computation of leaf A --> leaf B weight
% now we will compute leaf B --> leaf B weight
leafB_meets_leafB = (1-alphac)*(1-alpha2)^2*alpha1^k*alpha2^(N-k-3);
% group of size 2
for m=3:N-1
% group of size m
l1 = min([k m-3]); % max number of As in the group
a1 = max([0 m+k-N]); % need at least this many As
tempBB1 = 0; % this will accumulate the formula for m-groups when
% the center stays home
for j=a1:l1
tempBB1 = tempBB1 + nchoosek(k,j)*(1-alpha1)^j*alpha1^(k-j) ...
*nchoosek(N-k-3,m-3-j)*(1-alpha2)^(m-3-j)*alpha2^(N-k-m+j);
end
l2 = min([k m-2]); % max number of As in the group
a2 = max([0 m+1+k-N]); % need at least this many As
tempBB2 = 0; % this will accumulate the formula for m-groups when
% the center leaves
for j=a2:l2
tempBB2 = tempBB2 + nchoosek(k,j)*(1-alpha1)^j*alpha1^(k-j) ...
*nchoosek(N-k-3,m-2-j)*(1-alpha2)^(m-2-j)*alpha2^(N-k-m-1+j);
end
tempBB = (alphac*tempBB1+(1-alphac)*tempBB2)*(1-alpha2)^2;
leafB_meets_leafB = leafB_meets_leafB + tempBB/(m-1);
% accounted for m-group
end
leafB_meets_leafB = leafB_meets_leafB ...
+ alphac*(1-alpha1)^k*(1-alpha2)^(N-k-1)/(N-1); % N-group
% this concludes the computation of leaf B --> leaf B weight
% special formulas apply to the center
RWM(1,1,k,2) = center_alone;
for j=2:k+1
% leaf j is occupied by A
RWM(1,j,k,2) = center_meets_leafA;
RWM(j,1,k,2) = center_meets_leafA;
end
for j=k+2:N
% leaf j is occupied by B
RWM(1,j,k,2) = center_meets_leafB;
RWM(j,1,k,2) = center_meets_leafB;
end
for i=2:k+1
% leaf i is occupied by A
for j=2:k+1
% leaf j is occupied by A
RWM(i,j,k,2) = leafA_meets_leafA;
end
for j=k+2:N
% leaf j is occupied by B
RWM(i,j,k,2) = leafA_meets_leafB;
end
RWM(i,i,k,2) = leafA_alone; % correct formula for itself
end
for i=k+2:N
% leaf i is occupied by B
for j=2:k+1
% leaf j is occupied by A
RWM(i,j,k,2) = leafA_meets_leafB;
end
for j=k+2:N
% leaf j is occupied by B
RWM(i,j,k,2) = leafB_meets_leafB;
end
RWM(i,i,k,2) = leafB_alone;
end
end % went over all possible numbers of A individuals given B center
% 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 with one type A mutant in a random position
mutant_position = randi(N);
if mutant_position==1
% mutant goes to the center
IS(1,type)=A;
IS(2:N,type)=B*ones(N-1,1);
else
% mutant goes to a leaf
IS(1,type)=B;
IS(2,type)=A;
IS(3:N,type)=B*ones(N-2,1);
end
% 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,IS(1,type)))>=repr(2)*...
sum(RWM(parent,:,numberA,IS(1,type)))),1);
% update the state of the environment if there is a change
if IS(parent,type) ~= IS(dead,type)
if dead==1
% center gets replaced
IS(1,type)=IS(parent,type);
if IS(parent,type)==A
numberA = numberA+1; % one more A individual
else
numberA = numberA-1; % one more B individual
end
else
% leaf gets replaced
if IS(parent,type)==A
% one more A individual
% its position depends on who is in the center
% recall A=1 and B=2
IS(numberA+IS(1,type),type)=A;
numberA = numberA+1;
else
% one more B individual
IS(numberA-1+IS(1,type),type)=B;
numberA = numberA-1;
end
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('star_pg_bc_coopE_fine_N50.csv', Statistics, '-append');
end
%