Clade Size Distributions Under the Coalescent Diversification Model

Clade Size Coalescent and Yule Models Preliminaries and Notation Plotting Options Preliminaries Coalescent Model Simulations Simulating the (Kingman) coalescent up to the cut-off time. Simulating a coalescent genealogy with ni n lineages up to time T . The variable intS is dummy variable to index the stochastic replicate. In[144]:= Clear [simCoal] simCoal [nIn_, T_, intS_] :=simCoal[nIn, T, intS]=Block[{out, t, Δt, λ, l1, l2, n}, (*Initialize*) n=nIn; out =Table[{1},{i, 1, n}]; (*Print[out];*) t=0; λ=Binomial[n, 2]; Δt=RandomVariate[ExponentialDistribution[λ]]; While[t+Δt<T && n >1, t=t+Δt; (*Choose lineages to coalesce*) {l1, l2}=RandomSample[Table[i, {i, 1, n}], 2]; (*Coalescence, add joined lineages and delete parental lineages*) AppendTo[out, Join[out〚l1〛, out〚l2〛]]; out =Delete[out, {{l1},{l2}}]; (*Print[out];*) (*Update number of lineages and choose next coalescence time*) n=n-1; If[n≥2, λ=Binomial[n, 2]; Δt=RandomVariate[ExponentialDistribution[λ]]; ] ]; out ] Simulating nSim replicate coalescent simulations. In[146]:= Clear [simsC] simsC [n_, T_, nSims_] := simsC[n, T, nSims]=Table[simCoal[n, T, intS],{intS, 1, nSims}]; The size of the clades. In[364]:= simsSubCN [n_, T_, intS_] :=Map[Map[Length[#] &, #] &, {simCoal[n, T, intS]}] In[148]:= simsCN [n_, T_, nSims_] :=Map[Map[Length[#] &, #] &, simsC[n, T, nSims]] 2 CoalescentYule_11_21.nb Case 1: A randomly sampled clade In[149]:= (*Randomly sample a clade*) simCase1 [n_, T_, nSims_] :=Block{sims, samples}, sims =simsCN[n, T, nSims]; (*list of sample clade sizes*) samples =Flatten[Map[RandomSample[#, 1]&, sims]]; Tablem, Length[Select[samples, #m &]] nSims ,{m, 1, n}  In[151]:= SimPlot1 [n_, T_] :=ListPlot[simCase1[n, T, 50 000], PlotRange All, Filling Axis, PlotStyle Black, PlotMarkers {●, Medium}, FrameLabel {"Clade Size", "Prob."}] In[152]:= SimPlot1 [6, 1.0] Out[152]= ●●●●● ● 123456 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Clade Size Prob. Case 2: A randomly sampled lineage Randomly sampling a lineage is analogous to randomly sampling a clade proportional to its size. In[153]:= simCase2 [n_, T_, nSims_] :=Block  {sims, samples}, sims =simsCN[n, T, nSims]; (*list of sample clade sizes*) samples =Flatten[Map[RandomSample[# #, 1]&, sims]]; Tablem, Length[Select[samples, #m &]] nSims ,{m, 1, n}  CoalescentYule_11_21.nb 3 In[154]:= ListPlot[simCase2[6, 1, 10 000], PlotRange All, Filling Axis, FrameLabel {"Clade Size", "Prob."}] Out[154]= 123456 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Clade Size Prob. Case 3: Randomly sampling c clades (without replacement). In[155]:= simCase3 [n_, T_, c_, nSims_] :=Block  {sims, samples, i}, sims =Select[simsCN[n, T, nSims], Length[#] ≥c &]; (*Select the simulations with the minimum necessary number of clades*) samples =Map[Sort[RandomSample[#, c]] &, sims]; (*Randomly sample cclades and sort as ording does not matter*) TableM, Length[Select[samples, #M &]] nSims // N,{M, βList3[n, c]}  4 CoalescentYule_11_21.nb Case 4: Randomly sampling c lineages (without replacement). In[468]:= Clear [simSubCase4] simSubCase4 [n_, T_, c_, intS_] := simSubCase4[n, T, c, intS]=Block[{sim, samples, i, s, temp, index}, samples = {}; sim =simsSubCN[n, T, intS]〚1〛; (*Print[sim];*) index =Table[i, {i, 1, Length[sim]}]; (*Print[index];*) For[i=1, i ≤c, i++,(*Sampling of lineages*) temp =RandomSample[sim index, 1]〚1〛; (*Print[temp];*) sim〚temp〛--; (*Remove sampled lineage from clade for resampling*) AppendTo[samples, simsSubCN[n, T, intS]〚1, temp〛] ]; Sort[samples] ] In[470]:= temp4 [c_] :=Table[simSubCase4[6, 0.5, c, intS],{intS, 1, 20 000}]; test4 [c_] :=TableM, Length[Select[temp4[c],#M &]] Length [temp4[c]] // N,{M, βList4[6, c]} In[428]:= Clear [simCase4] simCase4 [n_, T_, c_, nSims_] :=Block{sim, samples, i, s, temp, index}, samples = {}; For[s=1, s ≤nSims, s++, AppendTo[samples, {}]; sim =simsCN[n, T, nSims]〚s〛; index =Table[i, {i, 1, Length[sim]}]; For[i=1, i ≤c, i++,(*Sampling of lineages*) temp =RandomSample[sim index, 1]〚1〛; sim〚temp〛--; (*Remove sampled lineage from clade for resampling*) AppendTo[samples〚-1〛, simsCN[n, T, nSims]〚s, temp〛] ]; ]; TableM, Length[Select[samples, #M &]] nSims // N,{M, βList4[n, c]}  CoalescentYule_11_21.nb 5 Clade Size Distributions Step 1 and 2 The probability of observing k ancestors at time T given that there are n samples. In[157]:= g [k_, n_, T_] :=If  k1, 1-SumExp[-Binomial[i, 2]T] (2 i -1) (-1)idF[n, i] aF[n, i],{i, 2, n}, SumExp[-Binomial[i, 2]T] (2 i -1) (-1)i-kaF[k, i -1] × dF[n, i] k! (i-k)! aF[n, i],{i, k, n}  An alternative expression for the probability above. In[158]:= g2 [k_, n_, T_] :=If  k1, 1-SumExp[-Binomial[i, 2]T]ProductBinomial[j, 2] Binomial[j, 2]-Binomial[i, 2], {j, DeleteCases[Table[x, {x, 2, n}], i]},{i, 2, n} ,1 Binomial[k, 2]Sum Binomial[i, 2]Exp[-Binomial[i, 2]T]ProductBinomial[j, 2] Binomial[j, 2]-Binomial[i, 2], {j, DeleteCases[Table[x, {x, k, n}], i]}  ,{i, k, n}   The number of ancestors as a function of time in the coalescent model 6 CoalescentYule_11_21.nb In[159]:= NumAncC = Show[Plot[{g[6, 6, t], g[5, 6, t], g[4, 6, t], g[3, 6, t], g[2, 6, t], g[1, 6, t]}, {t, 0, 2}, PlotRange All, PlotStyle MyCol2], (*Plot[Total[{g[6,6,t],g[5,6,t],g[4,6,t],g[3,6,t],g[2,6,t],g[1,6,t]}], {t,0,2},PlotStyleDirective[Black,Dashed]],*) ImageSize Medium, FrameLabel {"Time", "Prob"}] Out[159]= 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 Time Prob Step 1: Varying population size. In[]:= Let N (τ) be the population size at time τ where time is measured in units of N ref generations. The probability that n lineages have k parents in a population of changing size is given by: Three example population dynamic trajectories that have the same size at time τmax in the past. CoalescentYule_11_21.nb 7 In[160]:= Nref =5; Ndyn0 [τ_,Τmax_] :=5Exp[0(τ-Τmax)] Nref Ndyn1 [τ_,Τmax_] :=5Exp[-0.2 (τ-Τmax)] Nref Ndyn2 [τ_,Τmax_] :=5Exp[0.2 (τ-Τmax)] Nref In[164]:= Plot[{Ndyn0[τ2, 2.0], Ndyn1[τ2, 2.0], Ndyn2[τ2, 2.0]}, {τ2 , 0 , 2} , PlotStyle {{Black} , {Black , Dashed} , {Black , Dotted}}] Out[164]= 0.0 0.5 1.0 1.5 2.0 0.8 1.0 1.2 1.4 Calculating the harmonic mean size over time and representing the result as a interpolating function. In[165]:= Clear [Nh] Nh [Τmax_, Ndyn_] :=Nh[Τmax, Ndyn]=Block{out, Δt, temp}, Δt=Τmax 100. ; temp =TableΔt, 1 Ndyn[t, Τmax]Δt,{t, 0, Τmax, Δt}; temp =PrependTo[temp, {0, 0}]; Interpolation[Accumulate[temp]]  8 CoalescentYule_11_21.nb In[167]:= Clear [gTilde] gTilde [k_, n_, T_, Ndyn_, Tmax_] :=gTilde[k, n, T, Ndyn, Tmax]=Block{nRef}, Ifk1, 1SumExp[-Binomial[i, 2]Nh[Tmax, Ndyn][T]] (2 i -1) (-1)idF[n, i] aF[n, i],{i, 2, n}, Sum Exp[-Binomial[i, 2]Nh[Tmax, Ndyn][T]] (2 i -1) (-1)i-kaF[k, i -1] × dF[n, i] k! (i-k)! aF[n, i], {i, k, n}  In[]:= Clear [NumAncCDyn] In[169]:= NumAncCDyn [Ndyn_, sty_] := Show[Plot[Evaluate[Table[gTilde[k, 6, t, Ndyn, 2.0],{k, {6, 5, 4, 3, 2, 1}}]], {t, 0, 2}, PlotRange {0, 1}, PlotStyle Map[{#, sty}&, MyCol2]], ImageSize Medium , FrameLabel {"Time" , "Prob"}] In[170]:= Show [NumAncCDyn[Ndyn0, Automatic], NumAncCDyn [ Ndyn1 , Dashed] , NumAncCDyn [ Ndyn2 , Dotted]] Out[170]= 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 Time Prob Step 1: The SIS Model Calculating the dynamics of an SIRS model In[]:= Clear [pars] Pick two sets of parameters that have the same endemic outbreak size but very different transmission CoalescentYule_11_21.nb 9 In[200]:= leg2 =SwatchLegend[MyCol2, Table[m, {m, 1, 5}], Background Directive[White, Opacity[0.7]], LegendLabel Style["Mut. Count, m", Bold]]; Show PlotEvaluateTableTest1[m, 10, t] Tot[m, 10, 2.0]+Test1[10 -m, 10, t] Tot[10 -m, 10, 2.0],m, 1, 10 2, {t, 0, 2}, PlotRange All, PlotStyle MyCol2, ImageSize Medium, FrameLabel  "Time, T", "Pr(T m)",(*PlotLabel "Coalescent Model (Case 1)",*)Epilog Inset[leg2, Scaled[{0.85, 0.5}]], (*ListLinePlot[Table[{{T,0},{T,1.0}},{T,{0.5,1.0,1.5}}], PlotStyleLighter[Gray]],*)AspectRatio 0.5 Export[Dir <> "Test2. jpeg " , %] ; Out[200]= 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 2.5 Time, T Pr(T m) Mut. Count, m 1 2 3 4 5 Case 2: A randomly sampled lineage The probability of randomly sampling a lineage from a clade of size m at time T. In[202]:= Pr2[m_, n_, T_] := SumSumm n[x, n]〚m〛×g[n-kp, n, T] × Pr0[x, n, kp],{x, XTilde[n-kp, n]}, {kp, 0, n -1}  16 CoalescentYule_11_21.nb Case 2 versus Case 1 In[203]:= leg =SwatchLegend[MyCol2, Table[m, {m, 1, 6}], Background  Directive[White, Opacity[0.7]], LegendLabel Style["Clade Size", Bold]]; leg2 =LineLegend[ {Directive[Black, Thick, Opacity[0.3]], Directive[Black, Thick, Dashed]}, {"Case 1", "Case 2"}, Background Directive[White, Opacity[0.9]] (*,LegendLabelStyle["Clade Size",Bold]*)] Show [ Plot[Evaluate[Table[Pr2[m, 6, t],{m, 1, 6}]],{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#, Dashed]&, MyCol2], ImageSize Medium], Plot[Evaluate[Table[Pr1[m, 6, t],{m, 1, 6}]],{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#, Opacity[0.3]] &, MyCol2], ImageSize Medium], ListLinePlot[{{{0.5, 0},{0.5, 1}},{{1.0, 0},{1.0, 1}}}, PlotStyle Lighter[Gray]], Epilog Inset[leg2, Scaled[{0.8, 0.85}]]] Export[Dir <> "CoalCase2_A.jpeg", %]; Out[204]= Case 1 Case 2 Out[205]= 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 Case 1 Case 2 In[207]:= Case12Bar [T_, colx_, coly_] :=Show[BarChart[Table[Pr1[m, 6, T],{m, 1, 6}], ChartStyle Map[Directive[#, Opacity[0.3]] &, MyCol2], PlotRange {0, 0.5}, AspectRatio 1, ImageSize Small], Show[ListPlot[Table[{{m, Pr2[m, 6, T]}},{m, 1, 6}], PlotStyle MyCol2, PlotMarkers {●, Medium}], ListLinePlot[Table[{{m, 0},{m, Pr2[m, 6, T]}},{m, 1, 6}], PlotStyle Map[Directive[#, Dashed]&, MyCol2]]], FrameTicksStyle {{coly(*Directive[White,Opacity[0.0],12]*), False}, {colx (*Directive[Black,12]*), False}}] CoalescentYule_11_21.nb 17 In[208]:= Case12Bar [0.5, Directive[White, Opacity[0.0], 12], Directive[Black, 12]] Export[Dir <> "CoalCase2 _ B.jpeg", %]; Out[208]= ●● ●● ●● ●● ●● ●● 0.0 0.1 0.2 0.3 0.4 0.5 In[210]:= Case12Bar [1.0, Directive[Black, 12], Directive[Black, 12]] Export[Dir <> "CoalCase2_C.jpeg", %]; Out[210]= ●● ●● ●● ●● ●● ●● 123456 0.0 0.1 0.2 0.3 0.4 0.5 Validating the analytical results with the simulations In[212]:= valid2 [n_, T_] :=Show[BarChart[Table[Pr2[m, n, T],{m, 1, n}], ChartStyle Map[Directive[#, Opacity[0.75]] &, MyCol2], PlotRange {0, T}, AspectRatio 1, ImageSize Small], ListPlot[simCase2[n, T, 50 000], PlotRange All, Filling Axis, PlotStyle Black, PlotMarkers {●, Medium}]] 18 CoalescentYule_11_21.nb In[213]:= GraphicsRow [{valid2[6, 0.5], valid2[6, 1.0], valid2[6, 1.5]}] Export[Dir <> "CoalCase2 _ Valid.jpeg", %]; Out[213]= ● ●● ● ● ● 1 2 3 4 5 6 0.0 0.1 0.2 0.3 0.4 0.5 ●●●●● ● 1 2 3 4 5 6 0.0 0.2 0.4 0.6 0.8 1.0 ●●●●● ● 1 2 3 4 5 6 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Case 3: Randomly sampling c clades without replacement In[215]:= f3[x_, C_, n_, k_] := Product[dF[[x, n]〚i〛,ℬ[C, n]〚i〛],{i, 1, n}] dF[k , Length [C]] Length[C]! Product[ℬ[C , n]〚i〛! , {i , 1 , n}] In[216]:= Pr3[C_, n_, T_] := Sum[Sum[f3[x, C, n, n -kp] × g[n-kp, n, T] × Pr0[x, n, kp],{x, XTilde[n-kp, n]}], { kp , 0 , nLength [ C ] } ] Case 3 In[217]:= Plot[{Pr3[{1, 1, 2}, 6, t], Pr3[{2, 2}, 6, t], Pr3[{3, 3}, 6, t]}, { t , 0 , 2 } , PlotRange All , PlotStyle  CompCol ] Out[217]= 0.0 0.5 1.0 1.5 2.0 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 CoalescentYule_11_21.nb 19 Case 3 vs. Case 1 In[]:= ? LineLegend Out[ ] = Symbol LineLegend [{col1,…},{lbl1,…}] generates a legend that associates color col iwith label lbl i. LineLegend [{col1,…}, Automatic ]generates a legend with placeholder labels for the colors col i. LineLegend [{lbl1,…}] represents a legend with inherited colors within visualization functions. In[218]:= legC =LineLegendFlatten[ Map[{Directive[#, Dashed], Directive[#, Opacity[0.5], Thick]} &, CompCol]], "Pr((2, 2)6, T)", "Pr (2 6, T)2", "Pr((3, 3)6, T)", "Pr (3 6, T)2", "Pr((2, 2, 2)6, T)", "Pr (2 6, T)3", "Pr((2, 3)6, T)", "Pr (2 6, T)Pr (3 6, T)", Background Directive[{White, Opacity[0.95]}], ImageSize {110, Automatic} Export [ Dir <> "CoalCase13 _ Leg . jpeg " , % ] ; Out[218]= Pr((2, 2)6, T) Pr (2 6, T)2 Pr((3, 3)6, T) Pr (3 6, T)2 Pr((2, 2, 2)6, T) Pr (2 6, T)3 Pr((2, 3)6, T) Pr (2 6, T)Pr (3 6, T) 20 CoalescentYule_11_21.nb In[220]:= Case13 =Show Plot[{Pr3[{2, 2}, 6, t], Pr3[{3, 3}, 6, t], Pr3[{2, 2, 2}, 6, t], Pr3[{2, 3}, 6, t]}, {t, 0, 2}, PlotRange All, PlotStyle Map[Directive[{#, Dashed}] &, CompCol]], PlotPr1[2, 6, t]2, Pr1[3, 6, t]2, Pr1[2, 6, t]3, Pr1[2, 6, t] × Pr1[3, 6, t],{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[{#, Opacity[0.5]}] &, CompCol] (*,EpilogInset[legC,Scaled[{0.8,0.65}]]*) Export[Dir <> "CoalCase13.jpeg", %]; Out[220]= 0.0 0.5 1.0 1.5 2.0 0.00 0.02 0.04 0.06 0.08 0.10 CoalescentYule_11_21.nb 21 Validation with simulations In[222]:= temp =simCase3[6, 0.5, 2, 50 000]〚;; , 2〛; Valid3a =Show[BarChart[Map[Pr3[#, 6, 0.5]&, βList3[6, 2]], ChartStyle Directive[CompCol〚1〛, Opacity[0.5]]], ListPlot[Table[{j, temp〚j〛},{j, 1, Length[temp]}], PlotStyle CompCol〚1〛, PlotMarkers {●, Medium}], ListLinePlot[Table[{{j, 0},{j, temp〚j〛}},{j, 1, Length[temp]}], PlotStyle Directive[CompCol〚1〛, Dashed], PlotRange All], FrameTicks {{True, False},{tikList[βList3[6, 2]], False}}] Export[Dir <> "CoalCase3_Valida.jpeg", %]; Out[223]= ● ● ● ● ● ● ● ● ● (1, 1) (1, 2) (1, 3) (1, 4) (1, 5) (2, 2) (2, 3) (2, 4) (3, 3) 0.00 0.05 0.10 0.15 22 CoalescentYule_11_21.nb In[225]:= temp =simCase3[6, 0.5, 3, 50 000]〚;; , 2〛; Valid3b =Show[BarChart[Map[Pr3[#, 6, 0.5]&, βList3[6, 3]], ChartStyle Directive[CompCol〚2〛, Opacity[0.5]]], ListPlot[Table[{j, temp〚j〛},{j, 1, Length[temp]}], PlotStyle CompCol〚2〛, PlotMarkers {●, Medium}], ListLinePlot[Table[{{j, 0},{j, temp〚j〛}},{j, 1, Length[temp]}], PlotStyle Directive[CompCol〚2〛, Dashed], PlotRange All], FrameTicks {{True, False},{tikList[βList3[6, 3]], False}}] Export[Dir <> "CoalCase3_Validb.jpeg", %]; Out[226]= ● ●● ● ● ● ● (1, 1, 1) (1, 1, 2) (1, 1, 3) (1, 1, 4) (1, 2, 2) (1, 2, 3) (2, 2, 2) 0.00 0.05 0.10 0.15 0.20 0.25 CoalescentYule_11_21.nb 23 In[228]:= temp =simCase3[6, 0.5, 4, 50 000]〚;; , 2〛; Valid3c =Show[BarChart[Map[Pr3[#, 6, 0.5]&, βList3[6, 4]], ChartStyle Directive[CompCol〚3〛, Opacity[0.5]]], ListPlot[Table[{j, temp〚j〛},{j, 1, Length[temp]}], PlotStyle CompCol〚3〛, PlotMarkers {●, Medium}], ListLinePlot[Table[{{j, 0},{j, temp〚j〛}},{j, 1, Length[temp]}], PlotStyle Directive[CompCol〚3〛, Dashed], PlotRange All], FrameTicks {{True, False},{tikList[βList3[6, 4]], False}}] Export[Dir <> "CoalCase3_Validc.jpeg", %]; Out[229]= ● ● ● ● (1, 1, 1, 1) (1, 1, 1, 2) (1, 1, 1, 3) (1, 1, 2, 2) 0.00 0.02 0.04 0.06 0.08 0.10 Case 4: Randomly sampling c lineages without replacement In[]:= (*f4[x_,C_,n_]:=Block[{iList,c}, c=Length[C]; iList=Table[Length[Select[C,#C〚i〛&]],{i,1,c}]; Product[Binomial[[x,n]〚i〛,ℬ[C,n]〚i〛],{i,1,n}] Product[Binomial[C〚j〛,iList〚j〛],{j,1,c}]/Binomial[n,c] ]*) In[231]:= f4[x_, C_, n_] :=Block{iList, c}, c=Length[C]; Product[dF[i[x, n]〚i〛,ℬ[C, n]〚i〛],{i, 1, n}] Product[ℬ[C, n]〚i〛!,{i, 1, n}] Length[C]! dF[n, Length[C]]  In[234]:= Pr4[C_, n_, T_] :=Sum[ Sum[f4[x, C, n] × g[n-kp, n, T] × Pr0[x, n, kp],{x, XTilde[n-kp, n]}],{kp, 0, n}] 24 CoalescentYule_11_21.nb In[235]:= Show [ Plot[{Pr3[{1, 1, 2}, 6, t], Pr3[{2, 2}, 6, t], Pr3[{3, 3}, 6, t]},{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#, Dashed]&, CompCol]], Plot[{Pr4[{1, 1, 2}, 6, t], Pr4[{2, 2}, 6, t], Pr4[{3, 3}, 6, t]},{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#, Opacity[0.3], Thick]&, CompCol]]] Export[Dir <> "CoalCase34.jpeg", %]; Out[235]= 0.0 0.5 1.0 1.5 2.0 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 In[237]:= legD =LineLegendFlatten[Map[ {Directive[#, Dashed], Directive[#, Opacity[0.5], Thick]} &, CompCol〚1 ;; 3〛]], "Pr((1, 1, 2)6, T)", "Pr3((1, 1, 2)6, T)", "Pr((2, 2)6, T)", "Pr2((2, 2)6, T)", "Pr((3, 3)6, T)", "Pr2((3, 3)6, T)", Background Directive[{White, Opacity[0.95]}], ImageSize {110, Automatic} Export[Dir <> "CoalCase34 _ Leg . jpeg " , %] ; Out[237]= Pr((1, 1, 2)6, T) Pr3((1, 1, 2)6, T) Pr((2, 2)6, T) Pr2((2, 2)6, T) Pr((3, 3)6, T) Pr2((3, 3)6, T) In[239]:= CList = {{1, 1, 2},{2, 2},{3, 3}}; CoalescentYule_11_21.nb 25 In[301]:= leg3 =ShowPlotiDynτ, tMax iEqu /. pars1/. pars1, τ, 0, tMax iEqu /. pars1, PlotRange All, PlotStyle Black, PlotiDynτ,tMax iEqu /. pars2/. pars2, τ, 0, tMax iEqu /. pars2, PlotRange All, PlotStyle Directive[Black, Dashed], Background Directive[White, Opacity[0.8]], Axes False, FrameTicks {{{0.5, 1.0}, False},{{0, 1.0, 2.0}, False}}, FrameTicksStyle Directive[Black, 7], FrameLabel {Style["(Coalescent)Time", 8], Style["(Scaled)Prevalence", 8]}, ImageSize {150, 100} Out[301]= 0 1. 2. 0.5 1. (Coalescent)Time (Scaled)Prevalence Plot of the ratio of harmonic mean prevalence In[332]:= ListLinePlot   Tableτ2, 1 NIntegrate1 iDynτ,tMax iEqu /.pars1/.pars1 ,{τ,0,τ2} 1 NIntegrate1 iDynτ,tMax iEqu /.pars2/.pars2 ,{τ,0,τ2} ,τ2, 0.1, tMax iEqu /. pars1, 0.1, FrameLabel {"Time", "Ratio of Harmonic Mean Prevelances"}  Out[332]= 0.0 0.5 1.0 1.5 2.0 1.05 1.10 1.15 1.20 1.25 1.30 Time Ratio of Harmonic Mean Prevelances 32 CoalescentYule_11_21.nb In[302]:= Show  (*Large Outbreak*) ListLinePlotTableTable{t, Pr4SIR[m, 6, t, pars1]}, t, 0, tMax iEqu /. pars1, tMax iEqu 50 /. pars1,{m, 1, 6}, PlotRange All, PlotStyle Map[{#, Automatic}&, MyCol2], ImageSize Medium, (*Mild Outbreak*) ListLinePlotTableTable{t, Pr4SIR[m, 6, t, pars2]}, t, 0, tMax iEqu /. pars2, tMax iEqu 50 /. pars2,{m, 1, 6}, PlotRange All, PlotStyle Map[{#, Dashed}&, MyCol2], ImageSize Medium, (*Time points*) ListLinePlot[ {{{0.5, 0},{0.5, 1.0}},{{1.25, 0},{1.25, 1.0}}}, PlotStyle Gray], Epilog Inset[leg3, Scaled[{0.38, 0.75}]]  Export[Dir <> "CoalCase5_SIR.jpeg", %]; Out[302]= 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 0 1. 2. 0.5 1. (Coalescent)Time (Scaled)Prevalence CoalescentYule_11_21.nb 33 In[304]:= Case5SIRBar [T_, parsA_, parsB_, colx_, coly_] := Show[BarChart[Table[Pr4SIR[m, 6, T, parsA],{m, 1, 6}], ChartStyle Map[Directive[#, Opacity[0.3]] &, MyCol2], PlotRange {0, 1.0}, AspectRatio 1, ImageSize Small], Show[ListPlot[Table[{{m, Pr4SIR[m, 6, T, parsB]}},{m, 1, 6}], PlotStyle MyCol2, PlotMarkers {●, Medium}], ListLinePlot[Table[{{m, 0},{m, Pr4SIR[m, 6, T, parsB]}},{m, 1, 6}], PlotStyle Map[Directive[#, Dashed]&, MyCol2], PlotRange All]], FrameTicksStyle {{coly(*Directive[White,Opacity[0.0],12]*), False}, {colx (*Directive[Black,12]*), False}}] In[305]:= Case5SIRBar [0.5, pars1, pars2, Directive[White, Opacity[0.0], 12], Directive[Black, 12]] Export [ Dir <> "CoalCase5 _ SIR _ B. jpeg " , % ] ; Out[305]= ●● ●● ●● ●● ●● ●● 0.0 0.2 0.4 0.6 0.8 1.0 In[307]:= Case5SIRBar [1.25, pars1, pars2, Directive[Black, 12], Directive[Black, 12]] Export[Dir <> "CoalCase5 _ SIR _ C. jpeg ", %]; Out[307]= ●● ●● ●● ●● ●● ●● 123456 0.0 0.2 0.4 0.6 0.8 1.0 Yule 34 CoalescentYule_11_21.nb Deriving The Number of Ancestors at Time T Prob. of k ancestors at time T Deriving Master Equations for the number of ancestors In[]:= F[k_,τ_, n_] := If[kn , -λk P[k , τ] , If[k1 , λ(k+1)P[k+1 , τ] , -λk P[k , τ]+λ(k+1)P[k+1 , τ]]] In[]:= equs [n_] :=Join[Table[D[P[k, τ],τ]F[k, τ, n],{k, 1, n}], Table[P[k, 0]If[kn, 1, 0],{k, 1, n}]] In[]:= vars [n_] :=Table[P[k, τ],{k, 1, n}] Semi-analytical solution given n . In[]:= sol [n_] :=DSolve[equs[n], vars[n],τ] In[]:= sol [2] Out[ ] = P[1, τ] -2λ τ -1+2λ τ, P[2, τ] -2λ τ In[]:= sol [3] Out[  ] =   P[1, τ] -3λ τ  -1+λ τ  2  2+λ τ  , P[2, τ]3-3λ τ  -1+λ τ  , P[3, τ] -3λ τ   In[]:= sol [4] Out[ ] = P[1, τ] -4λ τ -1+λ τ 3 3+λ τ, P[2, τ]6-4λ τ  -1+λ τ  2, P[3, τ]4-4λ τ  -1+λ τ  , P[4, τ] -4λ τ   In[]:= sol [6] Out[  ] = P[1, τ] -6λ τ -1+λ τ 5 5+λ τ, P[2, τ]15 -6λ τ -1+λ τ4, P[3, τ]20 -6λ τ -1+λ τ3, P[4, τ]15 -6λ τ  -1+λ τ  2, P[5, τ]6-6λ τ  -1+λ τ  , P[6, τ] -6λ τ   Numerical Solution CoalescentYule_11_21.nb 35 In[]:= tempN =NDSolve[equs[6] /.λ1.0, vars[6],{τ, 0, 1.0}]〚1〛; NumPlot =Plot[{P[6, τ] /. tempN, P[5, τ] /. tempN, P[4, τ] /. tempN, P[3, τ] /. tempN, P[2, τ] /. tempN, P[1, τ] /. tempN, P[6, τ]+P[5, τ]+P[4, τ]+P[3, τ]+P[2, τ]+P[1, τ] /. tempN}, {τ, 0, 1}, PlotStyle Join[Map[Directive[#, Thickness[0.003]] &, MyCol], {Directive[Black, Dashed]}]] Out[  ] = 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Analytical (Sequential) Solution In[]:= DSolve [{D[P[n, τ],τ]-λn P[n, τ], P[n, 0]1},{P[n, τ]},τ] Out[ ] =   P[n, τ] -nλ τ   In[]:= DSolve D[P[n-1, τ],τ]-λ(n-1)P[n-1, τ]+λn-nλ τ, P[n-1, 0]0, {P[n-1, τ]},τ  // Simplify Out[  ] =   P[-1+n, τ] -nλ τ  -1+λ τ  n   In[]:= DSolve D[P[n-2, τ],τ]-λ(n-2)P[n-2, τ]+λ(n-1)-nλ τ -1+λ τn, P[n-2, 0]0  ,{P[n-2, τ]},τ  // Simplify Out[  ] = P[-2+n, τ]1 2 -nλ τ -1+λ τ2(-1+n)n General solution for k > n In[]:= PrTest[n_, kp_,λ_,τ_] := 1 kp !dF[n, kp]Exp[-nλ τ] (Exp[λ τ]-1)kp In[]:= PrTest[n, n -2, λ,τ]//Simplify Out[ ] = -nλ τ -1+λ τ- 2 +nn! 2 ( -2+n ) ! 36 CoalescentYule_11_21.nb In[]:= DSolve D[P[1, τ],τ] λ 2-nλ τ -1+λ τ- 2 +nn! 2(-2+n)! /.n6, P[1, 0]0, {P[1, τ]},τ  // Simplify Out[  ] =   P[1, τ] -6λ τ  -1+λ τ  5  5+λ τ    General solution for k 1 In[]:= PrTest2[n_, kp_,λ_,τ_] :=If  kp n-1, -nλ τ -1+λ τ(n-1)(n-1)+λ τ,1 kp !dF[n, kp]Exp[-nλ τ] (Exp[λ τ]-1)kp Numerical check In[]:= Show [NumPlot, Plot[{PrTest2[6, 0, 1.0, τ], PrTest2[6, 1, 1.0, τ], PrTest2[6, 2, 1.0, τ], PrTest2[6, 3, 1.0, τ], PrTest2[6, 4, 1.0, τ], PrTest2[6, 5, 1.0, τ]},{τ, 0, 1.0}, PlotStyle Map[Directive[#, Dashed, Thick]&, MyCol], PlotRange All]] Out[  ] = 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Simulation Test Simulating the Yule (coalescent) up to cut-off time T . Case 1: A randomly sampled clade Case 2: A randomly sampled lineage Case 3: Randomly sampling c clades (without replacement). Case 4: Randomly sampling c lineages (without replacement). CoalescentYule_11_21.nb 37 Clade Size Distributions Step 1 and 2: ◼ P r(mλ,T): The probability of a clade of size m descending from an ancestor at time T . In[]:= PrY[m_,λ_, T_] :=Exp[-λT] (1-Exp[-λT])m-1 In[]:= PlotEvaluate[Table[PrY[m, 1.0, t],{m, 1, 6}]],{t, 0, 2.0}, PlotRange All, FrameLabel  "Time", "Clade Size (k1)", PlotStyle MyCol2 Out[ ] = 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 1.0 Time Clade Size (k1) Pr (m n,T) x Pr (m x,n,T) k1 n P r (x k,n,T)P r (k n,T) ◼Pr(k n,T): The probability that the Yule process with n lineages at the present day has k ancestors at time T before the present day is: In[]:= PrYk[k_, n_,λ_, T_] :=If  k1, -nλT-1+λT(n-1)(n-1)+λT, 1 ( n-k ) !dF[n, (n-k)] Exp[-nλT] (Exp[λT]-1)(n-k) ◼Pr(x k,n,T): Next we wish to calculate the probability of observing a given partition x of the n taxa at time T given that there are k ancestors at that time. In[]:= Clear [PrY0] In[]:= PrY0[x_, k_, n_] := k! Binomial [ n-1 , k-1 ] Product [  [ x , n ] 〚 i 〛 ! , { i , 1 , Length [  [ x , n ] ] } ] Case 1: A randomly sampled clade 38 CoalescentYule_11_21.nb Case 2: A randomly sampled lineage Case 3: Randomly sampling multiple clades (without replacement) Case 4: Randomly sampling multiple lineages (without replacement) Coalescent Yule Comparison Time to the most recent common ancestor Comparing the times to common ancestry in the Yule versus Coalescent Models. In[]:= Plot{g[1, 6, t], PrYk[1, 6, 1.0, t]},{t, 0, 2.0}, PlotStyle {Gray, Black}, FrameLabel  "Time, T", "Pr(TMRCA ≤T)", Epilog InsetLineLegend{Gray, Black},  "Coalescent", "Yule (λ1.0)"   , Scaled[{0.25, 0.75}]   Out[  ] = 0.0 0.5 1.0 1.5 2.0 0.0 0.2 0.4 0.6 0.8 Time, T Pr(TMRCA ≤T) Coalescent Yule (λ  1.0) Comparison of Coalescent and Yule Clade-Size Distributions Density-Dependence In[]:= λ Comp =1.0; CoalescentYule_11_21.nb 39 In[]:= GraphicsRow [{Show[ Plot[{Pr3[{1, 1, 2}, 6, t], Pr3[{2, 2}, 6, t], Pr3[{3, 3}, 6, t]},{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#, Dashed]&, CompCol]], Plot[{Pr4[{1, 1, 2}, 6, t], Pr4[{2, 2}, 6, t], Pr4[{3, 3}, 6, t]}, {t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#, Opacity[0.3], Thick]&, CompCol]]], Show[ Plot[{PrY3[{1, 1, 2},6,λComp, t], PrY3[{2, 2},6,λComp, t], PrY3[{3, 3},6,λComp, t]},{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#, Dashed]&, CompCol]], Plot[{PrY4[{1, 1, 2},6,λComp, t], PrY4[{2, 2},6,λComp, t], PrY4[{3, 3},6,λComp, t]},{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#, Opacity[0.3], Thick]&, CompCol]]]}] Out[  ] = 0.0 0.5 1.0 1.5 2.0 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.0 0.5 1.0 1.5 2.0 0.00 0.05 0.10 0.15 0.20 0.25 0.30 40 CoalescentYule_11_21.nb In[]:= Show [ (*Coalescent Clade Sampling*) Plot[{Pr3[{1, 1, 2}, 6, t], Pr3[{2, 2}, 6, t], Pr3[{3, 3}, 6, t]},{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#, Opacity[0.3]] &, CompCol]], (*Coalescent Lineage Sampling*) Plot[{Pr4[{1, 1, 2}, 6, t], Pr4[{2, 2}, 6, t], Pr4[{3, 3}, 6, t]},{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#, Dashed, Thick]&, CompCol]], (*Yule Lineage Sampling*) Plot[{PrY4[{1, 1, 2},6,λComp, t], PrY4[{2, 2},6,λComp, t], PrY4[{3, 3},6,λComp, t]},{t, 0, 2}, PlotRange All, PlotStyle Map[Directive[#(* , Opacity[0.3]*) , Thick , Dotted]& , CompCol]]] Out[ ] = 0.0 0.5 1.0 1.5 2.0 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Simulation Test In[]:= simYB [6, 1.0, 1] Out[ ] = {{0, {1, 1, 1, 1, 1, 1}},{0.017626, {2, 1, 1, 1, 1}},{0.142699, {3, 1, 1, 1}}, {0.275207, {3, 2, 1}},{0.407731, {3, 3}},{0.45118, {6}},{3., {7}}} Time T belongs to what interval of the piecewise constant function? In[]:= part[n_,λ_, T_, intS_] := Floor[Interpolation[Transpose[{simYB[n, λ, intS]〚;; , 1〛, Table[i, {i, 0, n}]}], InterpolationOrder 0 ] [ T ] ] In[]:= Clear [randClade] randClade[n_,λ_, T_, intS_] :=randClade[n, λ, T, intS] = RandomSample[simYB[n, λ, intS]〚part[n, λ, T, intS], 2〛, 1]〚1〛 CoalescentYule_11_21.nb 41