Setting up

The data are in the file:  “DispersalData.csv”, found in the supplementary material. 
These should be put in the same directory, and the directory name set below, under “Import data”.
Then, “Evaluate Initialization” to set up all data and definitions.

Simulating a uniform spatial pedigree for 1000 generations with gaussian dispersal Kernel

Initializations for simulation

Large Area

nuni is the number of plants simulated, Rx and Ry define the border XY coordinates of the spatial area of the simulations, tm is the generation time being simulated,
startuni creates the initial locations, and inigenmeetuni creates the genetic starting condition where on one side of the midpoint we find only homozygots with {0, 0 } and on the other with {1, 1} resembling a secondary contact scenario

nuni=15000;
Rx={418800,434700};
Ry={4684700,4687800};

tm= 1000;

s = 0.01

start uni creates the starting locations of the first generation

startuni = makeN[nuni,Rx,Ry]

as a starting condition the two genotypes meet at the midpoint, re-creating a secondary contact scenario

midpoint = 426750;

inigenmeetuni = If[#<=midpoint,{0,0},{1,1}]&/@Flatten[startuni][[1;;-1;;2]] ;

Small Area

for a smaller area and N, choose  nsmall, and the small area and starting condition

nsmall = 5262

Rysmall = {4684700,4687025 };

Rxsmall ={418800,426750};

tm = 1000;
s = 0.15

startunismall = makeN[nsmall,Rxsmall,Rysmall];

midsmall = 422775;

inigenmeetsmall= If[#<=midsmall,{0,0},{1,1}]&/@Flatten[startunismall][[1;;-1;;2]] ;

Generate a spatial pedigree

Given the starting locations,  gg generates a spatial pedigree storing the genotypes, locations and parents of every individual for all generations .  Set a different value for  tm and s under "Initializations for simulation" to generation time and selection strength

gg= makeGenealogyRealuni[startunismall,tm,Rxsmall,Rysmall,inigenmeetsmall,s];//Timing

Export["Unifodmdist_normaldisp_15000inds_s001",gg[[1;;-1;;10]]]

Results

s = 0.1, n = 8000, uniform distribution, Gaussian dispersal

Importing

Rep1= ToExpression[Import["/Users/mjulseth/Documents/Barton/small/Rep1_2/smallarea_8000_01_allgen_lessdisp.csv"]][[-1]];

Rep2 = ToExpression[Import["/Users/mjulseth/Documents/Barton/small/Rep1_2/smallarea_8000_01_allgen_lessdisp_2.csv"]][[-1]];

Rep3= ToExpression[Import["/Users/mjulseth/Documents/Barton/small/Rep1_2/smallarea_8000_01_allgen_lessdisp_3.csv"]][[-1]];

Rep4= ToExpression[Import["/Users/mjulseth/Documents/Barton/small/Rep1_2/smallarea_8000_01_allgen_lessdisp_4.csv"]][[-1]];

Rep5= ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\Rep1_2\\smallarea_8000_01_allgen_lessdisp_5.csv"]][[-1]];

basic sigmoid

sig1= makeMetSNP[Rep1,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig2= makeMetSNP[Rep2,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig3= makeMetSNP[Rep3,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig4= makeMetSNP[Rep4,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig5= makeMetSNP[Rep5,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl1=mleMetSig[sig1,2500]

cl2=mleMetSig[sig2,2500]

cl3=mleMetSig[sig3,2500]

cl4=mleMetSig[sig4,2500]

cl5=mleMetSig[sig5,2500]

wi018 = {842.909,788.244,975.009,926.021,944.489}

center018 = {12283.5, 12662.4,12670.9,12461.2,12947.2}

avwi018 = Mean[wi018]

avcen = Mean[center018]

sig1= makeMetSNP[Rep1,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig2= makeMetSNP[Rep2,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig3= makeMetSNP[Rep3,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig4= makeMetSNP[Rep4,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig5= makeMetSNP[Rep5,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl1=mleMetSig[sig1,2500]

cl2=mleMetSig[sig2,2500]

cl3=mleMetSig[sig3,2500]

cl4=mleMetSig[sig4,2500]

cl5=mleMetSig[sig5,2500]

Basic Sigmoid Cline Plot

Cline1sigmoid =Plot[pCline[{842.909,12283.5}][x],{x,11000,14500}, PlotStyle->{Blue},PlotLabel->"Replicate 1 Sigmoid"] ;

Cline2sigmoid =Plot[pCline[{788.244,12662.4}][x],{x,11000,14500}, PlotStyle->{Green},PlotLabel->"Replicate 2 Sigmoid"] ;

Cline3sigmoid =Plot[pCline[{975.009,12670.9}][x],{x,11000,14500}, PlotStyle->{Orange},PlotLabel->"Replicate 3 Sigmoid"] ;

Cline4sigmoid =Plot[pCline[{926.021,12461.2}][x],{x,11000,14500}, PlotStyle->{Purple},PlotLabel->"Replicate 4 Sigmoid"] ;

Cline5sigmoid =Plot[pCline[{944.498,12947.2}][x],{x,11000,14500}, PlotStyle->{Pink},PlotLabel->"Replicate 5 Sigmoid"] ;

Average basic sigmoid cline plot

wi018sig = {842.909,788.244,975.009,926.021,944.489};

center018sig = {12283.5, 12662.4,12670.9,12461.2,12947.2};

avwi018 = Mean[wi018sig]

avcen018sig = Mean[center018sig]

StandardDeviation[center018sig]

StandardDeviation[wi018sig]

AveragClines01n8000 = Plot[pCline[{895.2,12605}][x],{x,11000,14500}, PlotStyle->{Black, Dashed},PlotLabel->"Average sigmoid cline s = 0.1. N = 8000"]

2 D Plot

tgg=Transpose[Rep1];

g00=Cases[tgg,{{0,0},_,_}][[All,3]];g01=Cases[tgg,{{0,1}|{1,0},_,_}][[All,3]];g11=Cases[tgg,{{1,1},_,_}][[All,3]];
Show[ListPlot[g00,PlotStyle->Black],ListPlot[g01,PlotStyle->Red],ListPlot[g11,PlotStyle->Blue],PlotRange->All,AspectRatio->0.19]

stepped

stepped1 =  makeMetasym[Rep1,Log[{2000,2000,0.05,0.05,2000,15000,0.1}],ConstantArray[0.1,7],350,1.01,3000];

stepped2 =  makeMetasym[Rep2,Log[{2000,2000,0.05,0.05,2000,15000,0.1}],ConstantArray[0.1,7],350,1.01,3000];

clstep1=mleMetasym[stepped1,2000]

clstep2=mleMetasym[stepped2,2000]

step1 =Plot[pbCline[{1.48026043062880292941`8.01031001176798*^38200062,0.`,2.068706527202383`*^272,0.`,830.1589647151136`,12285.994915741094`}][x],{x,8500,16450},PlotStyle->{Purple,Dashed}, PlotLegends->{"Stepped"}, PlotLabel->"Replicate 1 Stepped"];

step2 =Plot[pbCline[{4.6967600168668934170233570574`4.842499065392325*^56217951447,1.211810941967958032490842689855196`4.564599385030056*^106604293539,7.234191745262121`*^254,0.`,782.258498715769`,12652.939319710065`}][x],{x,11000,15000},PlotStyle->{Orange,Dashed}, PlotLegends->{"Stepped"}, PlotLabel->"Replicate 2 Stepped"];

Show[ Clinefreq1sigmoid,Clinefreq2sigmoid,stepfreq1,stepfreq2]

Frequency Dependent Selection

freq1 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\FreqDep\\smallarea_8000_01_allgen_FreqDep_1.csv"]][[-1]];

freq2 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\FreqDep\\smallarea_8000_01_allgen_FreqDep_2.csv"]][[-1]];

freq3 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\FreqDep\\smallarea_8000_01_allgen_FreqDep_1.csv"]][[-1]];

freq4 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\FreqDep\\smallarea_8000_01_allgen_FreqDep_2.csv"]][[-1]];

sigfreq1 = makeMetSNP[freq1,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sigfreq2 = makeMetSNP[freq2,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sigfreq3 = makeMetSNP[freq3,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

clfreq1=mleMetSig[sigfreq1,2500]

Clinefreq1sigmoid =Plot[pCline[{956.2891604726117`,12188.831466673682`}][x],{x,11000,14500}, PlotStyle->{Blue},PlotLabel->" Freq Replicate 1 Sigmoid"] ;

clfreq2=mleMetSig[sigfreq2,2500]

Clinefreq2sigmoid =Plot[pCline[{1196.54,12960.7}][x],{x,11000,14500}, PlotStyle->{Green},PlotLabel->" Freq Replicate 1 Sigmoid"]

clfreq3=mleMetSig[freq3,2500]

centerfreq = {12188.8,12960.7,12973.4}

avecenter = Mean[centerfreq]

wifreq = {956.289,1197.54,1159.11}

TTest[{wifreq,wi018sig}]

AveragClinesFreq= Plot[pCline[{1104.31,12707.6}][x],{x,11000,14500}, PlotStyle->{Red, Dashed},PlotLabel->"Average sigmoid cline s = 0.1. N = 8000"]

Average clines of frequency - dependent selection and selection against heterozygotes Plot

Show[AveragClines01n8000,AveragClinesFreq]

Width vs . selection, uniform distribution, Gaussian dispersal

s = 0.15, N = 8000

rep1158 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_8000_015_lastgen_new_1.csv"]];

rep2158= ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_8000_015_lastgen_new_2.csv"]];

rep3158 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_8000_015_lastgen_new_3.csv"]];

rep4158=ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_8000_015_allgen_lessdisp_new_3.csv"]][[-1]];

rep5158=ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_015_allgen_lessdisp_new_1.csv"]][[-1]];

sigmoidrep1=makeMetSNP[rep1158,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl0158=mleMetSig[sigmoidrep1,2500]

sigmoidrep2=makeMetSNP[rep2158,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl20158=mleMetSig[sigmoidrep2,2500]

sigmoidrep3=makeMetSNP[rep3158,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl30158=mleMetSig[sigmoidrep3,2500]

sigmoidrep4=makeMetSNP[rep4158,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl40158=mleMetSig[sigmoidrep4,2500]

sigmoidrep5=makeMetSNP[rep5158,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl50158=mleMetSig[sigmoidrep3,2500]

wi0158 = {725.35,860.46,618.329,650.716, 618.329}

av0158 = Mean[wi0158]

errorwi0158 = StandardDeviation[wi0158]/Sqrt[Length[wi0158]]

s = 0.2, N = 8000

s02n800rep1 = ToExpression[Import["/Users/mjulseth/Documents/Barton/small/smallarea_8000_02_allgen_lessdisp_1.csv"]];
s02n800rep2 = ToExpression[Import["/Users/mjulseth/Documents/Barton/small/smallarea_8000_02_allgen_lessdisp_2.csv"]];
rep1s028n = s02n800rep1[[-1]];
rep2s028n = s02n800rep2[[-1]];

sigmoid1028N=makeMetSNP[rep1s028n,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl1028n=mleMetSig[sigmoid1028N,2500]

sigmoid2018N=makeMetSNP[rep2s028n,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl2028n=mleMetSig[sigmoid2018N,2500]

s02rep3 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_8000_02_lastgen_new_3.csv"]];

s02rep4 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_8000_02_lastgen_new_2.csv"]];

s02rep5 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_8000_02_lastgen_new_1.csv"]];

sig0281= makeMetSNP[s02rep3,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig0282= makeMetSNP[s02rep4,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig0283= makeMetSNP[s02rep5,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

wi028 = {547.882,692,78, 694.608,663.687,804.507}

av028 =Mean[wi028]

errorwi028 = StandardDeviation[wi028]/Sqrt[Length[wi028]]

s = 0.1 , N = 5262

s01n5rep1 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5626_01_allgen_lessdisp_1.csv"]];

s01n5rep2 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5626_01_allgen_lessdisp_2.csv"]];

rep1015 = s01n5rep1[[-1]];

rep2015 = s01n5rep2[[-1]];

sigmoid1015N=makeMetSNP[rep1015,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl1015n=mleMetSig[sigmoid1015N,2500]

sigmoid1015N=makeMetSNP[rep2015,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl1015n=mleMetSig[sigmoid1015N,2500]

s015rep1 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_01_lastgen_new_1.csv"]];

s015rep2 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_01_lastgen_new_2.csv"]];

s015rep3 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_01_lastgen_new_3.csv"]];

sig015= makeMetSNP[s015rep1,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig0152= makeMetSNP[s015rep2,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig0153= makeMetSNP[s015rep3,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl0151=mleMetSig[sig015,2500]

cl0152=mleMetSig[sig0152,2500]

cl0153=mleMetSig[sig0153,2500]

wi015 = {887.60,961.447,1133.99,1148.59}

av015 = Mean[wi015]

err015 = StandardDeviation[wi015]/Sqrt[Length[wi015]]

s = 0.15, N = 5262

s015rep1 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_8000_015_allgen_lessdisp_new_3.csv"]]

rep2s015n5  =ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_015_lastgen_1.csv"]];

rep1s015n5= s015n5rep1[[-1]];

rep2s015n5 =s015n5rep2[[-1]];

sigmoid10155N=makeMetSNP[rep1s015n5,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl10155n=mleMetSig[sigmoid10155N,2500]

sigmoid20155N=makeMetSNP[rep2s015n5,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl20155n=mleMetSig[sigmoid20155N,2500]

s0155rep1 =ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_015_lastgen_new_1.csv"]];

s0155rep2 =ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_015_lastgen_new_2.csv"]];

s0155rep3 =ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_015_lastgen_new_3.csv"]];

sig01551= makeMetSNP[s0155rep1,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig01552= makeMetSNP[s0155rep2,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig01553= makeMetSNP[s0155rep3,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl01551=mleMetSig[sig01551,2500]

cl01552=mleMetSig[sig01552,2500]

cl01553=mleMetSig[sig01553,2500]

wi0155 = {719.02,908.21, 629.045,759.624,949.002}

avwi0155 = Mean[wi0155]

err0155 = StandardDeviation[wi0155]/Sqrt[Length[wi0155]]

s = 0.2, N = 5262

s02n5rep1 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_02_allgen_lessdisp_1.csv"]];

s02n5rep2 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_02_allgen_lessdisp_2.csv"]];

rep1s02n5 = s02n5rep1[[-1]];

rep2s02n5 = s02n5rep2[[-1]];

sigmoid1025N=makeMetSNP[rep1s02n5,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl1025n=mleMetSig[sigmoid1025N,2500]

sigmoid2025N=makeMetSNP[rep2s02n5,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl2025n=mleMetSig[sigmoid2025N,2500]

s025rep1 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_02_lastgen_new_1.csv"]];

s025rep2 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_02_lastgen_new_2.csv"]];

s025rep3 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\small\\smallarea_5262_02_lastgen_new_3.csv"]];

sig0251= makeMetSNP[s025rep1,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig0252= makeMetSNP[s025rep2,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

sig0253= makeMetSNP[s025rep3,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl251=mleMetSig[sig0251,2500]

cl252=mleMetSig[sig0252,2500]

cl253=mleMetSig[sig0253,2500]

wi025 = {712.29,810.25,776.612,671.978}

avwi025 = Mean[wi025]

err025 = StandardDeviation[wi025]/Sqrt[Length[wi025]]

Plot

sel = {0.1,0.15,0.2}

sigma1d = Sqrt[0.5*((28^2)+(0.5*(89^2)))];
thwi ={ sigma1d*Sqrt[8/(0.5*sel[[1]])],sigma1d*Sqrt[8/(0.5*sel[[2]])],sigma1d*Sqrt[8/(0.5*sel[[3]])]}

simwidths8n = {Around[895.33,34.58],Around[av0158,errorwi0158],Around[av028,errorwi028]};

simwidths5n = {Around[av015,err015],Around[avwi0155,err0155] ,Around[avwi025,err025]}

TTest[{wi0155,wi0158}]

TTest[{wi015,wi018}]

wvss= ListLogLogPlot[{Transpose[{sel, simwidths8n}],Transpose[{sel,thwi}],Transpose[{sel,simwidths5n}]},Joined ->True,PlotStyle ->{Green, Red, Orange}, PlotLegends->{"Simulation 8000 N","Theoretical Predictions", "Simulation 5262N"}, AxesLabel->{"Selection","Width"}]

Realized vs . Proposed Dispersal

dispcalcnew::usage="Given the spatial pedigree gg, dispcalcnew[gg,i] calculates seed and pollen dispersal distances in generation i";

dispcalcnew[gg_,i_]:=Module[{parind=Transpose[gg[[i,2]]],parloc,seed,poll},
parloc=gg[[i-1,3]][[#]]&/@parind;
seed=Norm[#]&/@(parloc[[1]]-gg[[i,3]]);
poll=Norm[#]&/@ Flatten[Differences[parloc],1];
{seed,poll}]

sdSeeds=Sqrt[0.5*(28^2)];
xSeeds=RandomVariate[NormalDistribution[0,sdSeeds],{1000+1,15000}];
ySeeds=RandomVariate[NormalDistribution[0,sdSeeds],{1000+1,15000}];
seeds=Sqrt[xSeeds^2+ySeeds^2];
sdPolls=Sqrt[0.5*(89^2)];
xPolls=RandomVariate[NormalDistribution[0,sdPolls],{1000+1,15000}];
yPolls=RandomVariate[NormalDistribution[0,sdPolls],{1000+1,15000}];
polls=Sqrt[xPolls^2+yPolls^2];

propseed=HistogramDistribution[seeds[[1]]];proppoll=HistogramDistribution[polls[[1]]];

realdisp=dispcalcnew[s015n800rep1,#]&/@Range[2,1001];

realseed=HistogramDistribution[Flatten[realdisp[[All,1]]]];

realpoll=HistogramDistribution[Flatten[realdisp[[All,2]]]];

Seed dispersal distribution

DiscretePlot[{PDF[realseed,x],PDF[propseed,x]},{x,0,200},PlotRange->All,PlotLegends->{"Simulation","Proposed"},AxesLabel->{"Distance (m)","Probability"}]

Export["SeedDisp_small.pdf",%97]

Pollen dispersal distribution

DiscretePlot[{PDF[realpoll,x],PDF[proppoll,x]},{x,0,200},PlotRange->All,PlotLegends->{"Simulation","Proposed"},AxesLabel->{"Distance (m)","Probability"}]

Export["polldisp_small.pdf",%98]

{{"","Mean","StandardDeviation"},{"Proposed seed",Mean[seeds[[1]]],StandardDeviation[seeds[[1]]]},{"Proposed pollen",Mean[polls[[1]]],StandardDeviation[polls[[1]]]},{"Realised seed- pedigree 1", Mean[realseed],StandardDeviation[realseed]},{"Realised pollen- pedigree 1",Mean[realpoll],StandardDeviation[realpoll]}}//TableForm

Export["TableDisp_small.jpg",%179]

Low Density Barrier

barrierrep1 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\Barrier\\barrier_015_8000_1.csv"]];

barrierrep2 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\Barrier\\barrier_015_8000_2.csv"]];

barrierrep3 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\Barrier\\barrier_015_8000_3.csv"]];

barrierrep4 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\Barrier\\barrier_015_8000_4.csv"]];

b1 = barrierrep1[[-1]];

b2 = barrierrep2[[-1]];

b3 = barrierrep3[[-1]];

b4 = barrierrep4[[-1]];

sigmoidb1=makeMetSNP[b1,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

clb1=mleMetSig[sigmoidb1,2500]

sigmoidb2=makeMetSNP[b2,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

clb2=mleMetSig[sigmoidb2,2500]

stepped

bstep1 =  makeMetasym[b1,Log[{2000,2000,0.05,0.05,2000,15000,0.1}],ConstantArray[0.1,7],350,1.01,2000];

bstep2 =makeMetasym[b2,Log[{2000,2000,0.05,0.05,2000,15000,0.1}],ConstantArray[0.1,7],350,1.01,2000];

bstep3 =  makeMetasym[b3,Log[{2000,2000,0.05,0.05,2000,15000,0.1}],ConstantArray[0.1,7],350,1.01,2000];

bstep4 =  makeMetasym[b4,Log[{2000,2000,0.05,0.05,2000,15000,0.1}],ConstantArray[0.1,7],350,1.01,2000];

clb1=mleMetasym[bstep1,1000]

clb2=mleMetasym[bstep2,1000]

clb3=mleMetasym[bstep3,1000]

clb4=mleMetasym[bstep4,1000]

centerbarrier = {12448.2,12409,12498.9,12424.9}

barrierwi = {742.485,506.391,699.825,668.826}

avwistepped = Mean[barrierwi]

errrorwistepped = StandardDeviation[barrierwi]/Sqrt[Length[barrierwi]]

TTest[{wi0158,barrierwi}]

wi0158

avecenteruni = Mean[centeruni]//N

avcenb = Mean[centerbarrier]

midsmall = 422775;

center = midsmall - X0
StandardDeviation[centerbarrier]

StandardDeviation[centeruni]//N

Long - Range - Dispersal

unilongdisp1 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\Dispersal\\smallarea_8000_01_allgen.csv"]];

unidisp1 = unilongdisp1[[-1]];

disp1 = ToExpression[Import["C:\\Users\\mjulseth\\Documents\\Barton\\Dispersal\\datadispersal_01_8000.csv"]];

steppeddisp1 =  makeMetasym[disp1,Log[{2000,2000,0.05,0.05,2000,15000,0.1}],ConstantArray[0.1,7],350,1.01,2000];

cldisp1stepped=mleMetasym[steppeddisp,1000]

sigmoid1unidisp1=makeMetSNP[unidisp1,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl1dispuni=mleMetSig[sigmoid1unidisp,2500]

sigmoid1disp=makeMetSNP[disp1,{Log[2000],Log[15000],Log[0.1]},ConstantArray[0.1,3],350,1.01,3000];

cl1disp=mleMetSig[sigmoid1disp,2500]

Code

Definitions

nearestF and nearest20

nearestF::usage="nearestF[inds] stores a function that gives the closest individuals to a given point. nearestF[{{(Subscript[x, 1]),(Subscript[y, 1])},…}] gives the same, but indexing individuals in numerical order corresponding to the list of positions";

nearestF[xy:{{_,_}...}]:=nearestF[xy]=Nearest[xy->Range[Length[xy]]];

nearest20::usage= 
"nearst20[g_] stores a function thate gives the closest 20 individuals to each individual of one generation"

nearest20[g_]:= Module[{inds,xy},
inds = Range[Length[Last@g]];
xy = Last@g;
Nearest[xy->Range[Length[inds]],xy[[#]],21]&/@inds]

findmom

findmom::usage="findmom[kids,pars, {j,k},seed] chooses k points on a circle at a distance seed[[i]] for kid i, and takes the j individuals nearest the circle." ;

findmom[kids:{{_,_}...},pars_List,{j_Integer,k_Integer},seed_List]:=findmom[kids,pars,{j,k},seed]=(findmom[#1[[1]],pars,#1[[2]],{j,k}]&)/@Transpose[{kids,seed}];
findmom[kidind:{_,_},pars_List,r_,{j_Integer,k_Integer}]:=Module[{θ=2 π (Range[0,1-1/k,1/k]+Random[]),y,rnf,dl},y=Transpose[r {Cos[θ],Sin[θ]}];rnf=Flatten[nearestF[pars][(#1+kidind&)/@y]];dl=(Abs[Norm[pars[[#1]]-kidind]-r]&)/@rnf;Last/@Take[Sort[Transpose[{dl,rnf}]],1]];

almostfitdad and fittestdad

almostfitdad::usage = "almostfitdad[moms,pars, {j,k},seed] chooses k points on a circle at a distance poll[[i]] for mom i, and stores the j individuals nearest the circle." ;

almostfitdad[moms_List, pars_List,{j_Integer,k_Integer},poll_List]:=
almostfitdad[moms,pars,{j,k},poll]=almostfitdad[#[[1]],pars,#[[2]],{j,k}]&/@Transpose[{moms,poll}];
almostfitdad[mom_Integer, pars_List,r_,{j_Integer,k_Integer}]:=
Module[{θ=2π (Random[]+Range[0,1-1/k,1/k]),y,rnf,fitlist},
y=Transpose[r{Cos[θ],Sin[θ]}];
rnf=DeleteCases[#,mom][[1]]&/@nearestF[pars][(pars[[mom]]+#)&/@y,2]]

fittestdad::usage ="fittestdad[dd,w] randomly chooses and individual weighted by it's fitnes value out of all possible inndividuals stored in almostfitdad"

fittestdad[dd_,w_]:=Module[{adad,fit,dads},
adad = DeleteDuplicates[dd[[#]]]&/@Range[Length[dd]];
fit = w[[#]]&/@adad;
dads = RandomChoice[fit[[#]]->adad[[#]]]&/@Range[Length[adad]]]

pfitnew

pfitnew::usage="pfitnew[kidsX,g,seed,poll,w] find the corresponding parents of the prior generation g to each kidX based on seed and pollen dispersal 
and relative fitness w, when choosing the fathers"

pfitnew[kidsX_List,g_List,seed_List,poll_List,w_]:=Module[{pars,moms,a,fitdad},
pars = Last@g;
moms=Flatten[findmom[kidsX,pars,{1,6},seed]];
a = almostfitdad[moms,pars,{1,12},poll];
fitdad = fittestdad[a,w];
Transpose[{moms,fitdad}]];

relfitfreq

reflitfreq::usage = "relfit[g,s] calculates ands stores the relative fitness of each individual in a generation g based on the genotype of it's 20 closes neightbours
 and positive frequency selection with the selection coefficent s
 afer Maller and Barton";

relfit[g_List,s_]:=Module[{inds,near,neighbours,l,gen,loci,p},
inds=Range[Length[Last[g]]];
near=nearest20[g];
neighbours=Drop[near,None,1];
gen=First[g];
loci=(gen[[#1]]&)/@neighbours;
loci=(#1/. {{0,0}->0,{0,1}->1,{1,0}->1,{1,1}->2}&)/@loci;
p=(((Total[#1])/(2 Length[#1]))&)/@loci;
l=(gen[[#1]]&)/@inds;
l=(#1/. {{0,0}->0,{0,1}->1,{1,0}->1,{1,1}->2}&)/@l;
MapThread[If[#1==1,1,If[#1==0,1-(#2-0.5) s,(#2-0.5) s+1]]&,{l,p}]]

relfit

reflit::usage = "relfit[g,s] calculates ands stores the relative fitness of each individual in a generation g with selection s against heterozygotes";

relfit[g_List,s_]:=Module[{inds,near,neighbours,l,gen,loci,p},
inds=Range[Length[Last[g]]];
gen = First@g;
loci=(gen[[#1]]&)/@inds;
loci=(#1/. {{0,0}->0,{0,1}->1,{1,0}->1,{1,1}->2}&)/@loci;
If[#1==1,1-s,1]&/@loci]

newoffgen

newoffgen ::usage = "newoffgen[pars,g]generates offpsring genotype from parental genotype after mendelian inheritance rule"

f[a_,b_]:=a[[b]]

newoffgen[pars_List,g_List]:=Module[{p,gen,pargen,genpos,kgen},
p = Flatten[pars];
gen = First[g];
pargen = gen[[#]]&/@p;
genpos= RandomInteger[{1,2},2 Length[pars]];
Partition[MapThread[f][{pargen,genpos}],2]]

newfititerFReal

newfititerFReal::usage="newfititerF[g,kidsX,seed,poll,s] iterates the population back through one generation, with seed dispersal values seed, pollen dispersal values poll and realtive fitness value w.  
For individuals with locations kidsX:{{x2,y2},{x2,y2}}, newfititerFReal finds the parents from the previous generation g, with locations given by {{x1,y1},{x1,y1}}. 
Returns a list giving the parents (par) and location of all indviduals that make up that generation."

newfititerFReal[g_List,kidsX_List,seed_List,poll_List,s_]:=
Module[{w,parents,kidsgen},
w = relfit[g,s];
parents = pfitnew[kidsX,g,seed,poll,w];
kidsgen = newoffgen[parents,g];
{kidsgen,parents,kidsX}]

makeN

makeN::usage="makeN[k,R1,R2] makes a list of coordinates of individuals, uniformly on {{0,R1},{0,R2}}";
makeN[k_Integer,R1_,R2_]:=Transpose[{RandomReal[{R1[[1]],R1[[2]]},k],RandomReal[{R2[[1]],R2[[2]]},k]}];

makeGenealogyRealuni

makeGenealogyReal::usage="[{(Subscript[x, 1,])...},t, lengthX, lengthY] iterates a pouplation forward in time through lengthX* lengthY habitat where individuals have a unifrom density. In each generation, n kids are uniformaly scattered over the space; two parents are chosen based on the seed and pollen dispersal distances and for fathers relativ fitness w with seleciton s ";

time=0;Dynamic[time]

makeGenealogyRealuni[x0_List,tm_Integer,Rx_,Ry_,startgen_,s_]:=
Module[{n=Length[x0],newg0,i,sdSeeds,xSeeds,ySeeds,xPolls,yPolls,sdPolls, seeds, polls,allruns,finalrun},
time= 0;
i=0;
sdSeeds=Sqrt[0.5*28^2];
xSeeds=RandomVariate[NormalDistribution[0,sdSeeds],{tm+1,n}];
ySeeds=RandomVariate[NormalDistribution[0,sdSeeds],{tm+1,n}];
seeds=Sqrt[xSeeds^2+ySeeds^2];
sdPolls=Sqrt[0.5*89^2];
xPolls=RandomVariate[NormalDistribution[0,sdPolls],{tm+1,n}];
yPolls=RandomVariate[NormalDistribution[0,sdPolls],{tm+1,n}];
polls=Sqrt[xPolls^2+yPolls^2];
newg0={startgen,ConstantArray[{},n],x0};
allruns=NestList[(time+=1;i++;newfititerFReal[#,makeN[n,Rx,Ry],seeds[[i]],polls[[i]],s])&,newg0,tm]]

makeNbarrier, makeGenealogyBarrier

makeNbarrier::usage="makeNBarrier, makes a list of coordinates of individuals with uniform distribution, but with a 10 fold lower density in a center area +/-80m from the midpoint ";

makeNbarrier[nsmall_,Rx_, Ry_, midpoint_] := Module[{area, originaldensity,
centerarea,outerarea,nout,left,center,right,nlowdensity,centerRx,leftbrx, rightbrx, ncenter,nleft,nright},
area= (Rx[[2]]-Rx[[1]])*(Ry[[2]]-Ry[[1]]);
originaldensity = nsmall/area;
nlowdensity = originaldensity*0.1;
centerRx = {midpoint-80,midpoint+80};
centerarea = (Ry[[2]]-Ry[[1]])*(centerRx[[2]]-centerRx[[1]]);
ncenter = centerarea*nlowdensity;
leftbrx = {Rx[[1]],Rx[[1]]+(centerRx[[1]]-Rx[[1]])};
rightbrx = {centerRx[[2]],Rx[[2]]};
outerarea = (Ry[[2]]-Ry[[1]])* (centerRx[[1]]-Rx[[1]]);
nout = (nsmall-IntegerPart[ncenter])/2;
left= makeN[IntegerPart[nout],leftbrx,Ry];
center = makeN[IntegerPart[ncenter],centerRx,Ry];
right = makeN[IntegerPart[nout],rightbrx,Ry];
Partition[Flatten[{left,center,right}],2]]

makeGenealogyRealbarrier::usage="[{(Subscript[x, 1,])...},t, lengthX, lengthY] iterates a pouplation forward in time through lengthX* lengthY habitat where individuals have a 10 fold lower density in the center area. In each generation, n kids are uniformaly scattered over the space; two parents are chosen based on the seed and pollen dispersal distances and for fathers relativ fitness w with seleciton s ";

makeGenealogyRealbarrier[x0_List,tm_Integer,Rx_,Ry_,startgen_,s_,midpoint_]:=
Module[{n = Length[x0],newg0,i,sdSeeds,xSeeds,ySeeds,xPolls,yPolls,sdPolls, seeds, polls,allruns,finalrun,nsmall},
time= 0;
nsmall = Length[x0];
i=0;
sdSeeds=Sqrt[0.5*38^2];
xSeeds=RandomVariate[NormalDistribution[0,sdSeeds],{tm+1,n}];
ySeeds=RandomVariate[NormalDistribution[0,sdSeeds],{tm+1,n}];
seeds=Sqrt[xSeeds^2+ySeeds^2];
sdPolls=Sqrt[0.5*119^2];
xPolls=RandomVariate[NormalDistribution[0,sdPolls],{tm+1,n}];
yPolls=RandomVariate[NormalDistribution[0,sdPolls],{tm+1,n}];
polls=Sqrt[xPolls^2+yPolls^2];
newg0={startgen,ConstantArray[{},n],x0};
allruns=NestList[(time+=1;i++;newfititerFReal[#,makeNbarrier[nsmall,Rx,Ry,midpoint],seeds[[i]],polls[[i]],s])&,newg0,tm]]

dispData,seeddips,polldisp,makeGenealogyRealdisp

Set Directory to where your data is stored

SetDirectory["***"];
dispData=Import["Dispersaldata.xlsx"];

disp::usage="disp lists seed & pollen dispersal distances from the pedigree, for 1463 trios";
disp=dispData[[1]][[2;;-1,{1,2}]];

seeddisp=First/@disp;

polldisp=Last/@disp;

makeGenealogyRealdisp[x0_List,tm_Integer,Rx_,Ry_,startgen_,s_]:=
Module[{n=Length[x0],newg0,i,sdSeeds,xSeeds,ySeeds,xPolls,yPolls,sdPolls, seeds, polls,allruns,finalrun},
time= 0;
i=0;
seeds= RandomVariate[HistogramDistribution[seeddisp],{tm+1,n}];
polls=RandomVariate[HistogramDistribution[polldisp],{tm+1,n}];
newg0={startgen,ConstantArray[{},n],x0};
allruns=NestList[(time+=1;i++;newfititerFReal[#,makeN[n,Rx,Ry],seeds[[i]],polls[[i]],s])&,newg0,tm];
finalrun = allruns[[-1]]]

Checks

Proposed vs Realised dispersal distance:

dispcalcnew

dispcalcnew::usage="Given the spatial pedigree gg, dispcalcnew[gg,i] calculates seed and pollen dispersal distances in generation i";

dispcalcnew[gg_,i_]:=Module[{parind=Transpose[gg[[i,2]]],parloc,seed,poll},
parloc=gg[[i-1,3]][[#]]&/@parind;
seed=Norm[#]&/@(parloc[[1]]-gg[[i,3]]);
poll=Norm[#]&/@ Flatten[Differences[parloc],1];
{seed,poll}]

propseed=HistogramDistribution[seeddispnorm[[1]]];proppoll=HistogramDistribution[polldispnorm[[1]]];

realdisp=dispcalcnew[gg2,#]&/@Range[2,1001];

realseed=HistogramDistribution[Flatten[realdisp[[All,1]]]];

realpoll=HistogramDistribution[Flatten[realdisp[[All,2]]]];

Seed dispersal distribution

DiscretePlot[{PDF[realseed,x],PDF[propseed,x]},{x,0,200},PlotRange->All,PlotLegends->{"Simulation","Proposed"},AxesLabel->{"Distance (m)","Probability"}]

Pollen dispersal distribution

DiscretePlot[{PDF[realpoll,x],PDF[proppoll,x]},{x,0,200},PlotRange->All,PlotLegends->{"Simulation","Proposed"},AxesLabel->{"Distance (m)","Probability"}]

{{"","Mean","StandardDeviation"},{"Proposed seed",Mean[seeddispnorm[[1]]],StandardDeviation[seeddispnorm[[1]]]},{"Proposed pollen",Mean[polldispnorm[[1]]],StandardDeviation[polldispnorm[[1]]]},{"Realised seed- pedigree 1", Mean[realseed],StandardDeviation[realseed]},{"Realised pollen- pedigree 1",Mean[realpoll],StandardDeviation[realpoll]}}//TableForm

Theoretical Cline

ThCline

ThCline = 1/2 (1+Tanh[((Sqrt[(0.5*s)/(2*(σ^2))])(x-clinecenter))]);

ThClinePlot

ListLinePlot[Transpose[{x,ThCline}],PlotStyle->{Red,Dashed}, PlotLegends->"Red Dashed, Theoretical Cline Prediction"]

Definitions Cline Fitting

clump::usage="clump[δ] stores a list of individuals for the entire dataset that are within a radius δ of the seed positions - by default, the positions of each individual.  Track progress using Dynamic[{i,Length[inds]}]. Individuals with missing positions not included. 
clump[δ, X, Y, Ninds] does the same with specified locations and number of individuals";

clump[δ_,X_,Y_,NInds_]:=clump[δ,X,Y,NInds]=Module[{XYI,inds,clxy,clxyInds,nf,clusterInds,nn, lengthInds},
XYI=Transpose[{X,Y, Range[NInds]}];
	clxyInds=XYI[[All,3]];(* individuals that correspond to each position *)
clxy=XYI[[All,{1,2}]];(* a list of positions *)
nf=Nearest[clxy->clxyInds];(* function that finds the individuals nearest a given position *)
inds=clxyInds;(* the individuals that remain to be allocated; should be 0 at the end *)
lengthInds=Length[inds];
clusterInds=ConstantArray[{},Length[clxy]];(* initial clusters*)
Do[nn=Intersection[nf[clxy[[i]],{∞,δ}],inds];clusterInds[[i]]=Join[clusterInds[[i]],nn];inds=Complement[inds,nn],{i,Length[clxy]}];
	Union/@DeleteCases[clusterInds,{}]];

zwjnTable

zwjnTable::usage"gives a table of Mean1D position of demes, totalf allelefreq per deme and freq of allele type 1 per deme"

zwjnTable[data_,δ_]:=Module[{deme,MeanXY,gen,genperdeme,n,j,z},
deme=clump[δ,data[[3,All,1]],data[[3,All,2]],Length[data[[3]]]];
MeanXY= Mean[data[[3]][[#]]]&/@deme;
gen=data[[1]]/.{{0,0}->0,{1,0}->1,{0,1}->1,{1,1}->2} ;
genperdeme=gen[[#]]&/@ deme;
n = 2Length[#] &/@genperdeme;
j= Total[#]&/@genperdeme;
z=MeanXY[[All,1]]-X0;
Transpose[{z,n,j}]]

makeMet

makeMet::usage="makeMet[L,α,δα,λ,n] makes a run of n Metropolis trials, starting with initial parameters α, with perturbations {δα,…}. The Metropolis chain uses log(parameters) for width, centre and Fst, uniformly perturbed in {-δ,δ}; perturbations are adjusted by a factor λ or (λ^(-1)), on acceptance or rejection. Log likelihood is given by L[α]. Perturbations can be set to zero, fixing the corresponding parameters. Returns a list of {L,α,δlog[α]} ";

makeMet[L_,α0_List,δα_List,λ_,n_Integer]:=NestList[iterMet[L,λ],{L[α0],α0,δα},n];

makeMet::usage="makeMetSNP stores a run of n Metropolis trials, based on data from zwjnTable[k,δ][[SNP,…]]. The barrier strengths are given an upper bound.  Tracked by {time,lik}. Perturbations can be set to zero, fixing the corresponding parameters. Returns a list of {L,α,δlog[α]}. makeMetSNP[{w,y,F},{δlog(w),…},SNP,k,δ,λ,n] fits a sigmoid cline";

makeMetSNP[data_,{w_, y_,f_},δα_List,δ_,λ_,n_Integer]:=
makeMetSNP[{w,y,f},δα,δ,λ,n]=Module[{L,znj},
znj=zwjnTable[data,δ];time=0;
L[{lw_,ly_,lf_}]:=likCline[znj,{0,1,Exp[lw],Exp[ly],Exp[lf]}];
makeMet[L,{w,y,f},δα,λ,n]];

makeMet[data_{p_,q_,w_, y_,f_},δα_List,δ_,λ_,n_Integer]:=
makeMetSNP[data,{p,q,w,y,f},δα,δ,λ,n]=Module[{L,znj},
znj=zwjnTable[data,δ];time=0;
L[{p0_,p1_,lw_,ly_,lf_}]:=likCline[znj,{p0,p1,Exp[lw],Exp[ly],Exp[lf]}];
makeMet[L,{p,q,w,y,f},δα,λ,n]];

makeMetasym[data_{b0_,b1_,θ0_,θ1_,w_,y_,f_},δα_List,δ_,λ_,n_Integer]:=
makeMetasym[data,{b0,b1,θ0,θ1,w,y,f},δα,δ,λ,n]=Module[{L,znj},
znj=zwjnTable[data,δ];time=0;
L[{lb0_,lb1_,lθ0_, lθ1_,lw_,ly_,lf_}]:=If[pb0[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]]>0.5∨ pb1[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]]<0.5∨pb0[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]]> pb1[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]] ,-∞,likbCline[znj,Exp[{lb0,lb1,lθ0, lθ1,lw,ly,lf}]]];
makeMet[L,{b0,b1,θ0,θ1,w,y,f},δα,λ,n]];

mleMet

mleMetSig::usage=" gives the MLE in the form {L,{(Subscript[B, 0]),…},{δlog((Subscript[B, 0])),…}}. Trials Δ;;-1 are used.";
mleMetSig[met_List,Δ_]:=Module[{met1=met[[Δ;;-1]]},
met1[[Position[met1,Max[met1[[All,1]]]][[1,1]]]]/.{L_,{w_,y_,f_},δα_}:>{L,{Exp[w],Exp[y],Exp[f]},δα}];

mleMetgenSig::usage="{L,{(Subscript[B, 0]),…},{…}},…},Δ] gives the MLE in the form {L,{(Subscript[B, 0]),…},{δlog((Subscript[B, 0])),…}}. Trials Δ;;-1 are used.";
mleMetgenSig[met_List,Δ_]:=Module[{met1=met[[Δ;;-1]]},
met1[[Position[met1,Max[met1[[All,1]]]][[1,1]]]]/.{L_,{p_,q_,w_,y_,f_},δα_}:>{L,{p,q,Exp[w],Exp[y],Exp[f]},δα}];

mleMetasym[met_List,Δ_]:=Module[{met1=met[[Δ;;-1]]},
met1[[Position[met1,Max[met1[[All,1]]]][[1,1]]]]/.{L_,{b0_,b1_,θ0_,θ1_,w_,y_,f_},δα_}:>{L,Exp[{b0,b1,θ0,θ1,w,y,f}],δα}];

likCline

likCline::usage="likCline[{{(Subscript[x, 1]),(Subscript[n, 1]),(Subscript[j, 1])},…},{(Subscript[p, 0]),(Subscript[p, 1]),w,y,F}] gives the log likelihood of a sigmoid cline. Using {(Subscript[p, 0]),(Subscript[p, 1]),λ0,λ1,w,y,F} gives exponential tails; (Subscript[λ, i]) must be larger than w.";
likCline[{x_,n_Integer,j_Integer},pars_List]:=Module[{},
Log[pr[Last[pars],pCline[Drop[pars,-1]][x],j,n]]];
likCline[xnj:{{_,_Integer,_Integer}...},pars_List]:=
Total[DeleteCases[likCline[#,pars]&/@xnj,Indeterminate]];

pCline

pCline::usage="pCline[\!\(TraditionalForm\`\(\({\*SubscriptBox[\(p\), \(0\)], \*SubscriptBox[\(p\), \(1\)], w, y}\)\(]\)\)[x]\\ is\\ a\\ sigmoid\\ cline, \\ with\\ allele\\ frequencies\\ \*SubscriptBox[\(p\), \(0\)], \*SubscriptBox[\(p\), \(1\)]\\ at\\ either\\ side . \\ \(pCline[{\*SubscriptBox[\(p\), \(0\)], \*SubscriptBox[\(p\), \(1\)], \*SubscriptBox[\(λ\), \(0\)], \*SubscriptBox[\(λ\), \(1\)], w, y}]\)[x]\\ has\\ exponential\\ tails\\ of\\ the\\ form\\ \*SubscriptBox[\(p\), \(i\)] exp[\(-x\)/\*SubscriptBox[\(λ\), \(i\)]]\\ on\\ either\\ side\)";
pCline[{p0_,p1_,w_,y_}][x_]:=Module[{},
p0+(((p1-p0))/(1+Exp[-4(x-y)/w]))];
pCline[{w_,y_}][x_]:=Module[{p0=0,p1=1},
p0+(((p1-p0))/(1+Exp[-4(x-y)/w]))];
pCline[{p0_,p1_,λ0_,λ1_,w_,y_}][x_]:=Module[{},
If[x≤ y,Max[p0 Exp[+(((x-y))/(λ0))],(1/(1+Exp[-4(x-y)/w]))],
Min[1-(1-p1) Exp[-(((x-y))/(λ1))],(1/(1+Exp[-4(x-y)/w]))]]];

pbCline

makeMetasym[data_{b0_,b1_,θ0_,θ1_,w_,y_,f_},δα_List,δ_,λ_,n_Integer]:=
makeMetasym[data,{b0,b1,θ0,θ1,w,y,f},δα,δ,λ,n]=Module[{L,znj},
znj=zwjnTable[data,δ];time=0;
L[{lb0_,lb1_,lθ0_, lθ1_,lw_,ly_,lf_}]:=If[pb0[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]]>0.5∨ pb1[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]]<0.5∨pb0[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]]> pb1[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]] ,-∞,likbCline[znj,Exp[{lb0,lb1,lθ0, lθ1,lw,ly,lf}]]];
makeMet[L,{b0,b1,θ0,θ1,w,y,f},δα,λ,n]];

pbCline::usage="pbCline\!\(TraditionalForm\`\(\(\\ \)\(is\\ a\\ better\\ representation\\ of\\ stepped\\ cline\\ based\\ on\\ rates\\ of\\ introgression\\ on\\ either\\ sides\)\)\)";

pb0[{B0_,B1_,θ0_,θ1_,w_,y_}]:=((B1 w (Sqrt[θ1]))/(B0 w (Sqrt[θ0])+B1 w (Sqrt[θ1])+4 B0 B1 (Sqrt[θ0]) (Sqrt[θ1])));
pb1[{B0_,B1_,θ0_,θ1_,w_,y_}]:=1-((B0 w (Sqrt[θ0]))/(B0 (Sqrt[θ0]) (w+4 B1 (Sqrt[θ1]))+B1 w (Sqrt[θ1])));

pbCline[{B0_,B1_,θ0_,θ1_,w_,y_}][x_]:=pCline[{pb0[{B0,B1,θ0,θ1,w,y}],pb1[{B0,B1,θ0,θ1,w,y}],(w/(4(Sqrt[θ0]))),(w/(4(Sqrt[θ1]))),w,y}][x];

pbCline[{B_,θ_,w_,y_}][x_]:=pCline[{pb0[{B,B,θ,θ,w,y}],pb1[{B,B,θ,θ,w,y}],(w/(4(Sqrt[θ]))),(w/(4(Sqrt[θ]))),w,y}][x];

pr

pr::usage="pr[F,p,j,n] gives the probability of sampling j from n, with binomial sampling from an underlying Beta";
pr[_,_,_,0|0.]:=1;
pr[f_,p_,j_Integer,n_Integer]:=Module[{α=1/f-1},
If[p≤0,
If[j==0,1,0],
If[p≥1,If[j==n,1,0],
Binomial[n,j]((Beta[α p+j,α (1-p)+n-j])/(Beta[α p,α (1-p)]))]]];
pr[0|0.,p_,j_Integer,n_Integer]:=If[p≤0,
If[j==0,1,0],
If[p≥1,If[j==n,1,0],
Binomial[n,j](p^j)(((1-p))^(n-j))]];

likCline

likCline::usage="likCline[{{(Subscript[x, 1]),(Subscript[n, 1]),(Subscript[j, 1])},…},{(Subscript[p, 0]),(Subscript[p, 1]),w,y,F}] gives the log likelihood of a sigmoid cline. Using {(Subscript[p, 0]),(Subscript[p, 1]),λ0,λ1,w,y,F} gives exponential tails; (Subscript[λ, i]) must be larger than w.";
likCline[{x_,n_Integer,j_Integer},pars_List]:=Module[{},
Log[pr[Last[pars],pCline[Drop[pars,-1]][x],j,n]]];
likCline[xnj:{{_,_Integer,_Integer}...},pars_List]:=
Total[DeleteCases[likCline[#,pars]&/@xnj,Indeterminate]];

likbCline

likbCline::usage="likCline[{{(Subscript[x, 1]),(Subscript[n, 1]),(Subscript[j, 1])},…},{(Subscript[p, 0]),(Subscript[p, 1]),w,y,F}] gives the log likelihood of a stepped cline. Using {(Subscript[p, 0]),(Subscript[p, 1]),λ0,λ1,w,y,F} gives exponential tails; (Subscript[λ, i]) must be larger than w.";
likbCline[{x_,n_Integer,j_Integer},pars_List]:=Module[{},
Log[pr[Last[pars],pbCline[Drop[pars,-1]][x],j,n]]];
likbCline[xnj:{{_,_Integer,_Integer}...},pars_List]:=
Total[DeleteCases[likbCline[#,pars]&/@xnj,Indeterminate]];

iterMet

iterMet::usage="iterMet[f,λ][{(OverHat[f]),α,δα}] iterates the Metropolis function through one step, making uniform random changes one by one to the parameters α in the range {-δα,δα}.  A change is accepted if Exp[(SuperStar[f])-f]>Random[]. If a change is accepted, the corresponding δα is increases by a factor λ>1, whilst if it is rejected, δα is decreased by λ. Tracks time, lik ";

iterMet[f_,λ_][{f0_,α_,δα_}]:=Module[{f1=f0,f2,α1=α,α2,δα1=δα,j,n=Length[α]},
lik=f0;Do[α2=α1;α2=ReplacePart[α2,j->α2[[j]]+δα1[[j]] RandomReal[{-1,1}]];f2=f[α2];(* make a random change to parameter j, yielding α2 *)
If[δα[[j]]>0,
If[Exp[f2-f1]>Random[],
δα1=ReplacePart[δα1,j-> δα1[[j]] λ];f1=f2;α1=α2,(* accept *)
δα1=ReplacePart[δα1,j->δα1[[j]]/λ]]],(* reject *)
{j,n}];time+=1;lik=f1;
{f1,α1,δα1}];

makeMet, makeMetSNP

makeMet::usage="makeMet[L,α,δα,λ,n] makes a run of n Metropolis trials, starting with initial parameters α, with perturbations {δα,…}. The Metropolis chain uses log(parameters) for width, centre and Fst, uniformly perturbed in {-δ,δ}; perturbations are adjusted by a factor λ or (λ^(-1)), on acceptance or rejection. Log likelihood is given by L[α]. Perturbations can be set to zero, fixing the corresponding parameters. Returns a list of {L,α,δlog[α]} ";

makeMet[L_,α0_List,δα_List,λ_,n_Integer]:=NestList[iterMet[L,λ],{L[α0],α0,δα},n];

makeMet::usage="makeMetSNP stores a run of n Metropolis trials, based on data from zwjnTable[k,δ][[SNP,…]]. The barrier strengths are given an upper bound.  Tracked by {time,lik}. Perturbations can be set to zero, fixing the corresponding parameters. Returns a list of {L,α,δlog[α]}. makeMetSNP[{w,y,F},{δlog(w),…},SNP,k,δ,λ,n] fits a sigmoid cline";

makeMetSNP[data_,{w_, y_,f_},δα_List,δ_,λ_,n_Integer]:=
makeMetSNP[{w,y,f},δα,δ,λ,n]=Module[{L,znj},
znj=zwjnTable[data,δ];time=0;
L[{lw_,ly_,lf_}]:=likCline[znj,{0,1,Exp[lw],Exp[ly],Exp[lf]}];
makeMet[L,{w,y,f},δα,λ,n]];

makeMet[data_,{p_,q_,w_, y_,f_},δα_List,δ_,λ_,n_Integer]:=
makeMet[data,{p,q,w,y,f},δα,δ,λ,n]=Module[{L,znj},
znj=zwjnTable[data,δ];time=0;
L[{p0_,p1_,lw_,ly_,lf_}]:=likCline[znj,{p0,p1,Exp[lw],Exp[ly],Exp[lf]}];
makeMet[L,{p,q,w,y,f},δα,λ,n]];

makeMetasym[data_{b0_,b1_,θ0_,θ1_,w_,y_,f_},δα_List,δ_,λ_,n_Integer]:=
makeMetasym[data,{b0,b1,θ0,θ1,w,y,f},δα,δ,λ,n]=Module[{L,znj},
znj=zwjnTable[data,δ];time=0;
L[{lb0_,lb1_,lθ0_, lθ1_,lw_,ly_,lf_}]:=If[pb0[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]]>0.5∨ pb1[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]]<0.5∨pb0[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]]> pb1[Exp[{lb0,lb1,lθ0, lθ1,lw,ly}]] ,-∞,likbCline[znj,Exp[{lb0,lb1,lθ0, lθ1,lw,ly,lf}]]];
makeMet[L,{b0,b1,θ0,θ1,w,y,f},δα,λ,n]];

mleMet

mleMetSig::usage=" gives the MLE in the form {L,{(Subscript[B, 0]),…},{δlog((Subscript[B, 0])),…}}. Trials Δ;;-1 are used.";
mleMetSig[met_List,Δ_]:=Module[{met1=met[[Δ;;-1]]},
met1[[Position[met1,Max[met1[[All,1]]]][[1,1]]]]/.{L_,{w_,y_,f_},δα_}:>{L,{Exp[w],Exp[y],Exp[f]},δα}];

mleMetgenSig::usage="{L,{(Subscript[B, 0]),…},{…}},…},Δ] gives the MLE in the form {L,{(Subscript[B, 0]),…},{δlog((Subscript[B, 0])),…}}. Trials Δ;;-1 are used.";
mleMetgenSig[met_List,Δ_]:=Module[{met1=met[[Δ;;-1]]},
met1[[Position[met1,Max[met1[[All,1]]]][[1,1]]]]/.{L_,{p_,q_,w_,y_,f_},δα_}:>{L,{p,q,Exp[w],Exp[y],Exp[f]},δα}];

mleMetasym[met_List,Δ_]:=Module[{met1=met[[Δ;;-1]]},
met1[[Position[met1,Max[met1[[All,1]]]][[1,1]]]]/.{L_,{b0_,b1_,θ0_,θ1_,w_,y_,f_},δα_}:>{L,Exp[{b0,b1,θ0,θ1,w,y,f}],δα}];

Usage:

clinefitsigmoid=makeMetSNP[{Log[w],Log[y],Log[f]},ConstantArray[Δu,3],l,k,δ,λ,n]; (*Use other makeMEt functions to fit stepped and general sigmoid*)

mleMetSig[clinefitsigmoid,500]

This will give you the MLE of parameters.

models[p0_,p1_,b0_,b1_,θ0_,θ1_,w_,y_,f_, l_,k_,δ_, λ_, n_ ,Δp_,Δu_]:= Module[{yys, yygenS, yysym,yyasym},
yys=makeMetSNP[{Log[w],Log[y],Log[f]},ConstantArray[Δu,3],l,k,δ,λ,n];
yygenS=makeMetSNP[{p0,p1,Log[w],Log[y],Log[f]},{Δp,Δp,Δu,Δu,Δu},l,k,δ,λ,n];
(*yysym= makeMetsymSNP[Log[{b0,θ0, w,y,f}], ConstantArray[Δu,5],l,k,δ,λ,n];*)
yyasym= makeMetasymSNP[Log[{b0,b1,θ0,θ1,w,y,f}], ConstantArray[Δu,7],l,k,δ,λ,n];
{yys,yygenS,(*yysym,*) yyasym}]
models[b0_,b1_,θ0_,θ1_,w_,y_,f_, l_,k_,δ_, λ_, n_ ,Δu_]:=makeMetasymSNP[Log[{b0,b1,θ0,θ1,w,y,f}], ConstantArray[Δu,7],l,k,δ,λ,n];

getmax[{sig_, gensig_,(*sym_,*)asym_}]:={mleMetSig[sig,500],mleMetgenSig[gensig,500],(*mleMetsym[sym,500], *)mleMetasym[asym,500]}

You dont need to use this, but to give an exampple, I have this function that fits all models .. to show an example of cline fitting in ROS,

rosall=models[0.1,0.9,20000,20000,0.05,0.05,2000,15000,0.1,187,6,25,1.01,3000,0.05,0.1];

getmax[rosall]

The former gives the likelihood for all runs, from which u can see how it changes with time , relation between paramaters and so on and the latter will give the paramaters and the best likelihood estimate