%This is the code used to simulate the dynamics of range expansion in "The
%effect of the recombination rate between adaptive loci on the capacity of
%a population to expand its range by M. Eriksson and M. Rafajlovic
rng('shuffle')
global M
global rm
global loci
global alleffect
global Vs
global sig
global self

%Parameters
rm = 1.025;                                         %Maximal growth rate
Vs = 1/2;                                           %Stabilizing selection
sig = 1/sqrt(2);                                    %Dispersal
mu = 1e-6;                                          %Mutation rate
recrate = 0.0001;                                   %Recombination rate
self = 1;                                           %Probability to accept selfing
b = 0.1;                                            %Slope of linear phenotypic optimum
alleffect = 0.25;                                   %Allelic effects sizes
M = 179;                                            %Number of habitat patches
[theta,ccenters,dtheta,loci] = DiploidHabitat(b);   %Phenotypic optimum, cline centers and number of loci

%Repetition parameters
g = 10;                                            %Number of generations
snap = 1;                                           %Number of generations between each saved snapshot
TL = snap*ones(1,g/snap);                           %Generate the measurement points

%Create initial conditions for the simulation
init_range = 1;                                     %Type of starting distribution 0 = uniform, 1 = localized
%N.B. only init_range = 1 was used in the submitted study!
K1 = 150;                                           %Maximum carrying capacity
InitVar = 1;                                        %0 = minimum, 1 = maximum
[input_pop,K,N,p] = InitConds(ccenters,K1,init_range,InitVar);

%Run simulation
[output_pop,N,Npm,Ntot,Vg,Pself,cline,selfpenalty,z,nohap,meanfitness] = DiploidMain(mu,recrate,K,N,theta,input_pop,TL);

%Calculate the critical gradient
Beff = dtheta*sig./((rm-dtheta*sig/(2*sqrt(Vs)))*sqrt(2*Vs));
Crit = 0.15*sig*sqrt((2*alleffect)^2/(2*Vs));
Nexp = K.*(1-sig*dtheta/(2*sqrt(Vs)*rm));

%Critical gradient is: Beff>Crit*Nexp