ODESystem = {
x'[t]==(Subscript[Λ, H])*(1-x[t])+(Subscript[β, AH])*y[t]*(1-x[t])+(Subscript[β, HH])*x[t]*(1-x[t])-(Subscript[μ, H])*x[t],
y'[t]==(Subscript[Λ, A])*(1-y[t])+(Subscript[β, HA])*x[t]*(1-y[t])+(Subscript[β, AA])*y[t]*(1-y[t])-(Subscript[μ, A])*y[t]
}

soln1=NDSolve[{ODESystem/.{(Subscript[β, HA])->0.1,(Subscript[β, AH])->0.1,(Subscript[β, HH])->0.1, (Subscript[β, AA])->0.1,(Subscript[μ, A])->0.1,(Subscript[μ, H])->0.1, (Subscript[Λ, A])->0.1,(Subscript[Λ, H])->0.1} , x[0]==0, y[0]==0}, {x, y}, {t,0,200}];

(*Figure 1A*)
Plot[{Evaluate[x[t]/.soln1[[1,1]]], Evaluate[y[t]/.soln1[[1,2]]]}, {t,0,200}, PlotRange->{{0,200},{0.0,0.72}}, PlotStyle->{Dashing[0.01], Dashing[0.03]},PlotLegends->None,PlotTheme->{"Detailed", "ThickLines"},PlotRange->All, ImageSize->1024, FrameLabel->{"Time", "Fraction of population"}, LabelStyle->Directive[Bold, Large]]

soln2=NDSolve[{ODESystem/.{(Subscript[β, HA])->0.001,(Subscript[β, AH])->0.1,(Subscript[β, HH])->0.1, (Subscript[β, AA])->0.1,(Subscript[μ, A])->0.1,(Subscript[μ, H])->0.1, (Subscript[Λ, A])->0.1,(Subscript[Λ, H])->0.1} , x[0]==0, y[0]==0}, {x, y}, {t,0,200}];

(*Figure 1B*)
Plot[{Evaluate[x[t]/.soln2[[1,1]]], Evaluate[y[t]/.soln2[[1,2]]]}, {t,0,200}, PlotRange->{{0,200},{0.0,0.72}}, PlotStyle->{Dashing[0.01], Dashing[0.03]}, PlotLegends->None,PlotTheme->{"Detailed", "ThickLines"},PlotRange->All, ImageSize->1024, FrameLabel->{"Time", "Fraction of population"}, LabelStyle->Directive[Bold, Large]]

(*Sensitivity Analysis, Figures 2 & 4 are created in R*)

ODESystem2 = {
lh*(1-x[t])+bah*y[t]*(1-x[t])+bhh*x[t]*(1-x[t])-mh*x[t],
la*(1-y[t])+bha*x[t]*(1-y[t])+baa*y[t]*(1-y[t])-ma*y[t]
}

equilibrium [bha1_, bah1_, bhh1_, baa1_, lh1_, la1_, mh1_, ma1_]:= NSolve[{ODESystem2[[1]]==0, ODESystem2[[2]]==0&& x[t]≥ 0&&y[t]≥0}/.{bha->bha1,bah->bah1,bhh->bhh1, baa->baa1,ma->ma1,mh-> mh1,la->la1,lh->lh1}, {x[t], y[t]}, Reals][[1,1,2]];

impact[bha1_, bah1_, bhh1_, baa1_, lh1_, la1_, mh1_, ma1_]:=1-(NSolve[{ODESystem2[[1]]==0, ODESystem2[[2]]==0&& x[t]≥ 0&&y[t]≥0}/.{bha->bha1,bah->bah1,bhh->bhh1, baa->baa1,ma->ma1,mh-> mh1,la->0.0,lh->lh1}, {x[t], y[t]}, Reals][[1,1,2]]/NSolve[{ODESystem2[[1]]==0, ODESystem2[[2]]==0&& x[t]≥ 0&&y[t]≥0}/.{bha->bha1,bah->bah1,bhh->bhh1, baa->baa1,ma->ma1,mh-> mh1,la->la1,lh->lh1}, {x[t], y[t]}, Reals][[1,1,2]])

parameters=Rest[Import[NotebookDirectory[]<>"paras_sens_analysis_manuscript.csv"]];

outcomeequilibrium = Table[equilibrium[parameters[[x,2]],parameters[[x,3]],parameters[[x,4]],parameters[[x,5]],parameters[[x,6]],parameters[[x,7]],parameters[[x,8]],parameters[[x,9]]], {x,1,Length[parameters]}];

Export[NotebookDirectory[]<>"outcomeequilibrium.csv", outcomeequilibrium];

outcomeimpact = Table[impact[parameters[[x,2]],parameters[[x,3]],parameters[[x,4]],parameters[[x,5]],parameters[[x,6]],parameters[[x,7]],parameters[[x,8]],parameters[[x,9]]], {x,1,Length[parameters]}];

Export[NotebookDirectory[]<>"outcomeimpact.csv", outcomeimpact];

(*Figure 3A*)
Plot3D[equilibrium[0.1,par1,0.1,0.1,0.1,par2,0.1,0.1],{par1,0,1},{par2,0,1}, AxesLabel->{"(Subscript[β, AH])","(Subscript[Λ, A])", "(Subscript[R, H]^( * ))"}, ColorFunction->Function[{x,y,z},Hue[.65(1-Rescale[z, {0.6180339887498953`, 0.9},{0,1}])]], ColorFunctionScaling->False,PlotRange->{{0,1},{0,1},{0.6,0.9}},
ImageSize->Full,AxesStyle->Directive[Bold, FontFamily->"Arial", FontSize->32], LabelStyle->Directive[Bold, FontFamily->"Arial"] ]
Plot3D[NSolve[{ODESystem2[[1]]==0, ODESystem2[[2]]==0&& x[t]≥ 0&&y[t]≥0}/.{bha->0.1,bah->par1,bhh->0.1, baa->0.1,ma->1.0,mh-> 0.1,la->par2,lh->0.1}, {x[t], y[t]}, Reals][[1,1,2]],{par1,0,1},{par2,0,1}, AxesLabel->{"(Subscript[β, AH])","(Subscript[Λ, A])", "(Subscript[R, H]^( * ))"}, ColorFunction->Function[{x,y,z},Hue[.65(1-Rescale[z, {0.6180339887498953`, 0.9},{0,1}])]], ColorFunctionScaling->False,PlotRange->{{0,1},{0,1},{0.6,0.9}},
ImageSize->Full,AxesStyle->Directive[Bold, FontFamily->"Arial", FontSize->32], LabelStyle->Directive[Bold, FontFamily->"Arial"] ]

(*Figure 3B*)
Plot3D[NSolve[{ODESystem2[[1]]==0, ODESystem2[[2]]==0&& x[t]≥ 0&&y[t]≥0}/.{bha->0.001,bah->par1,bhh->0.1, baa->0.1,ma->1.0,mh-> 0.1,la->par2,lh->0.1}, {x[t], y[t]}, Reals][[1,1,2]],{par1,0,1},{par2,0,1}, AxesLabel->{"(Subscript[β, AH])","(Subscript[Λ, A])", "(Subscript[R, H]^( * ))"}, ColorFunction->Function[{x,y,z},Hue[.65(1-Rescale[z, {0.6180339887498953`, 0.9},{0,1}])]], ColorFunctionScaling->False,PlotRange->{{0,1},{0,1},{0.6,0.9}},ImageSize->Full,AxesStyle->Directive[Bold, FontFamily->"Arial", FontSize->32], LabelStyle->Directive[Bold, FontFamily->"Arial"] ]

(*Legend for Figure 3*)
BarLegend[{Hue[.65(1-Rescale[#, {0.6, 0.9},{0,1}])]&,{0.6,0.9}},LegendLabel->"(Subscript[R, H]^( * ))", ColorFunctionScaling->False, LabelStyle->Directive[Bold, FontFamily->"Arial", FontSize->32], LegendMarkerSize->400]

(*Figure 5A*)
Plot3D[impact[0.1,par1,0.1,0.1,0.1,par2,0.1,0.1], {par1,0,1}, {par2,0,1},AxesLabel->{"(Subscript[β, AH])","(Subscript[Λ, A])", "ω"},ColorFunction->Function[{x,y,z},Hue[.65(1-z/0.23)]] ,ColorFunctionScaling->False,PlotRange->{{0,1.0},{0,1.0},{0,0.3}}, ImageSize->Full,AxesStyle->Directive[Bold, FontFamily->"Arial", FontSize->32], LabelStyle->Directive[Bold, FontFamily->"Arial"]]

(*Figure 5B *)
Plot3D[impact[0.001,par1,0.1,0.1,0.1,par2,0.1,0.1], {par1,0,1}, {par2,0,1},AxesLabel->{"(Subscript[β, AH])","(Subscript[Λ, A])", "ω"},ColorFunction->Function[{x,y,z},Hue[.65(1-z/0.23)]],ColorFunctionScaling->False, PlotRange->{{0,1.0},{0,1.0},{0,0.3}}, ImageSize->Full,AxesStyle->Directive[Bold, FontFamily->"Arial", FontSize->32], LabelStyle->Directive[Bold, FontFamily->"Arial"]]

(*Legend for Figure 5*)
BarLegend[{Hue[0.65*(1-#/0.23)]&,{0,0.23}},LegendLabel->"ω", ColorFunctionScaling->False, LabelStyle->Directive[Bold, FontFamily->"Arial", FontSize->32], LegendMarkerSize->400]