Recombination and the Undesirable Hitchhiker

Matthew Hartfield & Sarah P. Otto

Inputs [ENTER]

Off[General::Spell];
Off[General::spell1];

Clear["`*"];
Needs["ErrorBarPlots`"];

<<"MultivariateStatistics`";
<<"Combinatorica`"

Supplementary material 1:  Deriving the establishment probability for a recombinant

The probability of establishment for a recombinant

Let Π[t] stand for the fixation probability of the 10 haplotype appearing when the frequency of 11 is p[t].

Let selection for 10 be given by θ[t], which varies over time according to θ[t] = ( (Subscript[s, a]) - p[t] (Subscript[s, net])), where (Subscript[s, net]) = (Subscript[s, a]) - (Subscript[s, d]).

By standard branching process logic, the probability that the 10 haplotype is ultimately lost, 1-Π[t], is given by the sum over the distribution of possible offspring numbers multiplied by the probability that all offspring copies are lost:

Sum[PDF[PoissonDistribution[1+θ],i](((1-Π[t+1]))^i),{i,0,Infinity}]//Simplify

Thus, Π[t+1] = - Log[1-Π[t]]/(1+θ[t])].  The change in the fixation probability is then Π[t+1] - Π[t], which equals

-((Log[1-Π[t]])/(1+θ[t])) -Π[t];

Assuming weak selection and a correspondingly small fixation probability, this difference gives us ((d Π[t])/dt), which to leading order equals:

Normal[Series[%/.Π[t]->Π[t]*small/.θ[t]->θ[t]*small,{small,0,2}]]/.small->1//Simplify

Next, we use the chain rule to rewrite this derivative as ((d Π[p])/dp), noting that for the haploid model of weak selection (dp/dt)= (Subscript[s, net]) p (1-p):

(((1/2) Π[p] (Π[p]-2 θ[p]))/((Subscript[s, net]) p (1-p)))/.θ[p]-> (Subscript[s, a])-p  (Subscript[s, net])

Π[p]/.Flatten[DSolve[{D[Π[p],p]==((Π[p] (-2 ((Subscript[s, a])-p (Subscript[s, net]))+Π[p]))/(2 (1-p) p (Subscript[s, net])))},{Π[p]},p]]//Simplify

Finally, we assume all selection coefficients are weak:

Normal[Series[%/.(Subscript[s, a])->(Subscript[s, net])+(Subscript[s, d])/.(Subscript[s, net])->(Subscript[s, net])*small/.(Subscript[s, d])->(Subscript[s, d])*small,{small,0,1}]]/.small->1//Simplify

This is equation (2).

Note that it has the correct boundary behavior, approaching 2*(Subscript[s, a]) as p goes to zero and 2*(Subscript[s, d]) as p goes to one, where p is the frequency of the 11 haplotype when the recombinant 10 haplotype first appears.

Naively, we might have thought that the fixation probability of the recombinant 10 haplotype would be approximately twice its current selective advantage relative to the mean fitness of the population:  2( (Subscript[s, a])-p (Subscript[s, net])).  The following graphs show that the fixation probability derived above (solid curve) is substantially lower than this naive expectation (dashed curve).  That is, the fixation probability of the 10 recombinant is lower than expected based on its current fitness, because its future selective advantage drops as the 11 haplotype spreads.  Note that the discrepancy is most pronounced when the recombinant appears while the 11 haplotype is intermediate in frequency, which makes sense because p(t) changes most rapidly when near 1/2.

Show[Plot[((2 (Subscript[s, d]) ((Subscript[s, d])+(Subscript[s, net])))/((Subscript[s, d])+p (Subscript[s, net])))/.(Subscript[s, net])->0.04/.(Subscript[s, d])->0.01,{p,0,1},PlotRange->{{0,1},{0,0.1}}],
Plot[2( (Subscript[s, a])-p (Subscript[s, net]))/.(Subscript[s, a])->(Subscript[s, net])+(Subscript[s, d])/.(Subscript[s, net])->0.04/.(Subscript[s, d])->0.01,{p,0,1},PlotStyle->Dashing[{0.04}]],
AxesLabel->{"p","Π"}
]

Show[Plot[((2 (Subscript[s, d]) ((Subscript[s, d])+(Subscript[s, net])))/((Subscript[s, d])+p (Subscript[s, net])))/.(Subscript[s, net])->0.02/.(Subscript[s, d])->0.01,{p,0,1},PlotRange->{{0,1},{0,0.06}}],
Plot[2( (Subscript[s, a])-p (Subscript[s, net]))/.(Subscript[s, a])->(Subscript[s, net])+(Subscript[s, d])/.(Subscript[s, net])->0.02/.(Subscript[s, d])->0.01,{p,0,1},PlotStyle->Dashing[{0.04}]],
AxesLabel->{"p","Π"}
]

Solution to the integral in the semi-deterministic model (equation 4)

To solve equation (4), we take N times the chance of a successful recombinant:

part1=N*r/2*2*p*(1-p)*((2*(Subscript[s, d])*(Subscript[s, a]))/(((Subscript[s, a])-(Subscript[s, d]))*p+(Subscript[s, d])));

and divide this by dp/dt, which for the haploid model of selection is (assuming weak selection so that we can approximate the recursion by a differential equation):

part2=((Subscript[s, a])-(Subscript[s, d]))*p*(1-p);

Then integrating part1/part2:

Simplify[Integrate[Simplify[part1/part2],p],{(Subscript[s, a])>0,(Subscript[s, d])>0,0<p0<1,N>0}]

Technically the integral should be evaluated from p = p0 to p = 1-1/N:

((%/.p->1-1/N)-(%/.p->p0))//FullSimplify

But given that we assume a reasonably large population that N-1 is approximately N and p0 is small (relative to sd), we can equally take the integral from p = 0 to 1:

((%%/.p->1)-(%%/.p->0))//FullSimplify

The probability of a hitchhiking event can then be written as Exp[-((2 N r (Subscript[s, a]) (Subscript[s, d]) Log[((Subscript[s, a])/(Subscript[s, d]))])/((((Subscript[s, a])-(Subscript[s, d])))^2))].

Finally, we can solve for the critical rate of recombination, below which the hitchhiking to fixation occurs more than 50% of the time:

Simplify[Solve[Exp[-((2 N r (Subscript[s, a]) (Subscript[s, d]) Log[((Subscript[s, a])/(Subscript[s, d]))])/((((Subscript[s, a])-(Subscript[s, d])))^2))]==1/2,r]]

Taking the real solution:

rcritD=((Log[2] ((((Subscript[s, a])-(Subscript[s, d])))^2))/(2 N Log[((Subscript[s, a])/(Subscript[s, d]))] (Subscript[s, a]) (Subscript[s, d])));

which can be written in terms of scaled parameters as:

ρcritD=((Subscript[ψ, net])/(Subscript[ψ, d]))((Log[2])/(2 (1+((Subscript[ψ, d])/(Subscript[ψ, net])))Log[1+((Subscript[ψ, net])/(Subscript[ψ, d]))]));

The above has the form of  ((Subscript[ψ, net])/(Subscript[ψ, d]))C, where the C term hovers around 1/4 for values of interest:

Plot3D[((Log[2])/(2 (1+((Subscript[ψ, d])/(Subscript[ψ, net])))Log[1+((Subscript[ψ, net])/(Subscript[ψ, d]))])),{(Subscript[ψ, net]),2,100},{(Subscript[ψ, d]),2,100},AxesLabel->{(Subscript[ψ, net]),(Subscript[ψ, d]),C}]

Supplementary material 2:  Deriving the diffusion process with killing

In the standard one-locus diffusion model, the stochastic dynamics of an allele is approximately described  by a diffusion process, derived by calculating the mean change in allele frequency (= drift term) and the variance in the change in allele frequency (= diffusion term), as discussed in Box 15.2 of Otto and Day (2007) or Chapter 15 in Karlin and Taylor (1981).  This diffusion equation is then used to infer such things as the fixation probability of an allele.

Here, however, our initial "allele" is actually a 11 haplotype, and it can be broken apart by recombination, preventing fixation of the 11 type.

We thus use a diffusion process with killing, where killing occurs at rate κ and stops the diffusion process.  In the following, we first derive the probability that the process is killed (the 10 haplotype fixes carrying only the beneficial allele) versus the probability that it is not killed (11 or 00 fixes).

Note, in the following, we use (Subscript[ψ, i]) to stand in for the selection coefficient multiplied by the population size, N*(Subscript[s, i]), instead of (Subscript[S, i]) as used in the text.

Section 1.1:  Not conditioning on fixation of the (Subscript[A, 1]) allele

After one generation of selection, and assuming the process has not been killed, the expected frequency of the 11 haplotype is:

pprime=(((1+(Subscript[s, a])-(Subscript[s, d]))p)/(1+((Subscript[s, a])-(Subscript[s, d]))p));

But the actual number of 11 haplotypes (j) will depend on a draw from the Binomial Distribution:

part1a=(∑_(j=0)^N)((j-i)*PDF[BinomialDistribution[N, pnext ],j])//Simplify

Scaling the number of 11 haplotypes by (N^β) and scaling time by (N^α) (where α and β will be chosen to ensure that the process can be approximated by a diffusion model), we get the drift coefficient:

part1b=part1a*((N^α)/(N^β))/.i->(N^β)*iscaled

We also need the variance in the change in number of 11 haplotypes:

part2a=(∑_(j=0)^N)((((j-i))^2)*PDF[BinomialDistribution[N, pnext ],j])//Simplify

Scaling the number of 11 haplotypes by (N^β) and scaling time by (N^α), we get the diffusion coefficient:

part2b=part2a*((N^α)/(N^(2β)))/.i->(N^β)*iscaled

We must make choices for α and β so that the drift term is finite and the diffusion term is finite and positive, in the limit as N gets large.  In the absence of selection, α = β =1 works:

Simplify[part1b/.α->1/.β->1/.pnext->p/.iscaled->p,   {0<p<1,0<p (s-S)<1}]

FullSimplify[part2b/.α->1/.β->1/.pnext->p/.iscaled->p,   {0<p<1,0<p (s-S)<1}]

With selection, we must also assume that selection is weak, such that N*((Subscript[s, a])-(Subscript[s, d])) = (Subscript[ψ, net]) remains finite even as N gets large:

Simplify[part1b/.α->1/.β->1/.pnext->pprime/.iscaled->p,   {0<p<1,0<p (s-S)<1}];
Simplify[%/.(Subscript[s, a])->(Subscript[s, d])+(Subscript[ψ, net])/N]

Limit[%,N->∞]

FullSimplify[part2b/.α->1/.β->1/.pnext->pprime/.iscaled->p,   {0<p<1,0<p (s-S)<1}];
Simplify[%/.(Subscript[s, a])->(Subscript[s, d])+(Subscript[ψ, net])/N]

Limit[%,N->∞]

This gives us the standard drift and diffusion coefficients for the haploid model of selection:

	μ=(Subscript[ψ, net]) p(1-p)
(σ^2)=p(1-p)

Next, we figure out the killing term for a diffusion.  The probability that any one offspring is a 10 successful recombinant is:

	 κ = r p (1-p) Π 
	   = r p (1-p) ((2(Subscript[s, a])(Subscript[s, d]))/((Subscript[s, a])p + (Subscript[s, d]) (1-p))),

and the process will be killed if any one of the offspring is a successful recombinant:

part3=1-PDF[BinomialDistribution[N, κ ],0]//Simplify

We must assume that κ is sufficiently small that even N*κ is small, in which case, we can replace the above with the approximation:

part3a=N κ;

Scaling time by (N^α) and replacing for the κ gives, we get the killing term:

part3b=part3a*(N^α)/.κ->r p (1-p) ((2 (Subscript[s, a]) (Subscript[s, d]))/(p (Subscript[s, a])+(1-p) (Subscript[s, d])))

Simplify[%/.α->1/.(Subscript[s, a])->(Subscript[s, d])+(Subscript[ψ, net])/N/.(Subscript[s, d])->(Subscript[ψ, d])/N]

This will only be finite as N gets large if we assume that r is sufficiently weak that N r = ρ remains finite:

Simplify[%/.r->ρ/N]

The probability that the system reaches a boundary (fixation of 00 or 11) before being killed then solves (see Karlin & Taylor, p. 277):

notkilled[ψnet_,ψd_,ρ_]=P[p]/.Flatten[Simplify[DSolve[{((P''[p])/2)+ψnet*P'[p]-ρ*((2ψd(ψd+ψnet))/(p ψnet+ψd))*P[p]==0,P[0]==P[1]==1},P[p],p],{0<p<1,ψnet>0,ψd>0,ρ>0}]]

This is equation (B.5) in Appendix B.

(HypergeometricU[-ω,0,-2 (ψd+p ψnet)] (LaguerreL[ω,-1,-2 ψd]- LaguerreL[ω,-1,-2 (ψd+ψnet)])-LaguerreL[ω,-1,-2 (ψd+p ψnet)](HypergeometricU[-ω,0,-2 ψd] -HypergeometricU[-ω,0,-2 (ψd+ ψnet)]))/( (HypergeometricU[-ω,0,-2 (ψd+ψnet)] LaguerreL[ω,-1,-2 ψd]-HypergeometricU[-ω,0,-2 ψd] LaguerreL[ω,-1,-2 (ψd+ψnet)]))/.ω->((2 ρ ψd (ψd+ψnet))/((ψnet)^2));
%/%%//Simplify

Section 1.2:  Conditioning on fixation of the (Subscript[A, 1]) allele

Using the approach described on p. 262 and 263 of Karlin and Taylor (1982), we can condition on the probability that the system, if it remains unkilled, will reach the 11 not the 00 boundary.

For the drift term given by equation (9.5a) in Karlin and Taylor, we need:

conds[ξ_,ψ_]=Exp[-Integrate[2*ψ*x*(1-x)/(x*(1-x)),{x,0,ξ}]]

condS[ξ_,ψ_]=Integrate[conds[ξ ,ψ],{ξ,0,p}]

The new drift term is then:

(Subscript[ψ, net])*p*(1-p)+((conds[p,(Subscript[ψ, net])])/(condS[p,(Subscript[ψ, net])]))p*(1-p)//Simplify

The diffusion term is unaffected, as is the killing term (the probability of a recombinant occuring when the system is at frequency x is unaffected by what would have happened if the killing did not occur).

Accounting also for the changed boundary conditions, the probability that the system is not killed and ultimately reaches 11 is then given by the solution to:

notkilledC[ψnet_,ψd_,ρ_]=P[p]/.Flatten[Simplify[DSolve[{((P''[p])/2)+ψnet*(((1+(E^(-2 p ψnet)))/(1-(E^(-2 p ψnet)))))*P'[p]-ρ*((2ψd(ψd+ψnet))/(p ψnet+ψd))*P[p]==0,P'[0]==0,P[1]==1},P[p],p],{0<p<1,ψnet>0,ψd>0,ρ>0}]]

Mathematica does not solve this differential equation using the boundary condition P'[0], but it does if we assume that P'[c]== 0  for some small value of c and then let c go to 0:

P[p]/.Flatten[Simplify[DSolve[{((P''[p])/2)+ψnet*(((1+(E^(-2 p ψnet)))/(1-(E^(-2 p ψnet)))))*P'[p]-ρ*((2ψd(ψd+ψnet))/(p ψnet+ψd))*P[p]==0,P'[c]==0,P[1]==1},P[p],p],{0<p<1,ψnet>0,ψd>0,ρ>0}]]

notkilledC[ψnet_,ψd_,ρ_]=Simplify[%/.c->0]

This is equation (B.6) in Appendix B.

(((1-(E^(-2 ψnet))))/((1-(E^(-2 p ψnet))))) ( (HypergeometricU[-ω,0,-2 (ψd+p ψnet)] LaguerreL[ω,-1,-2 ψd]-HypergeometricU[-ω,0,-2 ψd] LaguerreL[ω,-1,-2 (ψd+p ψnet)]))/( (HypergeometricU[-ω,0,-2 (ψd+ψnet)] LaguerreL[ω,-1,-2 ψd]-HypergeometricU[-ω,0,-2 ψd] LaguerreL[ω,-1,-2 (ψd+ψnet)]))/.ω->((2 ρ ψd (ψd+ψnet))/((ψnet)^2));
%/%%//Simplify

Supplementary material 3:  Approximating the diffusion process with killing

Diffusion equations (from section 2) [ENTER]

notkilled[ψnet_,ψd_,ρ_]=P[p]/.Flatten[Simplify[DSolve[{((P''[p])/2)+ψnet*P'[p]-ρ*((2ψd(ψd+ψnet))/(p ψnet+ψd))*P[p]==0,P[0]==P[1]==1},P[p],p],{0<p<1,ψnet>0,ψd>0,ρ>0}]]

P[p]/.Flatten[Simplify[DSolve[{((P''[p])/2)+ψnet*(((1+(E^(-2 p ψnet)))/(1-(E^(-2 p ψnet)))))*P'[p]-ρ*((2ψd(ψd+ψnet))/(p ψnet+ψd))*P[p]==0,P'[c]==0,P[1]==1},P[p],p],{0<p<1,ψnet>0,ψd>0,ρ>0}]];

notkilledC[ψnet_,ψd_,ρ_]=Simplify[%/.c->0]

Approximation for the probability of hitchhiking (conditional)

ρ enters into the above in the same form in each term, so we define ω =((2 ρ (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d])))/((Subscript[ψ, net])^2)) to make the equations easier to read:

notkilledC[(Subscript[ψ, net]),(Subscript[ψ, d]),((ω ((Subscript[ψ, net])^2))/(2(Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))]

As we might expect, the conditional probability of hitchhiking declines roughly exponentially with the recombination rate (through ω):

plot1=Show[ListPlot[Table[{ω,Re[Block[{$MaxExtraPrecision=1000},N[Limit[notkilledC[40+10^(-16),10,((ω (40^2))/(2*10 (40+10)))],p->0],50]]]},{ω,1/10,5,1/4}],Joined->True,PlotRange->{{0,5},{0,1}}],
ListPlot[Table[{ω,Re[Block[{$MaxExtraPrecision=1000},N[Limit[notkilledC[20+10^(-16),10,((ω (20^2))/(2*10 (20+10)))],p->0],50]]]},{ω,1/10,5,1/4}],Joined->True,PlotRange->{{0,5},{0,1}}],
ListPlot[Table[{ω,Re[Block[{$MaxExtraPrecision=1000},N[Limit[notkilledC[10+10^(-16),10,((ω (10^2))/(2*10 (10+10)))],p->0],50]]]},{ω,1/10,5,1/4}],Joined->True,PlotRange->{{0,5},{0,1}}],
AxesLabel->{ω,(SuperStar[P])}
]

In the following, we assume that selection is strong relative to drift and find the value that the above approaches when (Subscript[ψ, d]) and (Subscript[ψ, net]) are large.

For large selection coefficients (so that c is always large in the HypergeometricU[a,0,c] terms), HypergeometricU[a,0,c] is approximately:

Simplify[Normal[Series[HypergeometricU[a,0,c/small],{small,0,0}]/.small->1],{Im[c]==0}]

while the LaguerreL is approximately:

Simplify[Normal[Series[LaguerreL[a,-1,c/small],{small,0,0}]],{Im[c]==0,Im[small]==0,c<0,small>0,a>0}]/.small->1//Simplify

Plugging in these limits into the conditional probability that the process is not killed gives:

notkilledC[(Subscript[ψ, net]),(Subscript[ψ, d]),((ω ((Subscript[ψ, net])^2))/(2(Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))]/.HypergeometricU[a_,0,c_]->c^(-a)/.LaguerreL[a_,-1,c_]->(((c^(-a)) ((((-(c^2)))^a) Gamma[-a]+(E^c) Gamma[a]))/(Gamma[-a] Gamma[1+a]))

We assume that the initial frequency of the 11 haplotype is 1/N and that this is very small.  We thus focus on the conditional probability given that the (Subscript[A, 1]) allele ultimately fixes, in the limit when p -> 0.

FullSimplify[Limit[%,p->0],{(Subscript[ψ, d])>0,(Subscript[ψ, net])>0,ω>0}]

This can be written as:

temp= (((1-(E^(-2 (Subscript[ψ, net]))))(1+(ω/(Subscript[ψ, d]))))/( ((((((Subscript[ψ, net])+(Subscript[ψ, d]))/(Subscript[ψ, d]))))^ω)(1- (E^(-2 (Subscript[ψ, net])))(((((Subscript[ψ, d]^2)/((((Subscript[ψ, d])+(Subscript[ψ, net])))^2))))^ω) )));
FullSimplify[temp/%%,{(Subscript[ψ, net])>0,(Subscript[ψ, d])>0,ω>0}]

This form compares well to the notkilledC[] function evaluated above (the worst fit occurs when selection is weakest, (Subscript[ψ, net])->10, top curves)

Show[plot1,
Plot[{temp/.(Subscript[ψ, net])->10/.(Subscript[ψ, d])->10,temp/.(Subscript[ψ, net])->20/.(Subscript[ψ, d])->10,temp/.(Subscript[ψ, net])->40/.(Subscript[ψ, d])->10},{ω,0,5},AxesLabel->{ω,P},PlotStyle->Dashing[{0.01}]]]

Assuming reasonably large (Subscript[ψ, net]) and recalling that  (Subscript[ψ, d]) is assumed smaller than (Subscript[ψ, a])=(Subscript[ψ, net])+(Subscript[ψ, d]), we drop (e^(-2 (Subscript[ψ, net]))) terms in the temp solution (each of which is small compared to one):

temp/.(E^(-2 (Subscript[ψ, net])))->0

To determine the best exponential approximation for the above, we take the Taylor series of the log of the above.  This gives us the following approximation for this exponential function:

FullSimplify[Exp[Simplify[Normal[Series[Log[(1+(ω/(Subscript[ψ, d]))) ((((((Subscript[ψ, d])+(Subscript[ψ, net]))/(Subscript[ψ, d]))))^(-ω))/.ω->ω*small],{small,0,1}]]]/.small->1],{(Subscript[ψ, net])>0,(Subscript[ψ, d])>0,ω>0}]

which is simply  ((((e^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω)), as reported in equation (12).  

Note that this approximation really requires that (ω/(Subscript[ψ, d])) be small (not ω), which is likely as long as the recombination rate is small enough relative to the strength of selection, which is true in all cases where there is a reasonable chance of hitchhiking to fixation.

Again, this performs well, with the worst fit occurring when selection is weaker (top curve):
[The long-dashed curve gives the semi-deterministic approximation for comparison.]

Show[plot1,
Plot[{((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, net])+(Subscript[ψ, d]))/(Subscript[ψ, d]))))^(-ω))/.(Subscript[ψ, net])->10/.(Subscript[ψ, d])->10,((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, net])+(Subscript[ψ, d]))/(Subscript[ψ, d]))))^(-ω))/.(Subscript[ψ, net])->20/.(Subscript[ψ, d])->10,((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, net])+(Subscript[ψ, d]))/(Subscript[ψ, d]))))^(-ω))/.(Subscript[ψ, net])->40/.(Subscript[ψ, d])->10},{ω,0,5},AxesLabel->{ω,P},PlotStyle->Dashing[{0.01}]],
Plot[{((((((Subscript[ψ, net])+(Subscript[ψ, d]))/(Subscript[ψ, d]))))^(-ω))/.(Subscript[ψ, net])->10/.(Subscript[ψ, d])->10,((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, net])+(Subscript[ψ, d]))/(Subscript[ψ, d]))))^(-ω))/.(Subscript[ψ, net])->20/.(Subscript[ψ, d])->10,((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, net])+(Subscript[ψ, d]))/(Subscript[ψ, d]))))^(-ω))/.(Subscript[ψ, net])->40/.(Subscript[ψ, d])->10},{ω,0,5},AxesLabel->{ω,P},PlotStyle->Dashing[{0.02}]]]

Solving for the critical value of recombination above which hitchhiking to fixation becomes unlikely

Using the above approximation, we find the recombination rate that gives a 50% chance of hitchhiking for the linked deleterious allele:

This can then be solved for the critical value of ρ:

ρ/.Simplify[Solve[(((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, d])+(Subscript[ψ, net]))/(Subscript[ψ, d]))))^(-ω))/.ω->((2 ρ (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d])))/((Subscript[ψ, net])^2)))==c,ρ][[2]],{0<c<1,(Subscript[ψ, d])>0,(Subscript[ψ, net])>0}]/.Log[(((E^(-(1/(Subscript[ψ, d])))) ((Subscript[ψ, d])+(Subscript[ψ, net])))/(Subscript[ψ, d]))]->-(1/(Subscript[ψ, d]))+Log[(((Subscript[ψ, d])+(Subscript[ψ, net]))/(Subscript[ψ, d]))]

This critical value for ρ can be written as:

ρcrit=(((Subscript[ψ, net]) Log[1/c])/(2 (1+((Subscript[ψ, d])/(Subscript[ψ, net]))) ((Subscript[ψ, d]) Log[1+((Subscript[ψ, net])/(Subscript[ψ, d]))]-1)));
Simplify[%%/ρcrit,{0<c<1,(Subscript[ψ, d])>0,(Subscript[ψ, net])>0}]

The above has the form of  ((Subscript[ψ, net])/(Subscript[ψ, d]))const, where the const term hovers around 1/4 when c = 1/2, for values of interest:

Plot3D[((Log[1/c])/(2 (1+((Subscript[ψ, d])/(Subscript[ψ, net]))) ( Log[1+((Subscript[ψ, net])/(Subscript[ψ, d]))]-1/(Subscript[ψ, d]))))/.c->1/2,{(Subscript[ψ, net]),10,100},{(Subscript[ψ, d]),10,100},AxesLabel->{(Subscript[ψ, net]),(Subscript[ψ, d]),C}]

As we have assumed that selection is not too small, the above is reasonably well approximated when (Subscript[ψ, d]) > 1 by:

ρcrit1=((Subscript[ψ, net])/(Subscript[ψ, d]))(( Log[1/c])/(2 (1+((Subscript[ψ, d])/(Subscript[ψ, net])))  Log[1+((Subscript[ψ, net])/(Subscript[ψ, d]))]));

which depends only on the ratio, ((Subscript[ψ, net])/(Subscript[ψ, d])).  The above has the form of  ((Subscript[ψ, net])/(Subscript[ψ, d]))const, where the const term on the left is ~0.35 when c = 1/2, unless (Subscript[ψ, d]) is small relative to (Subscript[ψ, net]):

Plot[ (( Log[1/c])/(2 (1+((Subscript[ψ, d])/(Subscript[ψ, net])))  Log[1+((Subscript[ψ, net])/(Subscript[ψ, d]))]))/.c->1/2/.(Subscript[ψ, d])->r*(Subscript[ψ, net]),{r,0.001,5},PlotRange->{Automatic,{0,1}}]

Figure 1: Contour plots

<<PlotLegends`
Needs["ErrorBarPlots`"];

ρcrit=(((Subscript[ψ, net]) Log[1/c])/(2 (1+((Subscript[ψ, d])/(Subscript[ψ, net]))) ((Subscript[ψ, d]) Log[1+((Subscript[ψ, net])/(Subscript[ψ, d]))]-1)));

ρ vs (Subscript[ψ, d])

Here we compare the results to a contour plot where the net selection rate is set to (Subscript[ψ, net]) = N ((Subscript[s, a])-(Subscript[s, d])) = 20.

tabρ1=Table[Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[(Subscript[ψ, net])+10^(-10),(Subscript[ψ, d]),ρ]/.p->1/10000/.(Subscript[ψ, net])->20,50]]],{ρ,1/200,3/2,1/20},{(Subscript[ψ, d]),2,100,3}]

The probability that the system is not killed (and so 11 fixes) is given in purple.
	Values of (Subscript[ψ, d]) are arrayed along the x-axis.  
	Values of ρ are arrayed along the y-axis.

ShowLegend[ListContourPlot[(1-tabρ1),DataRange->{{2.1,100},{1/200,3/2}},ColorFunctionScaling->False],
{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

As (Subscript[ψ, d]) increases (x-axis), the 10 haplotype becomes increasingly favored whenever it appears by recombination, making it less likely for the 11 haplotype to sweep to fixation.

Note that the approximation gives a very similar figure (discrepancies arise, however, if ψ values less than 2 are used):

tab1ρb=Table[((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, d])+(Subscript[ψ, net]))/(Subscript[ψ, d]))))^(-ω))/.ω->((2 ρ (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d])))/((Subscript[ψ, net])^2))/.(Subscript[ψ, net])->20,{ρ,1/200,3/2,1/20},{(Subscript[ψ, d]),2,100,3}]//N;

ShowLegend[ListContourPlot[(1-tab1ρb),DataRange->{{2.1,100},{1/200,3/2}},ColorFunctionScaling->False],{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

Furthermore, the critical value ρcrit lines up with the correct 10%, 50%, and 90% contour lines:

plotA=Show[ListContourPlot[(1-tabρ1),DataRange->{{2.1,100},{1/200,3/2}},ColorFunctionScaling->False],
ContourPlot[ρcrit-ρ/.c->{1/10,1/2,9/10}/.(Subscript[ψ, net])->20,{(Subscript[ψ, d]),2.1,100},{ρ,1/200,3/2},Contours->{{0,{Thick,Blue,Dashing[{0.03,0.04}]}}},ContourShading->None]]

Here we compare the results to a contour plot where the net selection is set to (Subscript[ψ, net]) = N ((Subscript[s, a])-(Subscript[s, d])) = 50.

tabρ2=Table[Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[(Subscript[ψ, net])+10^(-10),(Subscript[ψ, d]),ρ]/.p->1/10000/.(Subscript[ψ, net])->50,50]]],{ρ,1/20,3/2,1/20},{(Subscript[ψ, d]),2,100,3}]

The probability that the system is not killed (and so 11 fixes) is given in purple.
	Values of (Subscript[ψ, d]) are arrayed along the x-axis.  
	Values of ρ are arrayed along the y-axis.

ShowLegend[ListContourPlot[(1-tabρ2),DataRange->{{2.1,100},{1/20,3/2}},ColorFunctionScaling->False],{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

Note that the approximation gives a very similar figure (discrepancies arise, however, if ψ values less than 2 are used):

tabρ2b=Table[((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, d])+(Subscript[ψ, net]))/(Subscript[ψ, d]))))^(-ω))/.ω->((2 ρ (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d])))/((Subscript[ψ, net])^2))/.(Subscript[ψ, net])->50,{ρ,1/20,3/2,1/20},{(Subscript[ψ, d]),2,100,3}]//N;

ShowLegend[ListContourPlot[(1-tabρ2b),DataRange->{{2.1,100},{1/20,3/2}},ColorFunctionScaling->False],{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

Furthermore, the critical value ρcrit lines up with the correct 10%, 50%, and 90% contour lines:

plotB=Show[ListContourPlot[(1-tabρ2),DataRange->{{2.1,100},{1/20,3/2}},ColorFunctionScaling->False],
ContourPlot[ρcrit-ρ/.c->{1/10,1/2,9/10}/.(Subscript[ψ, net])->50,{(Subscript[ψ, d]),2.1,100},{ρ,1/20,3/2},Contours->{{0,{Thick,Blue,Dashing[{0.03,0.04}]}}},ContourShading->None]]

ρcritrough=1/4(Subscript[ψ, net])/(Subscript[ψ, d]);

plotA=Show[ListContourPlot[(1-tabρ1),DataRange->{{2.1,100},{1/200,3/2}},ColorFunctionScaling->False],
ContourPlot[ρcrit-ρ/.c->{1/10,1/2,9/10}/.(Subscript[ψ, net])->20,{(Subscript[ψ, d]),2.1,100},{ρ,1/200,3/2},Contours->{{0,{Thick,Blue,Dashing[{0.03,0.04}]}}},ContourShading->None],
ContourPlot[ρcritrough-ρ/.c->{1/10,1/2,9/10}/.(Subscript[ψ, net])->20,{(Subscript[ψ, d]),2.1,100},{ρ,1/200,3/2},Contours->{{0,{Thick,Red,Dashing[{0.03,0.04}]}}},ContourShading->None]]

plotB=Show[ListContourPlot[(1-tabρ2),DataRange->{{2.1,100},{1/20,3/2}},ColorFunctionScaling->False],
ContourPlot[ρcrit-ρ/.c->{1/10,1/2,9/10}/.(Subscript[ψ, net])->50,{(Subscript[ψ, d]),2.1,100},{ρ,1/20,3/2},Contours->{{0,{Thick,Blue,Dashing[{0.03,0.04}]}}},ContourShading->None],
ContourPlot[ρcritrough-ρ/.c->{1/10,1/2,9/10}/.(Subscript[ψ, net])->50,{(Subscript[ψ, d]),2.1,100},{ρ,1/20,3/2},Contours->{{0,{Thick,Red,Dashing[{0.03,0.04}]}}},ContourShading->None]]

ρ vs (Subscript[ψ, net])

Here we compare the results to a contour plot where the deleterious selection is set to (Subscript[ψ, d]) = N (Subscript[s, d]) = 20.

tabρ3=Table[Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[(Subscript[ψ, net])+10^(-10),(Subscript[ψ, d]),ρ]/.p->1/10000/.(Subscript[ψ, d])->20,50]]],{ρ,1/200,3/2,1/20},{(Subscript[ψ, net]),2,100,3}]

The probability that the system is not killed (and so 11 fixes) is given in purple.
	Values of (Subscript[ψ, net]) are arrayed along the x-axis.  
	Values of ρ are arrayed along the y-axis.

ShowLegend[ListContourPlot[(1-tabρ3),DataRange->{{2.1,100},{1/200,3/2}},ColorFunctionScaling->False],
{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

As (Subscript[ψ, d]) increases (x-axis), the 10 haplotype becomes increasingly favored whenever it appears by recombination, making it less likely for the 11 haplotype to sweep to fixation.

As (Subscript[ψ, net]) increases (y-axis), by contrast, it becomes easier for the 11 haplotype to fix because the sweeps are faster.

Note that the approximation gives a very similar figure (discrepancies arise, however, if ψ values less than 2 are used):

tab3ρb=Table[((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, d])+(Subscript[ψ, net]))/(Subscript[ψ, d]))))^(-ω))/.ω->((2 ρ (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d])))/((Subscript[ψ, net])^2))/.(Subscript[ψ, d])->20,{ρ,1/200,3/2,1/20},{(Subscript[ψ, net]),2,100,3}]//N;

ShowLegend[ListContourPlot[(1-tab3ρb),DataRange->{{2.1,100},{1/200,3/2}},ColorFunctionScaling->False],{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

Furthermore, the critical value ρcrit lines up with the correct 10%, 50%, and 90% contour lines:

plotC=Show[ListContourPlot[(1-tabρ3),DataRange->{{2.1,100},{1/200,3/2}},ColorFunctionScaling->False],
ContourPlot[ρcrit-ρ/.c->{1/10,1/2,9/10}/.(Subscript[ψ, d])->20,{(Subscript[ψ, net]),2.1,100},{ρ,1/200,3/2},Contours->{{0,{Thick,Blue,Dashing[{0.03,0.04}]}}},ContourShading->None]]

Here we compare the results to a contour plot where the net selection is set to (Subscript[ψ, d]) = N (Subscript[s, d]) = 50.

tabρ4=Table[Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[(Subscript[ψ, net])+10^(-10),(Subscript[ψ, d]),ρ]/.p->1/10000/.(Subscript[ψ, d])->50,50]]],{ρ,1/20,3/2,1/20},{(Subscript[ψ, net]),2,100,3}]

The probability that the system is not killed (and so 11 fixes) is given in purple.
	Values of (Subscript[ψ, net]) are arrayed along the x-axis.  
	Values of ρ are arrayed along the y-axis.

ShowLegend[ListContourPlot[(1-tabρ4),DataRange->{{2.1,100},{1/20,3/2}},ColorFunctionScaling->False],{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

Note that the approximation gives a very similar figure (discrepancies arise, however, if ψ values less than 2 are used):

tabρ4b=Table[((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, d])+(Subscript[ψ, net]))/(Subscript[ψ, d]))))^(-ω))/.ω->((2 ρ (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d])))/((Subscript[ψ, net])^2))/.(Subscript[ψ, d])->50,{ρ,1/20,3/2,1/20},{(Subscript[ψ, net]),2,100,3}]//N;

ShowLegend[ListContourPlot[(1-tabρ4b),DataRange->{{2.1,100},{1/20,3/2}},ColorFunctionScaling->False],{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

Furthermore, the critical value ρcrit lines up with the correct 10%, 50%, and 90% contour lines:

plotD=Show[ListContourPlot[(1-tabρ4),DataRange->{{2.1,100},{1/20,3/2}},ColorFunctionScaling->False],
ContourPlot[ρcrit-ρ/.c->{1/10,1/2,9/10}/.(Subscript[ψ, d])->50,{(Subscript[ψ, net]),2.1,100},{ρ,1/20,3/2},Contours->{{0,{Thick,Blue,Dashing[{0.03,0.04}]}}},ContourShading->None]]

Figure 1

The following plots ρ (y-axis) versus (Subscript[ψ, d]) (left) with (Subscript[ψ, net]) = 20 (top) or 50 (bottom) and versus (Subscript[ψ, net])(right) with (Subscript[ψ, d]) = 20 (top) or 50 (bottom).  The dashed curve is from ρcrit and it lies virtually on top of the 10%, 50%, and 90% contours (discrepancies arise if  (Subscript[ψ, d]) or (Subscript[ψ, net]) fall below 2).

GraphicsGrid[{{plotA,plotC},{plotB,plotD}}]

Export["contours.eps",%];

(Subscript[ψ, net]) vs (Subscript[ψ, d])

Here we compare the results to a contour plot where the recombination rate is set to ρ = N r = 1/10.

The table is arrayed such that the first row gives the values with ψ = (Subscript[ψ, net]):

tab1=Table[Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[(Subscript[ψ, net])+10^(-10),(Subscript[ψ, d]),ρ]/.p->1/10000/.ρ->1/10,50]]],{(Subscript[ψ, net]),2,100,3},{(Subscript[ψ, d]),2,100,3}]

The probability that the system is not killed (and so 11 fixes) is given in purple.
	Values of (Subscript[ψ, d]) are arrayed along the x-axis.  
	Values of (Subscript[ψ, net]) are arrayed along the y-axis.

ShowLegend[ListContourPlot[(1-tab1),DataRange->{{2.1,100},{2.1,100}},ColorFunctionScaling->False],
{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

As (Subscript[ψ, d]) increases (x-axis), the 10 haplotype becomes increasingly favored whenever it appears by recombination, making it less likely for the 11 haplotype to sweep to fixation.

As (Subscript[ψ, net]) increases (y-axis), by contrast, it becomes easier for the 11 haplotype to fix because the sweeps are faster.

Note that the approximation gives a very similar figure (discrepancies arise, however, if ψ values less than 2 are used):

tab1b=Table[((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, d])+(Subscript[ψ, net]))/(Subscript[ψ, d]))))^(-ω))/.ω->((2 ρ (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d])))/((Subscript[ψ, net])^2))/.ρ->1/10,{(Subscript[ψ, net]),2,100,3},{(Subscript[ψ, d]),2,100,3}]//N;

ShowLegend[ListContourPlot[(1-tab1b),DataRange->{{2.1,100},{2.1,100}},ColorFunctionScaling->False],{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

Furthermore, the critical value ρcrit lines up with the correct 50% contour line:

Show[ListContourPlot[(1-tab1),DataRange->{{2.1,100},{2.1,100}},ColorFunctionScaling->False],
ContourPlot[ρcrit-1/10/.c->1/2,{(Subscript[ψ, d]),2,100},{(Subscript[ψ, net]),2,100},Contours->{{0,{Thick,Blue,Dashing[{0.03,0.04}]}}},ContourShading->None]]

Here we compare the results to a contour plot where the recombination rate is set to ρ = N r = 1.

tab2=Table[Re[Block[{$MaxExtraPrecision=1200},N[notkilledC[(Subscript[ψ, net])+10^(-10),(Subscript[ψ, d]),ρ]/.p->1/10000/.ρ->1,50]]],{(Subscript[ψ, net]),2,100,3},{(Subscript[ψ, d]),2,100,3}]

The probability that the system is not killed (and so 11 fixes) is given in purple.
	Values of (Subscript[ψ, d]) are arrayed along the x-axis.  
	Values of (Subscript[ψ, net]) are arrayed along the y-axis.

tab2=tab2/.Indeterminate->0;

ShowLegend[ListContourPlot[(1-tab2),DataRange->{{2.1,100},{2.1,100}},ColorFunctionScaling->False],{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

Note that the approximation gives a very similar figure (discrepancies arise, however, if ψ values less than 2 are used):

tab2b=Table[((((E^(-(1/(Subscript[ψ, d]))))(((Subscript[ψ, d])+(Subscript[ψ, net]))/(Subscript[ψ, d]))))^(-ω))/.ω->((2 ρ (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d])))/((Subscript[ψ, net])^2))/.ρ->1,{(Subscript[ψ, net]),2,100,3},{(Subscript[ψ, d]),2,100,3}]//N;

ShowLegend[ListContourPlot[(1-tab2b),DataRange->{{2.1,100},{2.1,100}},ColorFunctionScaling->False],{ColorData["LakeColors"][#1]&,10,"1","0",LegendPosition->{1.1,-.4}}]

Furthermore, the critical value ρcrit lines up with the correct 50% contour line:

Show[ListContourPlot[(1-tab2),DataRange->{{2.1,100},{2.1,100}},ColorFunctionScaling->False],
ContourPlot[ρcrit-1/.c->1/2,{(Subscript[ψ, d]),2,100},{(Subscript[ψ, net]),2,100},Contours->{{0,{Thick,Blue,Dashing[{0.03,0.04}]}}},ContourShading->None]]

Figure 3: Comparison against simulations holding (Subscript[ψ, d]) constant

N=10000 with (Subscript[ψ, d])=10

Simulations with N=10000:
The following gives the proportion fixed for #11, #10, #00 (each immediately followed by 1.96 SE), here for 100,000 replicates.
The nine rows report (Subscript[s, a]) equal to 0.002, then 0.005, then 0.01 for ρ = 0.5, then repeats these values of (Subscript[s, a])for ρ =1 and then ρ =2.

raw10000={{0.0006700000,0.0001603792,0.0015500000,0.0002438288,0.9977800000,0.0002917091},{0.0051400000,0.0004432193,0.0031700000,0.0003484145,0.9916900000,0.0005626581},{0.0134100000,0.0007129167,0.0043500000,0.0004079003,0.9822400000,0.0008186276},{0.0002100000,0.0000898091,0.0022100000,0.0002910528,0.9975800000,0.0003045353},{0.0028500000,0.0003304142,0.0052700000,0.0004487598,0.9918800000,0.0005562418},{0.0108600000,0.0006423916,0.0077800000,0.0005445652,0.9813600000,0.0008382879},{0.0000600000,0.0000480086,0.0023900000,0.0003026463,0.9975500000,0.0003064125},{0.0012800000,0.0002216067,0.0065900000,0.0005014909,0.9921300000,0.0005476811},{0.0055300000,0.0004596365,0.0121900000,0.0006801342,0.9822800000,0.0008177218}};

The probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3):

means=Transpose[raw10000][[1]]/(Transpose[raw10000][[1]]+Transpose[raw10000][[3]])

fixed=(Transpose[raw10000][[1]]+Transpose[raw10000][[3]])*reps/.reps->100000//Round

Note that these fixation probabilities are all slightly too low relative to Kimura's fixation probability for N = 10000, as expected because of the hitchhiking deleterious allele:

(1-Exp[-2s])/(1-Exp[-2N s])*reps/.reps->100000/.N->10000/.s->{0.002,0.005,0.01}//Round

Below, we'll plot the probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3), so we need to adjust the confidence intervals (based on +/- 1.96 SE):

CIs=1.96*Sqrt[means*(1-means)/fixed]

scaled10sd10000=Transpose[{means,CIs}]

Graphs:

N1 = 10000;
initialp=(1/N1);

plot1=Show[ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/1000+10^(-10)),N1/1000,N1/20000]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/1000+10^(-10)),N1/1000,N1/20000]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[1,0,0],Thickness[0.005]}],
ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/1000+10^(-10)),N1/1000,N1/10000]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/1000+10^(-10)),N1/1000,N1/10000]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[0,0,1],Thickness[0.005]}],
ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/1000+10^(-10)),N1/1000,N1/5000]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/1000+10^(-10)),N1/1000,N1/5000]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[0,1,0],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/20000/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->10,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[1,0,0],Dashing[0.016],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/10000/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->10,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,0,1],Dashing[0.016],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/5000/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->10,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,1,0],Dashing[0.016],Thickness[0.005]}]
];

The following plots the semi-deterministic solution:

(Subscript[trys, d])=10/N1;
tryr1=1/(2N1);
tryr2=1/N1;
tryr3=2/N1;
plot1b=Show[
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr1/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[1,0,0],Thickness[0.001]}],
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr2/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[0,0,1],Thickness[0.001]}],
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr3/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[0,1,0],Thickness[0.001]}]

]

plot2=ErrorListPlot[{{{20,scaled10sd10000[[1,1]]},ErrorBar[scaled10sd10000[[1,2]]]},{{50,scaled10sd10000[[2,1]]},ErrorBar[scaled10sd10000[[2,2]]]},{{100,scaled10sd10000[[3,1]]},ErrorBar[scaled10sd10000[[3,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[1,0,0],Thickness[0.005],PointSize[0.015]}];
plot2b=ErrorListPlot[{{{20,scaled10sd10000[[4,1]]},ErrorBar[scaled10sd10000[[4,2]]]},{{50,scaled10sd10000[[5,1]]},ErrorBar[scaled10sd10000[[5,2]]]},{{100,scaled10sd10000[[6,1]]},ErrorBar[scaled10sd10000[[6,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,0,1],Thickness[0.005],PointSize[0.015]}];
plot2c=ErrorListPlot[{{{20,scaled10sd10000[[7,1]]},ErrorBar[scaled10sd10000[[7,2]]]},{{50,scaled10sd10000[[8,1]]},ErrorBar[scaled10sd10000[[8,2]]]},{{100,scaled10sd10000[[9,1]]},ErrorBar[scaled10sd10000[[9,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,1,0],Thickness[0.005],PointSize[0.015]}];


largeN=Show[plot1,plot1b,plot2,plot2b,plot2c,AxesLabel->{(Subscript[ψ, a]),(SuperStar[P])},PlotRange->{{0,120},{-0.01,1}},AxesOrigin->{0,0}]

Export["det10000.eps",%];

N=1000 with (Subscript[ψ, d])=10

Simulations with N=1000:
The following gives the proportion fixed for #11, #10, #00 (each immediately followed by 1.96 SE), here for 100,000 replicates.
The nine rows report (Subscript[s, a]) equal to 0.02, then 0.05, then 0.1 for ρ = 0.5, then repeats these values of (Subscript[s, a])for ρ =1 and then ρ =2.

raw1000={{0.0064300000,0.0004954055,0.0145000000,0.0007409150,0.9790700000,0.0008872532},{0.0472600000,0.0013151949,0.0282700000,0.0010272870,0.9244700000,0.0016378051},{0.1256000000,0.0020540261,0.0377500000,0.0011812954,0.8366500000,0.0022913295},{0.0021900000,0.0002897357,0.0186000000,0.0008374050,0.9792100000,0.0008843440},{0.0291900000,0.0010433746,0.0469300000,0.0013108221,0.9238800000,0.0016436648},{0.0953100000,0.0018200160,0.0650700000,0.0015287484,0.8396200000,0.0022744299},{0.0004500000,0.0001314512,0.0213100000,0.0008950976,0.9782400000,0.0009042910},{0.0127700000,0.0006959222,0.0639400000,0.0015163318,0.9232900000,0.0016494955},{0.0563200000,0.0014288932,0.1020100000,0.0018759152,0.8416700000,0.0022626043}};

The probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3):

means=Transpose[raw1000][[1]]/(Transpose[raw1000][[1]]+Transpose[raw1000][[3]])

The probability that (Subscript[A, 1]) fixes (columns 1 + 3):

fixed=(Transpose[raw1000][[1]]+Transpose[raw1000][[3]])*reps/.reps->100000//Round

Note that these fixation probabilities are all slightly too low relative to Kimura's fixation probability for N = 1000, as expected because of the hitchhiking deleterious allele:

(1-Exp[-2s])/(1-Exp[-2N s])*reps/.reps->100000/.N->1000/.s->{0.02,0.05,0.1}//Round

Below, we'll plot the probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3), so we need to adjust the confidence intervals (based on +/- 1.96 SE):

CIs=1.96*Sqrt[means*(1-means)/fixed]

scaled10sd1000=Transpose[{means,CIs}]

Graphs:

N1 = 1000;
initialp=(1/N1);

Solid curve based on full solution to the diffusion, dashed curve based on ((((e^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω)), which assumes selection reasonably strong, and thin curve based on semi-determistic model.

plot1=Show[ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/100+10^(-10)),N1/100,N1/2000]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/100+10^(-10)),N1/100,N1/2000]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[1,0,0],Thickness[0.005]}],
ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/100+10^(-10)),N1/100,N1/1000]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/100+10^(-10)),N1/100,N1/1000]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[0,0,1],Thickness[0.005]}],
ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/100+10^(-10)),N1/100,N1/500]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/100+10^(-10)),N1/100,N1/500]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[0,1,0],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/2000/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->10,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[1,0,0],Dashing[0.016],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/1000/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->10,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,0,1],Dashing[0.016],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/500/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->10,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,1,0],Dashing[0.016],Thickness[0.005]}]
];

The following plots the semi-deterministic solution:

(Subscript[trys, d])=10/N1;
tryr1=1/(2N1);
tryr2=1/N1;
tryr3=2/N1;
plot1b=Show[
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr1/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[1,0,0],Thickness[0.001]}],
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr2/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[0,0,1],Thickness[0.001]}],
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr3/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[0,1,0],Thickness[0.001]}]

]

plot2=ErrorListPlot[{{{20,scaled10sd1000[[1,1]]},ErrorBar[scaled10sd1000[[1,2]]]},{{50,scaled10sd1000[[2,1]]},ErrorBar[scaled10sd1000[[2,2]]]},{{100,scaled10sd1000[[3,1]]},ErrorBar[scaled10sd1000[[3,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[1,0,0],Thickness[0.005],PointSize[0.015]}];
plot2b=ErrorListPlot[{{{20,scaled10sd1000[[4,1]]},ErrorBar[scaled10sd1000[[4,2]]]},{{50,scaled10sd1000[[5,1]]},ErrorBar[scaled10sd1000[[5,2]]]},{{100,scaled10sd1000[[6,1]]},ErrorBar[scaled10sd1000[[6,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,0,1],Thickness[0.005],PointSize[0.015]}];
plot2c=ErrorListPlot[{{{20,scaled10sd1000[[7,1]]},ErrorBar[scaled10sd1000[[7,2]]]},{{50,scaled10sd1000[[8,1]]},ErrorBar[scaled10sd1000[[8,2]]]},{{100,scaled10sd1000[[9,1]]},ErrorBar[scaled10sd1000[[9,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,1,0],Thickness[0.005],PointSize[0.015]}];

smallN=Show[plot1,plot1b,plot2,plot2b,plot2c,AxesLabel->{(Subscript[ψ, a]),(SuperStar[P])},PlotRange->{{0,120},{-0.01,1}},AxesOrigin->{0,0}]

Export["det1000.eps",%];

N=10000 with (Subscript[ψ, d])=2

Simulations with N=10000:
The following gives the proportion fixed for #11, #10, #00 (each immediately followed by 1.96 SE), here for 100,000 replicates.
The nine rows report (Subscript[s, a]) equal to 0.002, then 0.005, then 0.01 for ρ = 0.5, then repeats these values of (Subscript[s, a])for ρ =1 and then ρ =2.

raw10000={{0.0029400000,0.0003355756,0.0008400000,0.0001795615,0.9962200000,0.0003803467},{0.0082100000,0.0005592906,0.0014600000,0.0002366548,0.9903300000,0.0006065398},{0.0182400000,0.0008294136,0.0013700000,0.0002292549,0.9803900000,0.0008593979},{0.0021900000,0.0002897357,0.0016300000,0.0002500320,0.9961800000,0.0003823461},{0.0070800000,0.0005196726,0.0023400000,0.0002994713,0.9905800000,0.0005987235},{0.0160500000,0.0007788971,0.0027700000,0.0003257569,0.9811800000,0.0008422485},{0.0011800000,0.0002127848,0.0022900000,0.0002962620,0.9965300000,0.0003644735},{0.0060600000,0.0004810303,0.0037400000,0.0003783365,0.9902000000,0.0006105631},{0.0143200000,0.0007363691,0.0057400000,0.0004682330,0.9799400000,0.0008690030}};

The probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3):

means=Transpose[raw10000][[1]]/(Transpose[raw10000][[1]]+Transpose[raw10000][[3]])

fixed=(Transpose[raw10000][[1]]+Transpose[raw10000][[3]])*reps/.reps->100000//Round

Note that these fixation probabilities are all slightly too low relative to Kimura's fixation probability for N = 10000, as expected because of the hitchhiking deleterious allele:

(1-Exp[-2s])/(1-Exp[-2N s])*reps/.reps->100000/.N->10000/.s->{0.002,0.005,0.01}//Round

Below, we'll plot the probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3), so we need to adjust the confidence intervals (based on +/- 1.96 SE):

CIs=1.96*Sqrt[means*(1-means)/fixed]

scaled10sd10000=Transpose[{means,CIs}]

Graphs:

N1 = 10000;
initialp=(1/N1);

plot1=Show[ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/5000+10^(-10)),N1/5000,N1/20000]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/5000+10^(-10)),N1/5000,N1/20000]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[1,0,0],Thickness[0.005]}],
ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/5000+10^(-10)),N1/5000,N1/10000]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/5000+10^(-10)),N1/5000,N1/10000]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[0,0,1],Thickness[0.005]}],
ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/5000+10^(-10)),N1/5000,N1/5000]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/5000+10^(-10)),N1/5000,N1/5000]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[0,1,0],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/20000/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->2,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[1,0,0],Dashing[0.016],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/10000/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->2,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,0,1],Dashing[0.016],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/5000/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->2,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,1,0],Dashing[0.016],Thickness[0.005]}]
];

The following plots the semi-deterministic solution:

(Subscript[trys, d])=2/N1;
tryr1=1/(2N1);
tryr2=1/N1;
tryr3=2/N1;
plot1b=Show[
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr1/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[1,0,0],Thickness[0.001]}],
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr2/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[0,0,1],Thickness[0.001]}],
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr3/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[0,1,0],Thickness[0.001]}]

]

plot2=ErrorListPlot[{{{20,scaled10sd10000[[1,1]]},ErrorBar[scaled10sd10000[[1,2]]]},{{50,scaled10sd10000[[2,1]]},ErrorBar[scaled10sd10000[[2,2]]]},{{100,scaled10sd10000[[3,1]]},ErrorBar[scaled10sd10000[[3,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[1,0,0],Thickness[0.005],PointSize[0.015]}];
plot2b=ErrorListPlot[{{{20,scaled10sd10000[[4,1]]},ErrorBar[scaled10sd10000[[4,2]]]},{{50,scaled10sd10000[[5,1]]},ErrorBar[scaled10sd10000[[5,2]]]},{{100,scaled10sd10000[[6,1]]},ErrorBar[scaled10sd10000[[6,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,0,1],Thickness[0.005],PointSize[0.015]}];
plot2c=ErrorListPlot[{{{20,scaled10sd10000[[7,1]]},ErrorBar[scaled10sd10000[[7,2]]]},{{50,scaled10sd10000[[8,1]]},ErrorBar[scaled10sd10000[[8,2]]]},{{100,scaled10sd10000[[9,1]]},ErrorBar[scaled10sd10000[[9,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,1,0],Thickness[0.005],PointSize[0.015]}];


largeN=Show[plot1,plot1b,plot2,plot2b,plot2c,AxesLabel->{(Subscript[ψ, a]),(SuperStar[P])},PlotRange->{{0,120},{-0.01,1}},AxesOrigin->{0,0}]

Export["det10000b.eps",%];

N=1000 with (Subscript[ψ, d])=2

Simulations with N=1000:
The following gives the proportion fixed for #11, #10, #00 (each immediately followed by 1.96 SE), here for 100,000 replicates.
The nine rows report (Subscript[s, a]) equal to 0.02, then 0.05, then 0.1 for ρ = 0.5, then repeats these values of (Subscript[s, a])for ρ =1 and then ρ =2.

raw1000={{0.0263700000,0.0009931346,0.0084200000,0.0005663384,0.9652100000,0.0011357799},{0.0789100000,0.0016709872,0.0114600000,0.0006596985,0.9096300000,0.0017770539},{0.1606400000,0.0022759203,0.0128800000,0.0006988742,0.8264800000,0.0023471833},{0.0207100000,0.0008826769,0.0150600000,0.0007548722,0.9642300000,0.0011510809},{0.0698100000,0.0015794312,0.0210900000,0.0008905652,0.9091000000,0.0017817380},{0.1491700000,0.0022080973,0.0247000000,0.0009619969,0.8261300000,0.0023490518},{0.0132700000,0.0007092359,0.0231100000,0.0009312770,0.9636200000,0.0011604871},{0.0545900000,0.0014080650,0.0354200000,0.0011456434,0.9099900000,0.0017738617},{0.1268800000,0.0020629544,0.0443100000,0.0012754560,0.8288100000,0.0023346552}};

The probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3):

means=Transpose[raw1000][[1]]/(Transpose[raw1000][[1]]+Transpose[raw1000][[3]])

The probability that (Subscript[A, 1]) fixes (columns 1 + 3):

fixed=(Transpose[raw1000][[1]]+Transpose[raw1000][[3]])*reps/.reps->100000//Round

Note that these fixation probabilities are all slightly too low relative to Kimura's fixation probability for N = 1000, as expected because of the hitchhiking deleterious allele:

(1-Exp[-2s])/(1-Exp[-2N s])*reps/.reps->100000/.N->1000/.s->{0.02,0.05,0.1}//Round

Below, we'll plot the probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3), so we need to adjust the confidence intervals (based on +/- 1.96 SE):

CIs=1.96*Sqrt[means*(1-means)/fixed]

scaled10sd1000=Transpose[{means,CIs}]

Graphs:

N1 = 1000;
initialp=(1/N1);

Solid curve based on full solution to the diffusion, dashed curve based on ((((e^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω)), which assumes selection reasonably strong, and thin curve based on semi-determistic model.

plot1=Show[ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/500+10^(-10)),N1/500,N1/2000]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/500+10^(-10)),N1/500,N1/2000]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[1,0,0],Thickness[0.005]}],
ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/500+10^(-10)),N1/500,N1/1000]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/500+10^(-10)),N1/500,N1/1000]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[0,0,1],Thickness[0.005]}],
ListPlot[Join[Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/500+10^(-10)),N1/500,N1/500]/.p->1/N1,50]]]},{Nsa,11,19,1}],Table[{Nsa,Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[N1*(Nsa/N1-1/500+10^(-10)),N1/500,N1/500]/.p->1/N1,50]]]},{Nsa,20,120,10}]],Joined->True,AxesOrigin->{0,0},PlotRange->{Automatic,{0,1}},PlotStyle->{RGBColor[0,1,0],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/2000/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->2,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[1,0,0],Dashing[0.016],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/1000/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->2,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,0,1],Dashing[0.016],Thickness[0.005]}],
Plot[((((E^(-(1/(Subscript[ψ, d]))))((Subscript[ψ, a])/(Subscript[ψ, d]))))^(-ω))/.ω->ρ/(( ((Subscript[ψ, net])^2))/(2 (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d]))))/.ρ->N1/500/. (Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->2,{(Subscript[ψ, a]),11.5,120},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,1,0],Dashing[0.016],Thickness[0.005]}]
];

The following plots the semi-deterministic solution:

(Subscript[trys, d])=2/N1;
tryr1=1/(2N1);
tryr2=1/N1;
tryr3=2/N1;
plot1b=Show[
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr1/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[1,0,0],Thickness[0.001]}],
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr2/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[0,0,1],Thickness[0.001]}],
Plot[Exp[-((2 N r sa sd Log[(sa/sd)])/(((sa-sd))^2))]/.N->N1/.r->tryr3/.sd->(Subscript[trys, d])/.sa->Nsa/N1,{Nsa,10,120},PlotStyle->{RGBColor[0,1,0],Thickness[0.001]}]

]

plot2=ErrorListPlot[{{{20,scaled10sd1000[[1,1]]},ErrorBar[scaled10sd1000[[1,2]]]},{{50,scaled10sd1000[[2,1]]},ErrorBar[scaled10sd1000[[2,2]]]},{{100,scaled10sd1000[[3,1]]},ErrorBar[scaled10sd1000[[3,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[1,0,0],Thickness[0.005],PointSize[0.015]}];
plot2b=ErrorListPlot[{{{20,scaled10sd1000[[4,1]]},ErrorBar[scaled10sd1000[[4,2]]]},{{50,scaled10sd1000[[5,1]]},ErrorBar[scaled10sd1000[[5,2]]]},{{100,scaled10sd1000[[6,1]]},ErrorBar[scaled10sd1000[[6,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,0,1],Thickness[0.005],PointSize[0.015]}];
plot2c=ErrorListPlot[{{{20,scaled10sd1000[[7,1]]},ErrorBar[scaled10sd1000[[7,2]]]},{{50,scaled10sd1000[[8,1]]},ErrorBar[scaled10sd1000[[8,2]]]},{{100,scaled10sd1000[[9,1]]},ErrorBar[scaled10sd1000[[9,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,1,0],Thickness[0.005],PointSize[0.015]}];

smallN=Show[plot1,plot1b,plot2,plot2b,plot2c,AxesLabel->{(Subscript[ψ, a]),(SuperStar[P])},PlotRange->{{0,120},{-0.01,1}},AxesOrigin->{0,0}]

Export["det1000b.eps",%];

Additional figures

In these additional simulations, we explored how rapidly the solutions derived in the paper break down when either selection becomes stronger (first two graphs) or the population size becomes smaller (second two graphs).

Simulations with N=10000 and stronger selection ((Subscript[ψ, d])=100, that is, (Subscript[s, d])=0.01):
The figure reports results with (Subscript[s, a]) equal to 0.011,  0.014, and 0.019 (dots along x-axis) for ρ = 0.5 (red), ρ =1 (blue), and ρ =2 (green).

Simulations with N=1000 and stronger selection ((Subscript[ψ, d])=100, that is, (Subscript[s, d])=0.1):
The figure reports results with (Subscript[s, a]) equal to 0.11,  0.14, and 0.19 (dots along x-axis) for ρ = 0.5 (red), ρ =1 (blue), and ρ =2 (green).

This figure demonstrates that the solutions all break down as selection gets too strong (here,  (Subscript[s, a])≥ 0.11 and (Subscript[s, d]) = 0.1)

Simulations with smaller population size N=100 ((Subscript[ψ, d])=10, that is, (Subscript[s, d])=0.1):
The figure reports results with (Subscript[s, a]) equal to 0.2,  0.5, and 1 (dots along x-axis) for ρ = 0.5 (red), ρ =1 (blue), and ρ =2 (green).

This figure demonstrates that the solutions all break down as selection gets too strong (here,  (Subscript[s, a]) ≥ 0.2 and (Subscript[s, d]) = 0.1)

Simulations with smaller population size N=100 ((Subscript[ψ, d])=2, that is, (Subscript[s, d])=0.02):
The figure reports results with (Subscript[s, a]) equal to 0.2,  0.5, and 1 (dots along x-axis) for ρ = 0.5 (red), ρ =1 (blue), and ρ =2 (green).

This figure demonstrates that the solutions still perform reasonably well in small population sizes, but again the solutions perform more poorly as selection becomes stronger (here,  (Subscript[s, a]) ≥ 0.2 and (Subscript[s, d]) = 0.02).

Supplementary material 4:  The diffusion process with killing for a linked neutral allele

Diffusion equations

When the (Subscript[B, 1])allele is neutral rather than deleterious, the above diffusion approximation becomes undefined.

While the drift and diffusion terms remain the same, the killing term must be updated.  Now, when the recombinant 10 haplotype arises, there are p copies of the beneficial (Subscript[A, 1]) allele, and we assume that each of these p copies is equally likely to be the ancestor of the population at some distant future point in time.  Thus, the recombinant 10 haplotype is assumed to have a 1/(N p) fixation probability.  

[Technically, this assumes that only the (Subscript[A, 1]) allele has appreciable chance of fixing and that the probability of (Subscript[A, 0]) fixing is negligible.  This is always appropriate for the conditional process, but it will only be appropriate for the unconditional process if the system starts with enough copies of (Subscript[A, 1]) that its fixation is assured.  Our focus is on the conditional process:  amongst those processes where (Subscript[A, 1])fixes, how often does the (Subscript[A, 1])(Subscript[B, 0]) recombinant appear and establish.]

Accounting for the probability that a recombinant is formed, the probability that any one offspring is a 10 successful recombinant is:

	 κ = r p (1-p) (1/(N p))
	 
and the process will be killed if any one of the offspring is a successful recombinant:

neutralpart3=1-PDF[BinomialDistribution[N, κ ],0]//Simplify

We must assume that κ is sufficiently small that even N*κ is small, in which case, we can replace the above with the approximation:

neutralpart3a=N κ;

Scaling time by (N^α) and replacing for the κ gives, we get the killing term:

neutralpart3b=neutralpart3a*(N^α)/.κ->r p (1-p) (1/(N p))

Scaling the selection and recombination terms as before:

Factor[%/.α->1/.s->S+ψ/N/.S->ψd/N/.r->ρ/N]

This is now placed in the diffusion equation (8), and each term is divided by p(1-p) to get a differential equation whose solution gives P[p], the probability of the process not being killed (ending up fixed for either the wildtype 00 or the doubly mutant 11 haplotype).

DSolve has trouble solving this differential equation directly, so we replace (ρ/p) with (ρ/(p+d)) and later let d go to zero .

P[p]/.Flatten[Simplify[DSolve[{((P''[p])/2)+ψnet*P'[p]-ρ*(1/(p+d))*P[p]==0,P[0]==P[1]==1},P[p],p],{0<p<1,ψnet>0,ρ>0}]]

neutralnotkilled[ψnet_,ρ_]=FullSimplify[%/.d->0,{ρ>0,ψnet>0}]

Similarly, we can adjust the killing term to get the conditional probability that the process is not killed (i.e., the probability that (Subscript[B, 1])fixes given that (Subscript[A, 1]) has fixed):

P[p]/.Flatten[Simplify[DSolve[{((P''[p])/2)+ψnet*(((1+(E^(-2 p ψnet)))/(1-(E^(-2 p ψnet)))))*P'[p]-ρ*(1/(p+d))*P[p]==0,P'[c]==0,P[1]==1},P[p],p],{0<p<1,ψnet>0,ρ>0}]]

neutralnotkilledC[ψnet_,ρ_]=FullSimplify[%/.c->0/.d->0,{ρ>0,ψnet>0}]

Approximation for the probability of hitchhiking (conditional)

As we might expect, the conditional probability of hitchhiking declines roughly exponentially with the recombination rate (through ρ):

plot1=Show[ListPlot[Table[{ρ,Re[Block[{$MaxExtraPrecision=1000},N[Limit[neutralnotkilledC[40,ρ],p->0],50]]]},{ρ,1/10,10,1/4}],Joined->True,PlotRange->{{0,10},{0,1}}],
ListPlot[Table[{ρ,Re[Block[{$MaxExtraPrecision=1000},N[Limit[neutralnotkilledC[20,ρ],p->0],50]]]},{ρ,1/10,10,1/4}],Joined->True,PlotRange->{{0,10},{0,1}}],
ListPlot[Table[{ρ,Re[Block[{$MaxExtraPrecision=1000},N[Limit[neutralnotkilledC[10,ρ],p->0],50]]]},{ρ,1/10,10,1/4}],Joined->True,PlotRange->{{0,10},{0,1}}],
AxesLabel->{ρ,(SuperStar[P])}
]

We assume that the initial frequency of the 11 haplotype is 1/N and that this is very small.  We thus focus on the conditional probability given that the (Subscript[A, 1]) allele ultimately fixes, in the limit when p -> 0.

FullSimplify[Limit[neutralnotkilledC[(Subscript[ψ, net]),ρ],p->0],{(Subscript[ψ, net])>0,ω>0}]

In the following, we assume that selection is strong relative to drift and find the value that the above approaches when (Subscript[ψ, net]) is large.

From Abramowitz and Stegun equation 13.1.4,  Hypergeometric1F1[a, b, z] -> (Gamma[b]/Gamma[a])*Exp[z]*z^(a - b) as z becomes large.

Simplify[%/.Hypergeometric1F1[a_,b_,z_]->(Gamma[b]/Gamma[a])*Exp[z]*z^(a-b)/.z->2 (Subscript[ψ, net])/.b->2/.a->1+(ρ/(Subscript[ψ, net])),{ρ>0,(Subscript[ψ, net])>0,(Subscript[ψ, net])>ρ}]

Given that this is a large ψ approximation, we take (e^(-2 (Subscript[ψ, net]))) (-1+(e^(2 (Subscript[ψ, net])))) to be very nearly one, leaving us with (((2(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net])))) Gamma[1+(ρ/(Subscript[ψ, net]))].

To determine the best exponential approximation for (((2(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net])))) Gamma[1+(ρ/(Subscript[ψ, net]))] as a function of ρ, we take the Taylor series of the log of the above.  This gives us the following approximation for this exponential function:

Exp[Simplify[Normal[Series[Log[(((2(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net])))) Gamma[1+(ρ/(Subscript[ψ, net]))]/.ρ->ρ*small],{small,0,1}]]]/.small->1]

which can be written as:

temp=(((2 (E^EulerGamma)(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net]))));

as reported in equation (14).  

This form compares well to the neutralnotkilledC[] function evaluated above:

Show[plot1,
Plot[{temp/.(Subscript[ψ, net])->10,temp/.(Subscript[ψ, net])->20,temp/.(Subscript[ψ, net])->40},{ρ,0,10},AxesLabel->{ρ,P},PlotStyle->Dashing[{0.01}]]]

Solving for the critical value of recombination above which hitchhiking to fixation becomes unlikely

We can then solve for the critical value of ρ, where the linked neutral allele hitchhikes to fixation 50% of the time, on average:

ρcritN=ρ/.Flatten[Simplify[Solve[temp==c,ρ]]]

The above has the form of  (Subscript[ψ, net])(Subscript[C, n]), where the (Subscript[C, n]) term is roughly 1/4 for (Subscript[ψ, net]) = 5 and declines slowly for larger values of  (Subscript[ψ, net]).

((Log[1/c])/(EulerGamma+Log[2 (Subscript[ψ, net])]))/.c->1/2/.(Subscript[ψ, net])->{5.,50.,500.}

Plot[((Log[1/c])/(EulerGamma+Log[2 (Subscript[ψ, net])]))/.c->1/2,{(Subscript[ψ, net]),2,100},PlotRange->{All,{0,0.5}}]

Figure 2:  Plot of ρcrit

The following plot varies (Subscript[ψ, net]) along the x-axis and shows the critical value of ρ on the y-axis for the neutral case (thick) and for cases with a deleterious hitchhiking allele (dashed lines from top to bottom show (Subscript[ψ, d]) = 1, 2, 5, 10, and 20 (note that discrepancies arise when selection is too weak):

ρcrit=(((Subscript[ψ, net]) Log[1/c])/(2 (1+((Subscript[ψ, d])/(Subscript[ψ, net]))) ((Subscript[ψ, d]) Log[1+((Subscript[ψ, net])/(Subscript[ψ, d]))]-1)));

ρcritN=((Log[1/c] (Subscript[ψ, net]))/(EulerGamma+Log[2 (Subscript[ψ, net])]));

plotE=Show[Plot[{ρcritN/.c->1/2,ρcrit/.(Subscript[ψ, d])->1/.c->1/2,ρcrit/.(Subscript[ψ, d])->2/.c->1/2,ρcrit/.(Subscript[ψ, d])->5/.c->1/2,ρcrit/.(Subscript[ψ, d])->10/.c->1/2,ρcrit/.(Subscript[ψ, d])->20/.c->1/2},{(Subscript[ψ, net]),2,100},PlotStyle->{{Thickness[0.01]},{Dashing[{0.07}]},{Dashing[{0.042}]},{Dashing[{0.02}]},{Dashing[{0.01}]},{Dashing[{0.005}]}},PlotRange->{Automatic,{0,12}}],Graphics[{Text["(Subscript[ψ, d])=0",{105,11.8}],Text["(Subscript[ψ, d])=1",{105,9.4}],Text["(Subscript[ψ, d])=2",{105,4.9}],Text["(Subscript[ψ, d])=5",{105,2.3}],Text["(Subscript[ψ, d])=10",{105,1.3}],Text["(Subscript[ψ, d])=20",{105,0.8}]}],PlotRange->{{0,109},{0,12.1}}]

Export["neutral.eps",%];

The following graph avoids most of the approximations but generates a very similar picture when evaluated numerically for the critical value of ρ (here set to 50%):

solveme[ψd_,ψnet_]:=ρ/.FindRoot[((((1-(E^(-2 ψnet)))(1+(ω/(ψd))))/( (((((ψnet+ψd)/(ψd))))^ω)(1- (E^(-2 ψnet))(((((ψd^2)/(((ψd+ψnet))^2))))^ω) )))/.ω->ρ/(( ((ψnet)^2))/(2 ψd (ψnet+ψd))))==1/2,{ρ,1/1000,1}]

neutralsolveme[ψnet_]:=ρ/.FindRoot[(((-1+(E^(2 ψnet)))/(2 Hypergeometric1F1[1+(ρ/(ψnet)),2,2 ψnet]ψnet)))==1/2,{ρ,1/1000,10}]

plotF=Plot[{neutralsolveme[(Subscript[ψ, net])],solveme[1,(Subscript[ψ, net])],solveme[2,(Subscript[ψ, net])],solveme[5,(Subscript[ψ, net])],solveme[10,(Subscript[ψ, net])],solveme[20,(Subscript[ψ, net])]},{(Subscript[ψ, net]),2,100},PlotStyle->{{Thickness[0.01]},{Dashing[{0.07}]},{Dashing[{0.042}]},{Dashing[{0.02}]},{Dashing[{0.01}]},{Dashing[{0.005}]}},PlotRange->{Automatic,{0,12}}]

Show[plotE,plotF]

Comparison against simulations

Simulations with N=1000:
The following gives the proportion fixed for #11, #10, #00 (each immediately followed by 1.96 SE), here for 100,000 replicates.
The nine rows report (Subscript[s, a]) equal to 0.02, then 0.05, then 0.1 for ρ = 0.5, then repeats these values of (Subscript[s, a])for ρ =1 and then ρ =2.

raw1000={{0.0351600000,0.0011415847,0.0039100000,0.0003868065,0.9609300000,0.0012009465},{0.0885300000,0.0017606478,0.0046100000,0.0004198587,0.9068600000,0.0018013343},{0.1700200000,0.0023283051,0.0054700000,0.0004571500,0.8245100000,0.0023576548},{0.0318400000,0.0010882187,0.0067800000,0.0005086202,0.9613800000,0.0011942899},{0.0857400000,0.0017353324,0.0094200000,0.0005987235,0.9048400000,0.0018187340},{0.1685100000,0.0023200504,0.0102300000,0.0006236789,0.8212600000,0.0023746920},{0.0265400000,0.0009962437,0.0132900000,0.0007097630,0.9601700000,0.0012120912},{0.0754300000,0.0016368091,0.0179900000,0.0008238148,0.9065800000,0.0018037613},{0.1554000000,0.0022454693,0.0189500000,0.0008450964,0.8256500000,0.0023516086}};

The probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3):

means=Transpose[raw1000][[1]]/(Transpose[raw1000][[1]]+Transpose[raw1000][[3]])

The probability that (Subscript[A, 1]) fixes (columns 1 + 3):

fixed=(Transpose[raw1000][[1]]+Transpose[raw1000][[3]])*reps/.reps->100000//Round

These fixation probabilities are now fairly accurate, relative to Kimura's fixation probability for N = 1000, as expected because the hitchhiking allele is neutral:

(1-Exp[-2s])/(1-Exp[-2N s])*reps/.reps->100000/.N->1000/.s->{0.02,0.05,0.1}//Round

Below, we'll plot the probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3), so we need to adjust the confidence intervals (based on +/- 1.96 SE):

CIs=1.96*Sqrt[means*(1-means)/fixed]

scaledn1000=Transpose[{means,CIs}]

Simulations with N=10000:
The following gives the proportion fixed for #11, #10, #00 (each immediately followed by 1.96 SE), here for 100,000 replicates.
The nine rows report (Subscript[s, a]) equal to 0.002, then 0.005, then 0.01 for ρ = 0.5, then repeats these values of (Subscript[s, a])for ρ =1 and then ρ =2.

raw10000={{0.0035300000,0.0003676000,0.0004000000,0.0001239365,0.9960700000,0.0003877906},{0.0097200000,0.0006080905,0.0004400000,0.0001299831,0.9898400000,0.0006215634},{0.0186800000,0.0008391698,0.0004400000,0.0001299831,0.9808800000,0.0008488051},{0.0030600000,0.0003423350,0.0008000000,0.0001752376,0.9961400000,0.0003843350},{0.0083500000,0.0005639992,0.0009300000,0.0001889276,0.9907200000,0.0005942997},{0.0185400000,0.0008360789,0.0011800000,0.0002127848,0.9802800000,0.0008617566},{0.0026700000,0.0003198388,0.0011500000,0.0002100657,0.9961800000,0.0003823461},{0.0080900000,0.0005552217,0.0018000000,0.0002627248,0.9901100000,0.0006133324},{0.0177100000,0.0008174952,0.0021300000,0.0002857478,0.9801600000,0.0008643217}};

The probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3):

means=Transpose[raw10000][[1]]/(Transpose[raw10000][[1]]+Transpose[raw10000][[3]])

The probability that (Subscript[A, 1]) fixes (columns 1 + 3):

fixed=(Transpose[raw10000][[1]]+Transpose[raw10000][[3]])*reps/.reps->100000//Round

These fixation probabilities are now fairly accurate, relative to Kimura's fixation probability for N = 10000, as expected because the hitchhiking allele is neutral:

(1-Exp[-2s])/(1-Exp[-2N s])*reps/.reps->100000/.N->1000/.s->{0.002,0.005,0.01}//Round

Below, we'll plot the probability of hitchhiking among those cases where (Subscript[A, 1]) fixes (columns 1 + 3), so we need to adjust the confidence intervals (based on +/- 1.96 SE):

CIs=1.96*Sqrt[means*(1-means)/fixed]

scaledn10000=Transpose[{means,CIs}]

Graphs:

N1 = 1000;
initialp=(1/N1);

Solid curve based on full solution, dashed curve based on (((2 (e^EulerGamma)(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net])))), assuming (Subscript[ψ, net]) reasonably large.

plot1=Show[ListPlot[Table[{s,Re[Block[{$MaxExtraPrecision=1000},N[neutralnotkilledC[N1*(s+10^(-10)),N1/2000]/.p->1/N1,50]]]},{s,1/200,1/10,1/200}],Joined->True,AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[1,0,0]}],
ListPlot[Table[{s,Re[Block[{$MaxExtraPrecision=1000},N[neutralnotkilledC[N1*(s+10^(-10)),N1/1000]/.p->1/N1,50]]]},{s,1/200,1/10,1/200}],Joined->True,AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,0,1]}],
ListPlot[Table[{s,Re[Block[{$MaxExtraPrecision=1000},N[neutralnotkilledC[N1*(s+10^(-10)),N1/500]/.p->1/N1,50]]]},{s,1/200,1/10,1/200}],Joined->True,AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,1,0]}],
Plot[(((2 (E^EulerGamma)(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net]))))/.ρ->N1/2000/.(Subscript[ψ, net])->N1*s,{s,0,1/10},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[1,0,0],Dashing[0.016]}],
Plot[(((2 (E^EulerGamma)(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net]))))/.ρ->N1/1000/.(Subscript[ψ, net])->N1*s,{s,0,1/10},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,0,1],Dashing[0.016]}],
Plot[(((2 (E^EulerGamma)(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net]))))/.ρ->N1/500/.(Subscript[ψ, net])->N1*s,{s,0,1/10},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,1,0],Dashing[0.016]}]
];

plot2=ErrorListPlot[{{{0.02,scaledn1000[[1,1]]},ErrorBar[scaledn1000[[1,2]]]},{{0.05,scaledn1000[[2,1]]},ErrorBar[scaledn1000[[2,2]]]},{{0.1,scaledn1000[[3,1]]},ErrorBar[scaledn1000[[3,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[1,0,0]}];
plot2b=ErrorListPlot[{{{0.02,scaledn1000[[4,1]]},ErrorBar[scaledn1000[[4,2]]]},{{0.05,scaledn1000[[5,1]]},ErrorBar[scaledn1000[[5,2]]]},{{0.1,scaledn1000[[6,1]]},ErrorBar[scaledn1000[[6,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,0,1]}];
plot2c=ErrorListPlot[{{{0.02,scaledn1000[[7,1]]},ErrorBar[scaledn1000[[7,2]]]},{{0.05,scaledn1000[[8,1]]},ErrorBar[scaledn1000[[8,2]]]},{{0.1,scaledn1000[[9,1]]},ErrorBar[scaledn1000[[9,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,1,0]}];

Show[plot1,plot2,plot2b,plot2c]

N1 = 10000;
initialp=(1/N1);

plot1=Show[ListPlot[Table[{s,Re[Block[{$MaxExtraPrecision=1000},N[neutralnotkilledC[N1*(s+10^(-10)),N1/20000]/.p->1/N1,50]]]},{s,1/2000,1/100,1/2000}],Joined->True,AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[1,0,0]}],
ListPlot[Table[{s,Re[Block[{$MaxExtraPrecision=1000},N[neutralnotkilledC[N1*(s+10^(-10)),N1/10000]/.p->1/N1,50]]]},{s,1/2000,1/100,1/2000}],Joined->True,AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,0,1]}],
ListPlot[Table[{s,Re[Block[{$MaxExtraPrecision=1000},N[neutralnotkilledC[N1*(s+10^(-10)),N1/5000]/.p->1/N1,50]]]},{s,1/2000,1/100,1/2000}],Joined->True,AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,1,0]}],
Plot[(((2 (E^EulerGamma)(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net]))))/.ρ->N1/20000/.(Subscript[ψ, net])->N1*s,{s,0,1/100},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[1,0,0],Dashing[0.016]}],
Plot[(((2 (E^EulerGamma)(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net]))))/.ρ->N1/10000/.(Subscript[ψ, net])->N1*s,{s,0,1/100},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,0,1],Dashing[0.016]}],
Plot[(((2 (E^EulerGamma)(Subscript[ψ, net])))^(-(ρ/(Subscript[ψ, net]))))/.ρ->N1/5000/.(Subscript[ψ, net])->N1*s,{s,0,1/100},
AxesOrigin->{0,0},PlotRange->{All,{0,1}},PlotStyle->{RGBColor[0,1,0],Dashing[0.016]}]
];

plot2=ErrorListPlot[{{{0.002,scaledn10000[[1,1]]},ErrorBar[scaledn10000[[1,2]]]},{{0.005,scaledn10000[[2,1]]},ErrorBar[scaledn10000[[2,2]]]},{{0.01,scaledn10000[[3,1]]},ErrorBar[scaledn10000[[3,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[1,0,0]}];
plot2b=ErrorListPlot[{{{0.002,scaledn10000[[4,1]]},ErrorBar[scaledn10000[[4,2]]]},{{0.005,scaledn10000[[5,1]]},ErrorBar[scaledn10000[[5,2]]]},{{0.01,scaledn10000[[6,1]]},ErrorBar[scaledn10000[[6,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,0,1]}];
plot2c=ErrorListPlot[{{{0.002,scaledn10000[[7,1]]},ErrorBar[scaledn10000[[7,2]]]},{{0.005,scaledn10000[[8,1]]},ErrorBar[scaledn10000[[8,2]]]},{{0.01,scaledn10000[[9,1]]},ErrorBar[scaledn10000[[9,2]]]}},AxesOrigin->{0,0},PlotStyle->{RGBColor[0,1,0]}];


Show[plot1,plot2,plot2b,plot2c]

Supplementary material 5:  Hitchhiking in a human genome with clustering of genes

For a haploid genome of (3.*10^9/1000) kb containing (42001/2) SNPs, the average distance to the closest of the two deleterious sites (2) on either side would be exponentially distributed with mean distance in kb of:

aveDIST = (3*10^9/1000)/(42001/2)/(2)

Multiplying by Ne (10000) and using 1cM/Mb = 1/(100*1000) M/kb, we get an Ner of:

Ner=aveDIST*10000/(1000*100)

%//N

With Ne*sd = 1:

If we assumed a constant distance of Ne*r, we would get the following probabilities of hitchhiking for Ne*sa = {5, 25, 100}:

Re[N[notkilledC[Nsa-1,1,Ner]/.p->1/20000/.Nsa->{5, 25, 100}]]

If we integrated over the exponentially distributed distance of Ne*r, we would get a MUCH higher probability of hitchhiking:

NIntegrate[PDF[ExponentialDistribution[1/Ner],x]Re[notkilledC[Nsa-1,1,x]/.p->1/20000/.Nsa->5],{x,0,1000}]

NIntegrate[PDF[ExponentialDistribution[1/Ner],x]Re[notkilledC[Nsa-1,1,x]/.p->1/20000/.Nsa->25],{x,0,1000}]

NIntegrate[PDF[ExponentialDistribution[1/Ner],x]Re[notkilledC[Nsa-1,1,x]/.p->1/20000/.Nsa->100],{x,0,1000}]

If distance between selected sites were 10 times closer:

NIntegrate[PDF[ExponentialDistribution[10/Ner],x]Re[notkilledC[Nsa-1,1,x]/.p->1/20000/.Nsa->5],{x,0,1000}]

NIntegrate[PDF[ExponentialDistribution[10/Ner],x]Re[notkilledC[Nsa-1,1,x]/.p->1/20000/.Nsa->25],{x,0,1000}]

NIntegrate[PDF[ExponentialDistribution[10/Ner],x]Re[notkilledC[Nsa-1,1,x]/.p->1/20000/.Nsa->100],{x,0,1000}]

(1000 was chosen as a limit of integration, both because most of the exponential falls below this and because the hitchhiking probability is negligible above this.)

Notice that accounting for the exponential distribution of distances gives a much higher chance of hitchhiking (even if we assume that selected sites are uniformly and randomly distributed across the genome).

With Ne*sd = 10:

If we assumed a constant distance of Ne*r, we would get the following probabilities of hitchhiking for Ne*sa = {25, 100}:

aveDIST = (3*10^9/1000)/(4201/2)/(2)

Ner=aveDIST*10000/(1000*100)

%//N

Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,Ner]/.p->1/20000/.Nsa->{25,100},50]]]

Integrating over the exponetial...

For Nsa = 25:

NIntegrate[PDF[ExponentialDistribution[1/Ner],x]Re[Block[{$MaxExtraPrecision=2000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->25,200]]],{x,0,10}]

Integration over the exponential now runs into numerical problems, so the distribution was discretized into bins of width 1/100 and summed:

width=1/100;Sum[width*PDF[ExponentialDistribution[1/Ner],x]*Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->25,50]]],{x,width,10,width}]

Note that little would be added by extending the sum further (here to 11 rather than 10):

width=1/100;Sum[width*PDF[ExponentialDistribution[1/Ner],x]*Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->25,50]]],{x,10+width,11,width}]

And little is changed if we used an interval width of 1/200 instead:

width=1/200;Sum[width*PDF[ExponentialDistribution[1/Ner],x]*Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->25,50]]],{x,width,10,width}]

For Nsa = 100:

width=1/100;Sum[width*PDF[ExponentialDistribution[1/Ner],x]*Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->100,50]]],{x,width,10,width}]

Note that little would be added by extending the sum further (here to 11 rather than 10):

width=1/100;Sum[width*PDF[ExponentialDistribution[1/Ner],x]*Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->100,50]]],{x,10+width,11,width}]

And little is changed if we used an interval width of 1/200 instead:

width=1/200;Sum[width*PDF[ExponentialDistribution[1/Ner],x]*Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->100,50]]],{x,width,10,width}]

If distance between selected sites were 10 times closer...

For Nsa = 25:

width=1/100;Sum[width*PDF[ExponentialDistribution[10/Ner],x]*Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->25,50]]],{x,width,10,width}]

And little is changed if we used an interval width of 1/200 instead:

width=1/200;Sum[width*PDF[ExponentialDistribution[10/Ner],x]*Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->25,50]]],{x,width,10,width}]

For Nsa = 100:

width=1/100;Sum[width*PDF[ExponentialDistribution[10/Ner],x]*Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->100,50]]],{x,width,10,width}]

And little is changed if we used an interval width of 1/200 instead:

width=1/200;Sum[width*PDF[ExponentialDistribution[10/Ner],x]*Re[Block[{$MaxExtraPrecision=1000},N[notkilledC[Nsa-10+10^(-10),10,x]/.p->1/20000/.Nsa->100,50]]],{x,width,10,width}]

Supplementary material 6:  The case of two interfering beneficial mutations

Here we rederive our model for the case considered by Yu and Etheridge (2010), where a beneficial allele was segregating at one locus when a second beneficial allele occurred at a second case.  They were particularly interested in cases where the first beneficial allele was weaker and could be displaced by the second.

The probability of establishment for a recombinant [Yu & Etheridge case]

Using the notation of Yu & Etheridge (2010), let 11 be the double beneficial mutant haplotype that arises in a population segregating for two beneficial mutations.  Let 10 represent the haplotype bearing the first be beneficial mutations (and 01 for the second).  Following Yu & Etheridge, let's assume that the population is dominated by 10 and 01 types, with the wildtype rapidly eliminated from the population.

On pages 38 & 42 (second columns), Y&E assume that the establishment probability of the 11 recombinant is proportional to fitness(11) - mean fitness(when 11 appears).  A more accurate establishment probability can be derived along the same lines as in the Appendix (see also Sup. Mat. 1).  Here we revise our calculations for the case considered by Y&E.

Let selection favoring the first 10 mutant be s*γ, and let selection favoring the second 01 mutant be s (as in Yu & Etheridge 2010).  After the initial period where types 10 and 01 come to dominate the population and before recombination has successfully established type 11, then the frequency, p01, of 01 will rise over time according to a standard haploid model of selection (equation 1) with selective advantage (Subscript[s, net]) = s*(1-γ).  Equation (2) then gives the fixation probability of the most fit 11 double mutant, whose advantage over 10 is given by s ((Subscript[s, a])=s) and whose advantage over 01 is given by s*γ ((Subscript[s, d])=s*γ):

((2 (Subscript[s, d]) ((Subscript[s, d])+(Subscript[s, net])))/((Subscript[s, d])+p01 (Subscript[s, net])))/.(Subscript[s, net])->(Subscript[s, a])-(Subscript[s, d])/.(Subscript[s, a])->s /.(Subscript[s, d])->s*γ//Simplify

Relative to a fixation probability of 2 s, this predicts a fixation probability that depends on p01, with fixation most likely (2 s) when p01~0 because the less fit first beneficial mutation dominates and least likely (2 s γ) when p01~1 because the more fit second beneficial mutation dominates.

This equation [solid curve] can be compared to Y&E's assumed fixation probability of 2θ[t] = 2(s - s*(1-γ)*p01) [dashed curve]:

Show[
Plot[(γ/(p01+γ-p01 γ))/.γ->0.8,{p01,0,1}],
Plot[(1-(1-γ)*p01)/.γ->0.8,{p01,0,1},PlotStyle->Dashing[0.01]],
PlotRange->{{0,1},{0,1}},
AxesOrigin->{0,0}
]

For the case that Y&E consider in Figure 4, there is little difference, because γ = 0.8, meaning that the mean fitness doesn't change much as 10 and 01 haplotypes compete for fixation.  But for smaller γ, the changing mean fitness is important because it occurs over the same time scale that the 11 recombinant becomes established:

Show[
Plot[(γ/(p01+γ-p01 γ))/.γ->0.3,{p01,0,1}],
Plot[(1-(1-γ)*p01)/.γ->0.3,{p01,0,1},PlotStyle->Dashing[0.01]],
PlotRange->{{0,1},{0,1}},
AxesOrigin->{0,0}
]

Thus, the fixation probability of the 11 recombinant is lower than expected based on its current fitness, because its future selective advantage drops as the second beneficial allele (01) spreads and outcompetes the original beneficial allele (10).   Note that the discrepancy is most pronounced when the recombinant appears while the 10 and 01 haplotypes are intermediate in frequency, which makes sense because p01 changes most rapidly when near 1/2.

Solutions

To solve the integral in equation (4), we take N times the chance of a successful recombinant that appears when the frequency of 01 is p01, assuming the rest of the population is nearly all type 10:

part1=N*r/2*2*p01*(1-p01)*((2 s γ)/(p01+γ-p01 γ));

and divide this by dq/dt, which for the haploid model of selection is (assuming weak selection so that we can approximate the recursion by a differential equation):

part2=(Subscript[s, net]) * p01 (1-p01)/.(Subscript[s, net])->s*(1-γ);

Then integrating part1/part2:

Simplify[Integrate[Simplify[part1/part2],p01],{0<γ<1}]

Technically the integral should be evaluated from p01 = p0 to p01 = 1-1/N:

Simplify[((%/.p01->1-1/N)-(%/.p01->p0)),{0<γ<1}]

But given that we assume a reasonably large population that N-1 is approximately N and p01 is initially small (relative to γ), we can equally take the integral from p = 0 to 1:

Simplify[((%%/.p01->1)-(%%/.p01->0)),{0<γ<1}]

Exp[((2 N r γ Log[γ])/(((1-γ))^2))] is thus the probability that the recombinant does not form, while the probability that it does form is 1-Exp[((2 N r γ Log[γ])/(((1-γ))^2))].  This should be comparable to Yu and Etheridge's Figure 4a, but first we must scale our calculations by the establishment probability of the second beneficial mutation.

We apply the same techniques in Appendix 1 to determine the establishment probability of the second beneficial mutation, that is, a single 01 type that arises in a population of 00 and 10 individuals, where the mean fitness is changing due to the spread of the 10 type, which is at frequency p10.

This is again given by equation (2) in our paper, but now (Subscript[s, net]) describes the selective advantage of type 10 spreading in a population of 00 individuals ((Subscript[s, net])=s*γ).  When the 01 individual appears, it has selective advantage s over type 00 ((Subscript[s, a])=s) and selective advantage s - γs over type 10 ((Subscript[s, d])=s*(1-γ)):

pfix=((2 (Subscript[s, d]) ((Subscript[s, d])+(Subscript[s, net])))/((Subscript[s, d])+p10 (Subscript[s, net])))/.(Subscript[s, net])->(Subscript[s, a])-(Subscript[s, d])/.(Subscript[s, a])->s /.(Subscript[s, d])->s*(1-γ) //Simplify

We then multiply this establishment probability by the probability that a successful recombinant does form (1-Exp[((2 N r γ Log[γ])/(((1-γ))^2))] ), for parameters equivalent to Figure 4a in Yu and Etheridge.

NOTE:  The recombination rate used in Figure 4a was r = 0.00001, not r = 0.00002 as reported (Feng Yu, personal communication).  Also, while Yu and Etheridge used a Moran model, it was standardized to give the same drift and diffusion terms as the Wright-Fisher model (Feng Yu, personal communication), so their simulation and analytical results should match up to those presented here, derived based on a Wright-Fisher model (simulations that we performed confirm this).

The semi-deterministic approximation (equivalent to our equation 5) is then:

plot1=LogLinearPlot[pfix*(1-Exp[((2 N r γ Log[γ])/(((1-γ))^2))])/.s->0.02/.γ->0.8/.r->0.00001/.p10->0.2,{N,300,50000},PlotRange->{Automatic,{0,0.0225}}]

The following takes the approximation for the diffusion with killing (our equation 12) and multiplies by the establishment probability:

	pfix*(1-(e^(ω/(Subscript[ψ, d]))) (((((Subscript[ψ, d])/((Subscript[ψ, d])+(Subscript[ψ, net])))))^ω))

plot2=Show[%,LogLinearPlot[pfix*(1-(E^(ω/(Subscript[ψ, d]))) (((((Subscript[ψ, d])/((Subscript[ψ, d])+(Subscript[ψ, net])))))^ω))/.ω->((2 ρ (Subscript[ψ, d]) ((Subscript[ψ, net])+(Subscript[ψ, d])))/((Subscript[ψ, net])^2))/.(Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->γ (Subscript[ψ, a])/.ρ->N r/.(Subscript[ψ, a])->N*s/.s->0.02/.γ->0.8/.r->0.00001/.p10->0.2,{N,300,50000},PlotRange->{Automatic,{0,0.0225}},PlotStyle->Dashing[{0.01}]]
]

Here's the more accurate solution (B.6) to the diffusion with killing, assuming the 01 type is at frequency p01 and the frequency of the 10 type is nearly at frequency 1-p01 (makes little difference):

P[p01]/.Flatten[Simplify[DSolve[{((P''[p01])/2)+ψnet*(((1+(E^(-2 p01 ψnet)))/(1-(E^(-2 p01 ψnet)))))*P'[p01]-ρ*((2ψd(ψd+ψnet))/(p01 ψnet+ψd))*P[p01]==0,P'[c]==0,P[1]==1},P[p01],p01],{0<p01<1,ψnet>0,ψd>0,ρ>0}]];

notkilledC[ψnet_,ψd_,ρ_]=Simplify[%/.c->0]

Multiplying by the probability of establishment of the 01 type in a population where the 10 type is initially 2/10:

table=Table[{10.^i,(pfix/.s->2/100/.γ->8/10/.p10->2/10)*Re[Block[{$MaxExtraPrecision=1000},N[(1-notkilledC[(Subscript[ψ, net]),(Subscript[ψ, d]),ρ])/.(Subscript[ψ, net])->(Subscript[ψ, a])-(Subscript[ψ, d])/.(Subscript[ψ, d])->γ (Subscript[ψ, a])/.ρ->num r/.(Subscript[ψ, a])->num*s/.s->2/100/.γ->8/10/.r->1/100000/.p01->1/num/.num->10^i,50]]]},{i,2+1/100,5,1/4}]

The following shows that equation (12) provides an excellent approximation to the more cumbersome solution (B.6) for this case (compare dots to dashed curve):

plot3=Show[plot2,ListLogLinearPlot[table,PlotRange->All]]

The following are the simulations and analytical curves presented in Figure 4a of Yu and Etheridge (2009).  We are grateful to Feng Yu for providing us with this data.

simsYE={{500,46/100000},{1000,121/100000},{2000,277/100000},{5000,725/100000},{10000,1266/100000},{20000,1784/100000}};

analysisYE10={{200,0.00157},{500,0.00143},{1000,0.00198},{2000,0.00336},{5000,0.00709},{10000,0.0119},{20000,0.0176},{50000,0.0215}};

analysisYE16={{200,0.00079},{500,0.00122},{1000,0.00198},{2000,0.00346},{5000,0.00729},{10000,0.0121},{20000,0.0178},{50000,0.0216}};

The semi-deterministic approximation and the approximation for the diffusion both provide an excellent fit to this simulation data:

plot4=Show[plot2,ListLogLinearPlot[simsYE,PlotStyle->PointSize[0.015]],PlotRange->All]

Export["Fig.eps",%]

At low population sizes, the diffusion (dashed) fits the simulation results better than the semi-deterministic approximation (solid), as expected.

Note that the algorithm outlined in Yu and Etheridge provides a similar fit:

plot5=Show[ListLogLinearPlot[{analysisYE10,analysisYE16},Joined->True],ListLogLinearPlot[simsYE,PlotStyle->PointSize[0.015]],PlotRange->{{Log[300],Log[50000]},{0,0.0225}},AxesOrigin->{Log[300],0}]

but this algorithm is more cumbersome to evaluate and performs less well than the diffusion for small population sizes (equation (9) in Yu and Etheridge assumes that the population size is relatively large).

Simulations (r = 1/100000)

Confirming that the simulations from Yu and Etheridge correspond to the Wright-Fisher model.  [Note:  There appears to be above binomial variance (evident when rerunning cases), suggesting a problem with Mathematica's random number generator, but the trend is the same.]

simsWF={{500,0.00036},{1000,0.00106},{2000,0.00324},{5000,0.0049},{10000,0.01254}};

The semi-deterministic approximation and the approximation for the diffusion both provide an excellent fit to this simulation data:

plot4=Show[plot2,ListLogLinearPlot[{simsYE,simsWF},PlotStyle->PointSize[0.02]],PlotRange->All]

Case:  s1=0.02*0.8, s2 = 0.02, r = 1/100000, N = 500

seed=RandomInteger[10000000]

sA=2/100*8/10;
sB=2/100;
sAB=sA+sB;
r=1/100000;
NUM=500;
SeedRandom[seed];
startfreq10=2/10;startfreq01=1/NUM;

gen=0; 
startvec={1-startfreq10-startfreq01,startfreq10,startfreq01,0};
vec=startvec;
finished=False;
lost=0;
fixed01=0;
fixed11=0;
timefixed=0;
timefixed2=0;
run=1;
numruns=50000;

While[run≤numruns,

gen++;
x1:=vec[[1]];
x2:=vec[[2]];
x3:=vec[[3]];
x4:=vec[[4]];
newx1:=(x1-r*(x1*x4-x2*x3));
newx2:=(1+sA)*(x2+r*(x1*x4-x2*x3));
newx3:=(1+sB)*(x3+r*(x1*x4-x2*x3));
newx4:=(1+sAB)*(x4-r*(x1*x4-x2*x3));

total:=newx1+newx2+newx3+newx4;

newvec:={newx1,newx2,newx3,newx4}/total;

freq=Random[MultinomialDistribution[NUM,newvec]];

vec=freq/NUM;

If[vec[[3]]+vec[[4]]<10^(-16),gen=0; vec=startvec; lost=lost+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
If[vec[[3]]>1-10^(-16),timefixed=timefixed+gen;timefixed2=timefixed2+gen^2;gen=0; vec=startvec;fixed01=fixed01+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
If[vec[[4]]>1-10^(-16),timefixed=timefixed+gen;timefixed2=timefixed2+gen^2;gen=0; vec=startvec; fixed11=fixed11+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
]

{lost,fixed01,fixed11}

%/(%.{1,1,1})//N

%%/(%%.{0,1,1})//N

Mean time to fixation:

timefixed/(fixed01+fixed11)//N

SD time to fixation:

Sqrt[(timefixed2/(fixed01+fixed11)-%^2)/(fixed01+fixed11-1)]//N

SE time to fixation:

%/Sqrt[fixed01+fixed11]

Case:  s1=0.02*0.8, s2 = 0.02, r = 1/100000, N = 1000

seed=RandomInteger[10000000]

sA=2/100*8/10;
sB=2/100;
sAB=sA+sB;
r=1/100000;
NUM=1000;
SeedRandom[seed];
startfreq10=2/10;startfreq01=1/NUM;

gen=0; 
startvec={1-startfreq10-startfreq01,startfreq10,startfreq01,0};
vec=startvec;
finished=False;
lost=0;
fixed01=0;
fixed11=0;
timefixed=0;
timefixed2=0;
run=1;
numruns=50000;

While[run≤numruns,

gen++;
x1:=vec[[1]];
x2:=vec[[2]];
x3:=vec[[3]];
x4:=vec[[4]];
newx1:=(x1-r*(x1*x4-x2*x3));
newx2:=(1+sA)*(x2+r*(x1*x4-x2*x3));
newx3:=(1+sB)*(x3+r*(x1*x4-x2*x3));
newx4:=(1+sAB)*(x4-r*(x1*x4-x2*x3));

total:=newx1+newx2+newx3+newx4;

newvec:={newx1,newx2,newx3,newx4}/total;

freq=Random[MultinomialDistribution[NUM,newvec]];

vec=freq/NUM;

If[vec[[3]]+vec[[4]]<10^(-16),gen=0; vec=startvec; lost=lost+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
If[vec[[3]]>1-10^(-16),timefixed=timefixed+gen;timefixed2=timefixed2+gen^2;gen=0; vec=startvec;fixed01=fixed01+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
If[vec[[4]]>1-10^(-16),timefixed=timefixed+gen;timefixed2=timefixed2+gen^2;gen=0; vec=startvec; fixed11=fixed11+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
]

{lost,fixed01,fixed11}

%/(%.{1,1,1})//N

%%/(%%.{0,1,1})//N

Mean time to fixation:

timefixed/(fixed01+fixed11)//N

SD time to fixation:

Sqrt[(timefixed2/(fixed01+fixed11)-%^2)/(fixed01+fixed11-1)]//N

SE time to fixation:

%/Sqrt[fixed01+fixed11]

Case:  s1=0.02*0.8, s2 = 0.02, r = 1/100000, N = 2000

seed=RandomInteger[10000000]

sA=2/100*8/10;
sB=2/100;
sAB=sA+sB;
r=1/100000;
NUM=2000;
SeedRandom[seed];
startfreq10=2/10;startfreq01=1/NUM;

gen=0; 
startvec={1-startfreq10-startfreq01,startfreq10,startfreq01,0};
vec=startvec;
finished=False;
lost=0;
fixed01=0;
fixed11=0;
timefixed=0;
timefixed2=0;
run=1;
numruns=50000;

While[run≤numruns,

gen++;
x1:=vec[[1]];
x2:=vec[[2]];
x3:=vec[[3]];
x4:=vec[[4]];
newx1:=(x1-r*(x1*x4-x2*x3));
newx2:=(1+sA)*(x2+r*(x1*x4-x2*x3));
newx3:=(1+sB)*(x3+r*(x1*x4-x2*x3));
newx4:=(1+sAB)*(x4-r*(x1*x4-x2*x3));

total:=newx1+newx2+newx3+newx4;

newvec:={newx1,newx2,newx3,newx4}/total;

freq=Random[MultinomialDistribution[NUM,newvec]];

vec=freq/NUM;

If[vec[[3]]+vec[[4]]<10^(-16),gen=0; vec=startvec; lost=lost+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
If[vec[[3]]>1-10^(-16),timefixed=timefixed+gen;timefixed2=timefixed2+gen^2;gen=0; vec=startvec;fixed01=fixed01+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
If[vec[[4]]>1-10^(-16),timefixed=timefixed+gen;timefixed2=timefixed2+gen^2;gen=0; vec=startvec; fixed11=fixed11+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
]

{lost,fixed01,fixed11}

%/(%.{1,1,1})//N

%%/(%%.{0,1,1})//N

Mean time to fixation:

timefixed/(fixed01+fixed11)//N

SD time to fixation:

Sqrt[(timefixed2/(fixed01+fixed11)-%^2)/(fixed01+fixed11-1)]//N

SE time to fixation:

%/Sqrt[fixed01+fixed11]

Case:  s1=0.02*0.8, s2 = 0.02, r = 1/100000, N = 5000

seed=RandomInteger[10000000]

sA=2/100*8/10;
sB=2/100;
sAB=sA+sB;
r=1/100000;
NUM=3000;
SeedRandom[seed];
startfreq10=2/10;startfreq01=1/NUM;

gen=0; 
startvec={1-startfreq10-startfreq01,startfreq10,startfreq01,0};
vec=startvec;
finished=False;
lost=0;
fixed01=0;
fixed11=0;
timefixed=0;
timefixed2=0;
run=1;
numruns=50000;

While[run≤numruns,

gen++;
x1:=vec[[1]];
x2:=vec[[2]];
x3:=vec[[3]];
x4:=vec[[4]];
newx1:=(x1-r*(x1*x4-x2*x3));
newx2:=(1+sA)*(x2+r*(x1*x4-x2*x3));
newx3:=(1+sB)*(x3+r*(x1*x4-x2*x3));
newx4:=(1+sAB)*(x4-r*(x1*x4-x2*x3));

total:=newx1+newx2+newx3+newx4;

newvec:={newx1,newx2,newx3,newx4}/total;

freq=Random[MultinomialDistribution[NUM,newvec]];

vec=freq/NUM;

If[vec[[3]]+vec[[4]]<10^(-16),gen=0; vec=startvec; lost=lost+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
If[vec[[3]]>1-10^(-16),timefixed=timefixed+gen;timefixed2=timefixed2+gen^2;gen=0; vec=startvec;fixed01=fixed01+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
If[vec[[4]]>1-10^(-16),timefixed=timefixed+gen;timefixed2=timefixed2+gen^2;gen=0; vec=startvec; fixed11=fixed11+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
]

{lost,fixed01,fixed11}

%/(%.{1,1,1})//N

%%/(%%.{0,1,1})//N

Mean time to fixation:

timefixed/(fixed01+fixed11)//N

SD time to fixation:

Sqrt[(timefixed2/(fixed01+fixed11)-%^2)/(fixed01+fixed11-1)]//N

SE time to fixation:

%/Sqrt[fixed01+fixed11]

Case:  s1=0.02*0.8, s2 = 0.02, r = 1/100000, N = 10000

seed=RandomInteger[10000000]

sA=2/100*8/10;
sB=2/100;
sAB=sA+sB;
r=1/100000;
NUM=10000;
SeedRandom[seed];
startfreq10=2/10;startfreq01=1/NUM;

gen=0; 
startvec={1-startfreq10-startfreq01,startfreq10,startfreq01,0};
vec=startvec;
finished=False;
lost=0;
fixed01=0;
fixed11=0;
timefixed=0;
timefixed2=0;
run=1;
numruns=50000;

While[run≤numruns,

gen++;
x1:=vec[[1]];
x2:=vec[[2]];
x3:=vec[[3]];
x4:=vec[[4]];
newx1:=(x1-r*(x1*x4-x2*x3));
newx2:=(1+sA)*(x2+r*(x1*x4-x2*x3));
newx3:=(1+sB)*(x3+r*(x1*x4-x2*x3));
newx4:=(1+sAB)*(x4-r*(x1*x4-x2*x3));

total:=newx1+newx2+newx3+newx4;

newvec:={newx1,newx2,newx3,newx4}/total;

freq=Random[MultinomialDistribution[NUM,newvec]];

vec=freq/NUM;

If[vec[[3]]+vec[[4]]<10^(-16),gen=0; vec=startvec; lost=lost+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
If[vec[[3]]>1-10^(-16),timefixed=timefixed+gen;timefixed2=timefixed2+gen^2;gen=0; vec=startvec;fixed01=fixed01+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
If[vec[[4]]>1-10^(-16),timefixed=timefixed+gen;timefixed2=timefixed2+gen^2;gen=0; vec=startvec; fixed11=fixed11+1; run=run+1;If[Mod[run,500]==0,Print[run]]];
]

{lost,fixed01,fixed11}

%/(%.{1,1,1})//N

%%/(%%.{0,1,1})//N

Mean time to fixation:

timefixed/(fixed01+fixed11)//N

SD time to fixation:

Sqrt[(timefixed2/(fixed01+fixed11)-%^2)/(fixed01+fixed11-1)]//N

SE time to fixation:

%/Sqrt[fixed01+fixed11]