Supplemental files for: Computing Tree Size Under Dynamical Models of Diversification

The Interpretation of the Congruence Class Supplementary File for MacPherson and Pennell “Computing Tree Size Under Dynamical Models of Diversification” Colours and Figure Options In[]:= PlotOptions = {Frame True, FrameTicks {{True, False},{True, False}}, FrameStyle Directive[Black, 12], LabelStyle Directive[Black, 13]}; PlotTypes = {Plot, ListPlot, ListLogPlot, ListLinePlot, DiscretePlot, BarChart}; Do [Map[SetOptions[x, #] &, PlotOptions],{x, PlotTypes}]; In[]:= PCols = {PIndigo, PCyan, PTeal, PGreen, POlive, PSand, PRose, PWine, PPurple} = {RGBColor[{51, 34, 136}/255], RGBColor[{136, 204, 238}/255], RGBColor[{68, 170, 153}/255], RGBColor[{17, 119, 51}/255], RGBColor[{153, 153, 51}/255], RGBColor[{221, 204, 119}/255], RGBColor[{204, 102, 119}/255], RGBColor[{136, 34, 85}/255], RGBColor[{170, 68, 153}/255]} Out[  ] = {,,,,,,,,} In[]:= Dir =NotebookDirectory[]<>"../Figures/"; Sampling from the congruence class In[]:= parsTest = {λ05, μ03, ψ01, ρ0.1}; pTest =1; Constant In[]:= Clear [μConst] μ Const[Δμ_, pars_] :=μConst[Δμ, pars]=Block[{out}, If[Δμ0, out =Interpolation[Table[{x, μ0/. pars},{x, 0, 1, 0.01}]], out =Interpolation[Table[{x, Δμ+μ0/. pars},{x, 0, 1, 0.01}]] ]; out] In[]:= f[μs_, pars_,τ_] :=λs[τ] ((λ0-μ0-ψ0)+μs[τ]-λs[τ])+λ0*ψ0/. pars Clear [λConst] λ Const[Δμ_, pars_] :=λConst[Δμ, pars]=Block[{out}, If[Δμ0, out =Interpolation[Table[{x, λ0/. pars},{x, 0, 1, 0.01}]], out = (λs[τ2] /. NDSolve[{D[λs[τ2],τ2]f[μConst[Δμ, pars], pars, τ2], λs[0] λ0/. pars},λs[τ2],{τ2, 0, 1}]〚1〛)〚0〛 ]; out] In[]:= Clear [ψConst] ψ Const[Δμ_, pars_] :=ψConst[Δμ, pars]=Block{out}, IfΔμ  0, out =Interpolation[Table[{x, ψ0/. pars},{x, 0, 1, 0.01}]], out =InterpolationTablex, λ0ψ0/. pars λConst[Δμ, pars][x],{x, 0, 1, 0.01} ; out  In[]:= rPConst[Δμ_, pars_, t_] :=λConst[Δμ, pars][t]-μConst[Δμ, pars][t] - ψConst[Δμ, pars][t] + 1 λConst[Δμ , pars][t](D[λConst[Δμ, pars][t2], t2] /.t2t) 2 SupMatB_Congruence.nb In[]:= λ Plot =PlotEvaluate[Table[λConst[μ, parsTest][x],{μ,-3, 3, 0.5}]],{x, 0, 1}, PlotStyle colours[13], Frame True, FrameTicks {{True, False},{True, False}} (*,FrameLabel"Time (Reverse),τ","Birth Rate, λ"*) (*Export[NotebookDirectory[]<>"Figures/ModelSearch_LambdaPlot.jpeg",%];*) Out[ ] = 0.0 0.2 0.4 0.6 0.8 1.0 3 4 5 6 7 In[]:= μ Plot =PlotEvaluate[Table[μConst[μ, parsTest][x],{μ,-3, 3, 0.5}]],{x, 0, 1}, PlotStyle colours[13], Frame True, FrameTicks {{True, False},{True, False}} (*,FrameLabel"Time (Reverse),τ","Death Rate, μ"*) (*Export[ NotebookDirectory []<>" Figures /ModelSearch _ MuPlot. jpeg " , %] ; *) Out[ ] = 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 6 SupMatB_Congruence.nb 3 In[]:= ψ Plot =PlotEvaluate[Table[ψConst[μ, parsTest][x],{μ,-3, 3, 0.5}]],{x, 0, 1}, PlotStyle colours[13], Frame True, FrameTicks {{True, False},{True, False}} (*,FrameLabel"Time (Reverse),τ","Death Rate, μ"*) (*Export[NotebookDirectory[]<>"Figures/ModelSearch_PsiPlot.jpeg",%];*) Out[ ] = 0.0 0.2 0.4 0.6 0.8 1.0 0.6 0.8 1.0 1.2 1.4 1.6 1.8 In[]:= GraphicsRow [{λPlot, μPlot, ψPlot}, ImageSize Full] Out[  ] = 0.0 0.2 0.4 0.6 0.8 1.0 3 4 5 6 7 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 6 0.0 0.2 0.4 0.6 0.8 1.0 0.6 0.8 1.0 1.2 1.4 1.6 1.8 In[]:= PlotEvaluate[Table[rPConst[μ, parsTest, x],{μ,-3, 3, 0.5}]],{x, 0, 1}, PlotStyle colours[5], Frame True, FrameTicks {{True, False},{True, False}} (*,PlotRange{(λ0-μ0-ψ0)*0.75/.parsTest,(λ0-μ0-ψ0)*1.25/.parsTest}, FrameLabel"Time (Reverse),τ","Death Rate, μ"*) (* Export [ NotebookDirectory [ ] <>" Figures /ModelSearch _ PsiPlot. jpeg " , % ] ; *) Out[ ] = 0.0 0.2 0.4 0.6 0.8 1.0 0.999994 0.999996 0.999998 1.000000 1.000000 1.000000 1.000010 4 SupMatB_Congruence.nb p-trajectories (Anreoletti & Morlon 2022) First, let’s find the value of λ such that the In[]:= λ MinSol =Solve01 λ min (λ min -λ Old ) Δ t -rP +λmin -λ 0 ψ 0 λ min ,λmin Out[  ] = λmin -1-rP Δt+ (-1+rP Δt)2+4Δt(λOld +Δtλ0ψ0) 2Δt, λmin -1+rP Δt+ (-1+rP Δt)2+4Δt(λOld +Δtλ0ψ0) 2 Δ t  In[]:= Reduce[{(λmin /.λMinSol〚1〛)>0, λOld >0, Δt>0, μOld >0, λ0>0, μ0>0, ψ0>0}] Out[  ] = False In[]:= Reduce[{(λmin /.λMinSol〚2〛)>0, λOld >0, Δt>0, μOld >0, λ0>0, μ0>0, ψ0>0}] Out[ ] = rP ∈&& μOld >0&&μ0>0&&Δt>0&&λOld >0&&ψ0>0&&λ0>0 The second of these solutions is the one that is non-negative and real. In[]:= λ Min[rP_,λOld_,λ0_,ψ0_,Δt_] :=-1+rP Δt+ (-1+rP Δt)2+4Δt(λOld +Δtλ0ψ0) 2 Δ t Then we solve for the value of λi such that μ i=μi-1 In[]:= λ StarSol =SolveμOld  1 λ ( λ - λ Old ) Δt-rP +λλ 0 ψ 0 λ,λ Out[ ] = λ  -1-rP Δt-ΔtμOld + (-1+rP Δt+ΔtμOld)2+4Δt(λOld +Δtλ0ψ0) 2Δt, λ  -1+rP Δt+ΔtμOld + (-1+rP Δt+ΔtμOld)2+4Δt(λOld +Δtλ0ψ0) 2 Δ t  In[]:= Collect (-1+rP Δt+ΔtμOld) 2 +4Δt(λOld +Δtλ0ψ0) Δ t 2,{4, Δt} Out[  ] = (-1+rP Δt+ΔtμOld) 2 +4Δt(λOld +Δtλ0ψ0) Δ t 2 In[]:= Reduce[{(λ/.λStarSol〚1〛)>0, λOld >0, Δt>0, μOld >0, λ0>0, μ0>0, ψ0>0}] Out[ ] = False In[]:= Reduce[{(λ/.λStarSol〚2〛)>0, λOld >0, Δt>0, μOld >0, λ0>0, μ0>0, ψ0>0}] Out[  ] = rP ∈&& μ0>0&&μOld >0&&Δt>0&&λOld >0&&ψ0>0&&λ0>0 SupMatB_Congruence.nb 5 The second of these solutions is the one that is non-negative and real. In[]:= λ Star[rP_,λOld_,μOld_,λ0_,ψ0_,Δt_] := -1+rP Δt+ΔtμOld + (-1+rP Δt+ΔtμOld)2+4Δt(λOld +Δtλ0ψ0) 2Δ t In[]:= μ New[rP_,λNew_,λOld_,λ0_,ψ0_,Δt_] :=1 λNew λ New -λ Old Δt-rP +λNew -λ 0 ψ 0 λNew ψ New[λNew_,λ0_,ψ0_] :=λ0ψ0 λ New In[]:= Clear [Ptraj, λP, μP, ψP] Ptraj[p_,Δμ_, pars_] := Ptraj[p, Δμ, pars]=Block[{rP =λ0-μ0-ψ0/. pars, Δt=0.0001, out = {{0, λ0, μ0+Δμ,ψ0} /. pars},t,λS, λM, λi, μi, ψi}, For[t=Δt, t ≤1, t =t+Δt, λS=λStar[rP, out〚-1, 2〛, out〚-1, 3〛,λ0/. pars, ψ0/. pars, Δt]; λM=λMin[rP, out〚-1, 2〛,λ0/. pars, ψ0/. pars, Δt]; λi=out〚-1, 2〛+p(λS-out〚-1, 2〛); If[λi<λM, μi=0; ψi=ψNew[λM, λ0/. pars, ψ0/. pars];, μi=μNew[rP, λi, out〚-1, 2〛,λ0/. pars, ψ0/. pars, Δt]; ψi=ψNew[λi, λ0/. pars, ψ0/. pars]; ]; AppendTo[out, {t, λi, μi, ψi}] ]; {Interpolation[out〚;; , {1, 2}〛], Interpolation[out〚;; , {1, 3}〛], Interpolation[out〚;; , {1, 4}〛]} ] In[]:= λ P[p_,Δμ_, pars_] :=Ptraj[p, Δμ, pars]〚1〛 μ P[p_,Δμ_, pars_] :=Ptraj[p, Δμ, pars]〚2〛 ψ P[p _ , Δμ _ , pars _ ]:= Ptraj [p , Δμ , pars]〚3〛 In[]:= rPP[p_,Δμ_, pars_, t_] :=λP[p, Δμ, pars][t]-μP[p, Δμ, pars][t] - ψP[p, Δμ, pars][t] + 1 λP [ p , Δμ , pars ] [ t ] (D[λP[p, Δμ, pars][t2], t2] /.t2t) In[]:= 6 SupMatB_Congruence.nb In[]:= leg0 =LineLegendMap[Directive[Thickness[0.01],#] &, colours[5]], Table"Δμ "<> ToString[Δμ],{Δμ,-1, 1, 0.5}, Background Directive[Gray, Opacity[0.5]] Export[NotebookDirectory[] <> "Figures/BDS_Legend.jpeg", %]; Out[ ] = LineLegendcolours[Directive[Thickness[0.01], 5]], Δμ -1., Δμ -0.5, Δμ 0., Δμ 0.5, Δμ 1., Background Directive[, Opacity[0.5]]  In[]:= parsTest Out[  ] = {λ05, μ03, ψ01, ρ0.1} In[]:= λ Plot =Plot[Evaluate[Table[λP[pTest, Δμ, parsTest][t],{Δμ,-1, 1, 0.5}]], {t, 0, 1}, PlotRange {λ0*0.75 /. parsTest, λ0*1.25 /. parsTest}, PlotStyle colours[5](*,EpilogInset[leg0,Scaled[{0.8,0.8}]]*)]; Export[NotebookDirectory[] <> "Figures/BDS_LambdaPlot.jpeg", %]; μ Plot = Plot[Evaluate[Table[μP[pTest, Δμ, parsTest][t],{Δμ,-1, 1, 0.5}]],{t, 0, 1}, PlotRange {μ0*0.5 /. parsTest, μ0*1.5 /. parsTest}, PlotStyle colours[5]]; Export[NotebookDirectory[] <> "Figures/BDS_MuPlot.jpeg", %]; ψ Plot = Plot[Evaluate[Table[ψP[pTest, Δμ, parsTest][t],{Δμ,-1, 1, 0.5}]],{t, 0, 1}, PlotRange {ψ0*0.75 /. parsTest, ψ0*1.25 /. parsTest}, PlotStyle colours[5]]; Export[NotebookDirectory[] <> "Figures/BDS_PsiPlot.jpeg", %]; In[]:= GraphicsRow [{λPlot, μPlot, ψPlot}] Out[  ] = 0.0 0.2 0.4 0.6 0.8 1.0 4.0 4.5 5.0 5.5 6.0 0.0 0.2 0.4 0.6 0.8 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 0.0 0.2 0.4 0.6 0.8 1.0 0.8 0.9 1.0 1.1 1.2 Double checking the the pulled birth rates of these models are equal. SupMatB_Congruence.nb 7 In[]:= Plot[Evaluate[Table[rPP[pTest, Δμ, parsTest, t],{Δμ,-1, 1, 0.5}]],{t, 0, 1}, PlotRange {(λ0-μ0-ψ0)*0.75 /. parsTest, (λ0-μ0-ψ0)*1.5 /. parsTest}, PlotStyle colours[5]] Out[ ] = 0.0 0.2 0.4 0.6 0.8 1.0 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 Model Likelihood In[]:= xRand =Tablex, x, 1 50 , 1 - 1 50 , 1 50 ; yRand =Tablex, x, 1 50 , 1, 1 50 ; In[]:= dEdt [λ_,μ_,ψ_,τ2_] := -(λ[τ2]+μ[τ2]+ψ[τ2]) e[τ2]+λ[τ2]e[τ2]2+μ[τ2] In[]:= Clear [nsolE] nsolE[λ_,μ_,ψ_, pars_] :=nsolE[λ,μ,ψ, pars]=NDSolve[ {D[e[τ2],τ2]dEdt[λ,μ,ψ,τ2], e[0]1-ρ/. pars},{e[τ2]},{τ2, 0, 1}]〚1〛 esol [λ _ , μ _ , ψ _ , pars _ , τ _ ]:=e[τ2] /. nsolE[λ , μ , ψ , pars] /.τ2 τ In[]:= Clear [ϕArg] ϕ Arg[λ_,μ_,ψ_, pars_] :=ϕArg[λ,μ,ψ, pars]=Block[{out, tab}, tab =Table[{τ, NIntegrate[-(λ[τ2]+μ[τ2]+ψ[τ2]) + 2λ[τ2] × esol[λ,μ,ψ, pars, τ2],{τ2, 0, τ}]},{τ, 0, 1, 0.05}]; Interpolation[tab] ] In[]:= ϕ sol[λ _ ,μ _ ,ψ _ , pars _ ,τ _ ]:=Exp[ϕArg[λ,μ,ψ, pars][τ]] 8 SupMatB_Congruence.nb In[]:= GraphicsRow [ {Plot[{esol[λConst[0, parsTest],μConst[0, parsTest],ψConst[0, parsTest], {ρ0.12},τ], esol[λConst[1, parsTest],μConst[1, parsTest], ψConst[1, parsTest],{ρ0.12},τ]},{τ, 0, 1}], Plot[{ϕsol[λConst[0, parsTest],μConst[0, parsTest], ψConst[0, parsTest],{ρ0.12},τ],ϕsol[λConst[1, parsTest], μConst[1 , parsTest ] , ψConst[1 , parsTest ] , {ρ0.12} , τ]} , {τ , 0 , 1}]}] Out[ ] = 0.0 0.2 0.4 0.6 0.8 1.0 0.5 0.6 0.7 0.8 0.9 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 In[]:= LogLike[λ_,μ_,ψ_, pars_, xVec_, yVec_] :=Block[{n, out}, n=Length[xVec]-Length[yVec]+1; out =0; out += n Log[(ρ/. pars)]; out -= Log[1-esol[λ,μ,ψ, pars, 1]]; out += Log[ϕsol[λ,μ,ψ, pars, 1]]; out += Sum[Log[λ[x]]+Log[ϕsol[λ,μ,ψ, pars, x]],{x, xVec}]; out += Sum[Log[ψ[y]]-Log[ϕsol[λ,μ,ψ, pars, y]],{y, yVec}]; out] Clade Size Master Equations We can obtain this probability from the solution to the master equation: d Pn(τ) dτ-(λ(τ)+μ(τ)+ψ(τ)) Pn(τ)+λ(τ) j0 n Pj(τ)Pn-j(τ)+μδ0,n+ψδ1, n P 0(0)1ρ and P 1(0) ρ In[]:= Clear [dPdt] In[]:= dPndt [n_,λ_,μ_,ψ_,τ2_] := -(λ[τ2]+μ[τ2]+ψ[τ2]) P[n, τ2] + λ[τ2] × Sum[P[ j , τ2] × P[nj , τ2] , { j , 0 , n}]+If[n0 , μ[τ2]+ψ[τ2] , 0] SupMatB_Congruence.nb 9