Data Notebook for “Coexistence in Periodic Environments”

Model:

The productive (P) and less-productive versions of this model are separated as they will or can have unique parameter combinations reflecting trade-offs or periodic biology. We use the classical LV competition model:

f1P=x1(rP1(1-αP11 x1-αP12 x2 ));
f2P=x2(rP2(1-αP22 x2-αP21 x1));
f1LP=x1(rLP1(1-αLP11 x1-αLP12 x2 ));
f2LP=x2(rLP2(1-αLP22 x2-αLP21 x1 ));

Parameterization:

Parameter set for Figure 2b

paramP1=Dispatch[{
αP21->0.53,
αP12->0.25,
αP11->0.44,
αP22->0.66,
rP1->1.70,
rP2->1.20
}];
paramLP1=Dispatch[{
αLP21->0.48,
αLP12->1.83,
αLP11->1.11,
αLP22->0.83,
rLP1->0.3,
rLP2->1
}];
Seaslength1=Dispatch[{
p ->0.5
}];

Seasonally Mediate Coexistence: Parameter set for Figures 3a, 4, S2.1a, S2.2a, S3.1.1, S3.1.2, S4.1,

paramP2=Dispatch[{
αP21->0.53,
αP12->0.25,
αP11->0.5,
αP22->0.66,
rP1->1.70,
rP2->1.20
}];
paramLP2=Dispatch[{
αLP21->0.583,
αLP12->1.22,
αLP11->1.11,
αLP22->0.83,
rLP1->0.3,
rLP2->1
}];
Seaslength2=Dispatch[{
p ->0.5
}];

Seasonally Mediated Competitive Exclusion: Parameter set for Figures 3b, 5, S2.1b, S2.2b, S3.2.1, S3.2.2, S4.2,

paramP3=Dispatch[{
αP21->0.4,
αP12->0.3125,
αP11->0.625,
αP22->0.66,
rP1->4,
rP2->1.20
}];
paramLP3=Dispatch[{
αLP21->1.14,
αLP12->0.8,
αLP11->2,
αLP22->1.43,
rLP1->0.3,
rLP2->1
}];
Seaslength3=Dispatch[{
p ->0.5
}];

Seasonally Mediated Contingent Coexistence: Parameter set for Figures 3c, 6, S2.1ci-cii, S2.2ci-cii, S3.3.1, S3.3.2, S4.3

paramP4=Dispatch[{
αP21->0.83,
αP12->0.5,
αP11->0.5,
αP22->0.55,
rP1->1.7,
rP2->1.20
}];
paramLP4=Dispatch[{
αLP21->2,
αLP12->2.86,
αLP11->2.5,
αLP22->1,
rLP1->0.3,
rLP2->1
}];
Seaslength4=Dispatch[{
p ->0.5
}];

Isocline Approximation

Here, for each species, the productive period (P) is scaled by p and the less-productive period (LP) is scaled by (1-p). The two scaled periods are set equal and opposite to each other such that they “cancel” each other out and the net change over a full time period is 0.

Yearly Isocline for species 1

s1=Solve[p f1P==-(1-p)f1LP,x1];

ss1=s1[[All,1,2]];

ss1[[2]];

fy1=ss1[[2]];

x1y=x1(-x1+fy1);

Yearly Isocline for species 2

s2=Solve[p f2P==-(1-p)f2LP,x1];

ss2=s2[[All,1,2]];

ss2[[1]];

fy2=ss2[[1]];

x2y=x2(-x1+fy2);

Figure 2b)

Yearly isocline plot for Figure 2bi when p = 0.5

Plot[{fy1/.paramP1/.paramLP1/.p->0.5,
fy2/.paramP1/.paramLP1/.p->0.5}, {x2, -0.5, 3},
PlotRange->{{0,2.1},{0,2.1}},FrameLabel->{Style["x2 Density",FontSize->14,Bold],Style["x1 Density",FontSize->14,Bold]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline plot for Figure 2bii when p = 0.6

Plot[{fy1/.paramP1/.paramLP1/.p->0.6,
fy2/.paramP1/.paramLP1/.p->0.6}, {x2, -0.5, 3},
PlotRange->{{0,2.1},{0,2.1}},FrameLabel->{Style["x2 Density",FontSize->14,Bold],Style["x1 Density",FontSize->14,Bold]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline plot for Figure 2biii when p = 0.7

Plot[{fy1/.paramP1/.paramLP1/.p->0.7,
fy2/.paramP1/.paramLP1/.p->0.7}, {x2, -0.5, 3},
PlotRange->{{0,2.1},{0,2.1}},FrameLabel->{Style["x2 Density",FontSize->14,Bold],Style["x1 Density",FontSize->14,Bold]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Figure 3: Seasonally Mediated Outcomes

The full model is given with a set of initial conditions, to facilitate use with NDSolve we make the following functions:

LVcompP[{x10_,x20_}]:=Module[{transform={x1->x1[t],x2->x2[t]}},
{x1'[t]==(f1P/.transform),
x2'[t]==(f2P/.transform),
x1[0]==x10,
x2[0]==x20}
]

LVcompLP[{x10_,x20_}]:=Module[{transform={x1->x1[t],x2->x2[t]}},
{x1'[t]==(f1LP/.transform),
x2'[t]==(f2LP/.transform),
x1[0]==x10,
x2[0]==x20}
]

Figure 3a) Seasonally Mediated Coexistence

For the model experiments, we table over time and integrate in steps of 0-frs, frs-1 iterated to total run time, e.g.,do t=1,1000 with two forced step sizes that alter parameters that alter in calling above L-V productive and less-productive equations (i.e. different parameter settings).
Function, solSeas, when given the productive length as a fraction of the year, frs, initial values and stepsize, stitches together a productive solution with a less-productive solution {x1,x2} over 1 year in increments of stepsize years (e.g., 0.001).

solSeas1[fr_,stepSize_,x10_,x20_] :=
Module[{solyear},
soltestP1=NDSolve[LVcompP[{x10,x20}]/.paramP2,{x1,x2},{t,0,fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
soltestLP1=NDSolve[LVcompLP[{x1[fr]/.soltestP1,x2[fr]/.soltestP1}]/.paramLP2,{x1,x2},{t,0,1.0-fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
sol1=Join[Flatten[First[x1/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x1/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
sol2=Join[Flatten[First[x2/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x2/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
solyear=Transpose[{sol1,sol2}]
]

Function to stitch all periodic solutions together into longer term solutions. We have moved away from “Table” to “Do Loops” with conditionals as needed local memory of previous solution (0 to 1) for next starting value and “table” does not like this (no local memory in previous function call). 
Note, there’s a small error that’s repeatedly joining last 1 (enters as 0) so comes in twice, therefore, remove first item in list from holdadd in join and will be fine.
Note, the second “do loop” accounts for any productive lengths (p’s) that are smaller than 0.01. Therefore, when finding the mean equilibrium in wintershrink (see below), each experiment needs to be run for an extra 100 integers than the max number of integers used for the mean (e.g., to get the mean from integer 900-1000, experiment needs to be run for 1100 time steps).

StitchedSol1[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas1[fr,stpSz,0.25,0.25]];
If[i>0,holdadd=solSeas1[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

wintshrink1=Table[{frs,Mean[StitchedSol1[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol1[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

Show[{
ListLinePlot[{Transpose[{wintshrink1[[All,1]],wintshrink1[[All,2]]}],Transpose[{wintshrink1[[All,1]],wintshrink1[[All,3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"},
Epilog->{LightGray,Opacity[0.5],Rectangle[{0.0,0},{0.11,2.2}],Rectangle[{0.692,0},{1,2.1}],Black,Line[{{0.11,0},{0.11,2.2}}],Line[{{0.692,0},{0.692,2.2}}]},
FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400]}]

Figure 3b) Seasonally Mediated Competitive Exclusion

solSeas2[fr_,stepSize_,x10_,x20_] :=
Module[{solyear},
soltestP1=NDSolve[LVcompP[{x10,x20}]/.paramP3,{x1,x2},{t,0,fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
soltestLP1=NDSolve[LVcompLP[{x1[fr]/.soltestP1,x2[fr]/.soltestP1}]/.paramLP3,{x1,x2},{t,0,1.0-fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
sol1=Join[Flatten[First[x1/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x1/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
sol2=Join[Flatten[First[x2/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x2/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
solyear=Transpose[{sol1,sol2}]
]

StitchedSol2[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas2[fr,stpSz,0.25,0.25]];
If[i>0,holdadd=solSeas2[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

wintshrink2=Table[{frs,Mean[StitchedSol2[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol2[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

Show[{
ListLinePlot[{Transpose[{wintshrink2[[All,1]],wintshrink2[[All,2]]}],Transpose[{wintshrink2[[All,1]],wintshrink2[[All,3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"},
Epilog->{LightGray,Opacity[0.5],Rectangle[{0.18,0},{0.51,2.2}],Black,Line[{{0.18,0},{0.18,2.2}}],Line[{{0.51,0},{0.51,2.2}}]},
FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400]}]

Figure 3c) Seasonally Mediated Contingent Coexistence

solSeas3[fr_,stepSize_,x10_,x20_] :=
Module[{solyear},
soltestP1=NDSolve[LVcompP[{x10,x20}]/.paramP4,{x1,x2},{t,0,fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
soltestLP1=NDSolve[LVcompLP[{x1[fr]/.soltestP1,x2[fr]/.soltestP1}]/.paramLP4,{x1,x2},{t,0,1.0-fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
sol1=Join[Flatten[First[x1/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x1/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
sol2=Join[Flatten[First[x2/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x2/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
solyear=Transpose[{sol1,sol2}]
]

StitchedSol3[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas3[fr,stpSz,20,0.25]];
If[i>0,holdadd=solSeas3[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

StitchedSol4[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas3[fr,stpSz,0.25,20]];
If[i>0,holdadd=solSeas3[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

wintshrink3=Table[{frs,Mean[StitchedSol3[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol3[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

wintshrink4=Table[{frs,Mean[StitchedSol4[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol4[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

Show[{
ListLinePlot[{Transpose[{wintshrink3[[Range[1,7],1]],wintshrink3[[Range[1,7],2]]}],Transpose[{wintshrink3[[Range[1,7],1]],wintshrink3[[Range[1,7],3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},

Epilog->{LightGray,Opacity[0.5],Rectangle[{0,0},{0.065,2.1}],Rectangle[{0.68,0},{1,2.1}],
Gray,Opacity[0.3],Rectangle[{0.065,0},{0.68,2.1}],
Black,Line[{{0.065,0},{0.065,2.2}}],Line[{{0.68,0},{0.68,2.2}}]},

PlotStyle->{{Black,Dashed,Thickness[0.008]},{Black,Thickness[0.008]}},PlotLegends->{"x1","x2"},FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400],
ListLinePlot[{Transpose[{wintshrink3[[Range[8,101],1]],wintshrink3[[Range[8,101],2]]}],Transpose[{wintshrink3[[Range[8,101],1]],wintshrink3[[Range[8,101],3]]}]},Frame->True,PlotRange->{{0,1},{-0.1,2.1}},PlotStyle->{{Black,Dashed,Thickness[0.008]},{Black,Thickness[0.008]}},ImageSize->400],

ListLinePlot[{Transpose[{wintshrink4[[Range[1,69],1]],wintshrink4[[Range[1,69],2]]}],Transpose[{wintshrink4[[Range[1,69],1]],wintshrink4[[Range[1,69],3]]}]},Frame->True,PlotRange->{{0,1},{-0.1,2.1}},PlotStyle->{{Gray,Dashed,Thickness[0.005]},{Gray,Thickness[0.005]}},PlotLegends->{"x1","x2"},ImageSize->400],
ListLinePlot[{Transpose[{wintshrink4[[Range[70,101],1]],wintshrink4[[Range[70,101],2]]}],Transpose[{wintshrink4[[Range[70,101],1]],wintshrink4[[Range[70,101],3]]}]},Frame->True,PlotRange->{{0,1},{-0.1,2.1}},PlotStyle->{{Gray,Dashed,Thickness[0.005]},{Gray,Thickness[0.005]}},ImageSize->400]
}]

Two-Dimensional Bifurcation Diagrams of Seasonally Mediated Outcomes

Figure 4a: Seasonally Mediated Coexistence

Bif1[param1_,Seasonality_]:=Table[If[((p rP1+(1-p)rLP1)/(p rP1 αP11 + (1-p)rLP1 αLP11 ))<((p rP2+(1-p)rLP2)/(p rP2 αP21 + (1-p)rLP2 αLP21 ))/.paramP2/.paramLP2/.Seaslength2,If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP2/.paramLP2/.Seaslength2,0,1],If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP2/.paramLP2/.Seaslength2,2,3]],{ param1,5,0,-0.01},{ Seasonality,0,1,0.01}];

ArrayPlot[Bif1[rP1,p],DataRange->{{0,1},{0,5}},AspectRatio->1/2,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},
Epilog->{Black,Thin,Line[{{0,1.7},{1,1.7}}]},FrameLabel->{Style["rP1",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure 5a: Seasonally Mediated Competitive Exclusion

Bif2[param1_,Seasonality_]:=Table[If[((p rP1+(1-p)rLP1)/(p rP1 αP11 + (1-p)rLP1 αLP11 ))<((p rP2+(1-p)rLP2)/(p rP2 αP21 + (1-p)rLP2 αLP21 ))/.paramP3/.paramLP3/.Seaslength3,If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP3/.paramLP3/.Seaslength3,0,1],If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP3/.paramLP3/.Seaslength3,2,3]],{ param1,10,0,-0.01},{ Seasonality,0,1,0.01}];

ArrayPlot[Bif2[rP1,p],DataRange->{{0,1},{0,10}},AspectRatio->1/2,FrameTicks->{{Range[0,10,2],None},{Range[0,1,0.2],None}},
Epilog->{Black,Thin,Line[{{0,4},{1,4}}]},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rP1",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure 6a: Seasonally Mediated Contingent Coexistence

Bif3[param1_,Seasonality_]:=Table[If[((p rP1+(1-p)rLP1)/(p rP1 αP11 + (1-p)rLP1 αLP11 ))<((p rP2+(1-p)rLP2)/(p rP2 αP21 + (1-p)rLP2 αLP21 ))/.paramP4/.paramLP4/.Seaslength4,If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP4/.paramLP4/.Seaslength4,0,1],If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP4/.paramLP4/.Seaslength4,2,3]],{ param1,5,0,-0.01},{ Seasonality,0,1,0.01}];

ArrayPlot[Bif3[rLP2,p],DataRange->{{0,1},{0,5}},AspectRatio->1/2,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},
Epilog->{Black,Thin,Line[{{0,1},{1,1}}]},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rLP2",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Isocline Approximations for Bifurcation Diagrams of Seasonally Mediated Outcomes

Here, we calculate the seasonal equilibrium points (SNequ). The table generated for SNequ along with the different coloured circles, which represent either stable (ST) or unstable (UNST), will be used in the following isocline approximation plots.

equ=Solve[{x1y==0,x2y==0},{x1,x2}][[2]];

snequ1[sparam_]:=Table[{x2,x1}/.equ,{sparam,0,1,0.001}]

Seasonally Mediated Coexistence: Figures 4b and S4.1

SNequ=snequ1[p]/.paramP2/.paramLP2/.Seaslength2;

blueUNST={Blue,Thickness[0.3],Circle[]};
blueST={Blue,Disk[],Blue,Thick,Circle[]};
greenUNST={Green,Thickness[0.3],Circle[]};
greenST={Green,Disk[],Green,Thick,Circle[]};
yellowUNST={Yellow,Thickness[0.3],Circle[]};
yellowST={Yellow,Disk[],Yellow,Thick,Circle[]};
orangeUNST={Orange,Thickness[0.3],Circle[]};
orangeST={Orange,Disk[],Orange,Thick,Circle[]};
redUNST={Red,Thickness[0.3],Circle[]};
redST={Red,Disk[],Red,Thick,Circle[]};
blackUNST={Black,Thick,Circle[]};
blackST={Black,Disk[],Black,Thick,Circle[]};

Yearly isocline approximation plot for Figures 4bi and S4.1a, when p=0.1

Seaslength2=Dispatch[{
p ->0.1
}];

Plot[{fy1/.paramP2/.paramLP2/.Seaslength2,
fy2/.paramP2/.paramLP2/.Seaslength2}, {x2, -0.5, 3},
PlotRange->{{-0.5,2.0},{-0.5,2.0}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueUNST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blackUNST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeUNST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeUNST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redUNST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redUNST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0,1.14},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{1.25,0},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figure S4.1b, when p=0.115

Seaslength2=Dispatch[{
p ->0.115
}];

Plot[{fy1/.paramP2/.paramLP2/.Seaslength2,
fy2/.paramP2/.paramLP2/.Seaslength2}, {x2, -0.5, 3},
PlotRange->{{-0.5,2.0},{-0.5,2.0}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueUNST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueUNST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeUNST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeUNST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redUNST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redUNST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0,1.18},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{SNequ[[116]][[1]],SNequ[[116]][[2]]},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figures 4bii and S4.1c, when p=0.4

Seaslength2=Dispatch[{
p ->0.4
}];

Plot[{fy1/.paramP2/.paramLP2/.Seaslength2,
fy2/.paramP2/.paramLP2/.Seaslength2}, {x2, -0.5, 3},
PlotRange->{{-0.5,2.0},{-0.5,2.0}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueUNST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueUNST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeUNST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeUNST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redUNST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redUNST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0,1.6},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{1.33,0},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figure S4.1d, when p=0.692

Seaslength2=Dispatch[{
p ->0.692
}];

Plot[{fy1/.paramP2/.paramLP2/.Seaslength2,
fy2/.paramP2/.paramLP2/.Seaslength2}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueUNST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueUNST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeUNST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeUNST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redUNST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redUNST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{1.415,0},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blackUNST],{SNequ[[693]][[1]],SNequ[[693]][[2]]},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figures 4biii and S4.1d, when p=0.8

Seaslength2=Dispatch[{
p ->0.8
}];

Plot[{fy1/.paramP2/.paramLP2/.Seaslength2,
fy2/.paramP2/.paramLP2/.Seaslength2}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueUNST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueUNST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeUNST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeUNST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redUNST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redUNST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{0,1.9},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{1.45,0},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Seasonally Mediated Competitive Exclusion: Figures 5b and S4.2

SNequ=snequ1[p]/.paramP3/.paramLP3/.Seaslength3;

Yearly isocline approximation plot for Figures 5bi and S4.2a, when p=0.1

Seaslength3=Dispatch[{
p ->0.1
}];

Plot[{fy1/.paramP3/.paramLP3/.Seaslength3,
fy2/.paramP3/.paramLP3/.Seaslength3}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0,0.85},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0.75,0},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figure S4.2b, when p=0.1695

Seaslength3=Dispatch[{
p ->0.1695
}];

Plot[{fy1/.paramP3/.paramLP3/.Seaslength3,
fy2/.paramP3/.paramLP3/.Seaslength3}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0.78,0},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0,1.0055},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figures 5bii and S4.2c, when p=0.4

Seaslength3=Dispatch[{
p ->0.4
}];

Plot[{fy1/.paramP3/.paramLP3/.Seaslength3,
fy2/.paramP3/.paramLP3/.Seaslength3}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{0,1.31},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0.92,0},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figure S4.2d, when p=0.539

Seaslength3=Dispatch[{
p ->0.539
}];

Plot[{fy1/.paramP3/.paramLP3/.Seaslength3,
fy2/.paramP3/.paramLP3/.Seaslength3}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{1.025,0},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0,1.4126},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figures 5biii and S4.2e, when p=0.8

Seaslength3=Dispatch[{
p ->0.8
}];

Plot[{fy1/.paramP3/.paramLP3/.Seaslength3,
fy2/.paramP3/.paramLP3/.Seaslength3}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0,1.54},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{1.26,0},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Seasonally Mediated Contingent Coexistence: Figures 6b and S4.3

SNequ=snequ1[p]/.paramP4/.paramLP4/.Seaslength4;

Yearly isocline approximation plot for Figures 6bi and S4.3a, when p=0

Seaslength4=Dispatch[{
p ->0.0
}];

Plot[{fy1/.paramP4/.paramLP4/.Seaslength4,
fy2/.paramP4/.paramLP4/.Seaslength4}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueUNST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0,0.4},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{1,0},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figure S4.3b, when p=0.07

Seaslength4=Dispatch[{
p ->0.07
}];

Plot[{fy1/.paramP4/.paramLP4/.Seaslength4,
fy2/.paramP4/.paramLP4/.Seaslength4}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueUNST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{0,0.52},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{1.04,0},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figures 6bii and S4.3c, when p=0.4

Seaslength4=Dispatch[{
p ->0.4
}];

Plot[{fy1/.paramP4/.paramLP4/.Seaslength4,
fy2/.paramP4/.paramLP4/.Seaslength4}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueUNST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{0,1.09},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{1.25,0},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figure S4.3d, when p=0.692

Seaslength4=Dispatch[{
p ->0.692
}];

Plot[{fy1/.paramP4/.paramLP4/.Seaslength4,
fy2/.paramP4/.paramLP4/.Seaslength4}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueUNST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{0,1.55},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{1.4885,0},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Yearly isocline approximation plot for Figures 6biii and S4.3e, when p=0.8

Seaslength4=Dispatch[{
p ->0.8
}];

Plot[{fy1/.paramP4/.paramLP4/.Seaslength4,
fy2/.paramP4/.paramLP4/.Seaslength4}, {x2, -0.5, 3},
PlotRange->{{-0.5,2},{-0.5,2}},AspectRatio->1/1,
Epilog->{Black,Thick,Line[{{0,-1},{0,5}}],Line[{{-1,0},{4,0}}],
Inset[Graphics[blueST],{SNequ[[1]][[1]],SNequ[[1]][[2]]},{0,0},Scaled[{0.035,0.035}]],Inset[Graphics[blueUNST],{SNequ[[101]][[1]],SNequ[[101]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[201]][[1]],SNequ[[201]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[greenUNST],{SNequ[[301]][[1]],SNequ[[301]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[401]][[1]],SNequ[[401]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[501]][[1]],SNequ[[501]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[yellowUNST],{SNequ[[601]][[1]],SNequ[[601]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[701]][[1]],SNequ[[701]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[orangeST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[901]][[1]],SNequ[[901]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[redST],{SNequ[[1001]][[1]],SNequ[[1001]][[2]]},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{0,1.71},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackUNST],{1.6,0},{0,0},Scaled[{0.035,0.035}]],
Inset[Graphics[blackST],{SNequ[[801]][[1]],SNequ[[801]][[2]]},{0,0},Scaled[{0.035,0.035}]]},
FrameLabel->{Style["x2 Density",FontSize->14,Black],Style["x1 Density",FontSize->14,Black]},Frame->True ,Axes -> True,PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1","x2"}]

Supplementary Material for “Species Coexistence in Periodic Environments”

S2: Isocline Approximation and Numerical Simulation Comparison

Here, we calculate each or the equilibria that may be found within the isocline approximation

equ1=Solve[{x1y==0,x2y==0},{x1,x2}][[2]];
equ2=Solve[{x1y==0,x2y==0},{x1,x2}][[1]];
equ3=Solve[{x1y==0,x2y==0},{x1,x2}][[4]];

Figure S2.1a: Seasonally Mediated Coexistence

Show[{
ListLinePlot[{Transpose[{wintshrink1[[All,1]],wintshrink1[[All,2]]}],Transpose[{wintshrink1[[All,1]],wintshrink1[[All,3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},
Epilog->{LightGray,Opacity[0.5],Rectangle[{0.0,0},{0.11,2.2}],Rectangle[{0.692,0},{1,2.1}],Black,Line[{{0.11,0},{0.11,2.2}}],Line[{{0.692,0},{0.692,2.2}}]},
FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP2/.paramLP2,{p,0.115,0.692},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}],Plot[{x1,x2}/.equ2/.paramP2/.paramLP2,{p,0,0.11},PlotStyle->{Red,Dotted}],Plot[{x1,x2}/.equ3/.paramP2/.paramLP2,{p,0.692,1},PlotStyle->{Red,Dotted}]}]

Figure S2.1b: Seasonally Mediated Competitive Exclusion

Show[{
ListLinePlot[{Transpose[{wintshrink2[[All,1]],wintshrink2[[All,2]]}],Transpose[{wintshrink2[[All,1]],wintshrink2[[All,3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},
Epilog->{LightGray,Opacity[0.5],Rectangle[{0.18,0},{0.51,2.2}],Black,Line[{{0.18,0},{0.18,2.2}}],Line[{{0.51,0},{0.51,2.2}}]},
FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP3/.paramLP3,{p,0,0.18},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}],
Plot[{x1,x2}/.equ3/.paramP3/.paramLP3,{p,0.18,0.51},PlotStyle->{Red,Dotted}],
Plot[{x1,x2}/.equ1/.paramP3/.paramLP3,{p,0.51,1},PlotStyle->{Red,Dotted}]}]

Figure S2.1c: Seasonally Mediated Contingent Coexistence

Show[{
ListLinePlot[{Transpose[{wintshrink3[[Range[1,7],1]],wintshrink3[[Range[1,7],2]]}],Transpose[{wintshrink3[[Range[1,7],1]],wintshrink3[[Range[1,7],3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},
Epilog->{LightGray,Opacity[0.5],Rectangle[{0,0},{0.065,2.1}],Rectangle[{0.68,0},{1,2.1}],
Gray,Opacity[0.3],Rectangle[{0.065,0},{0.68,2.1}],
Black,Line[{{0.065,0},{0.065,2.2}}],Line[{{0.68,0},{0.68,2.2}}]},
PlotStyle->{{Black,Dashed,Thickness[0.008]},{Black,Thickness[0.008]}},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400],
ListLinePlot[{Transpose[{wintshrink3[[Range[8,101],1]],wintshrink3[[Range[8,101],2]]}],Transpose[{wintshrink3[[Range[8,101],1]],wintshrink3[[Range[8,101],3]]}]},Frame->True,PlotRange->{{0,1},{-0.1,2.1}},PlotStyle->{{Black,Dashed,Thickness[0.008]},{Black,Thickness[0.008]}},ImageSize->400],
ListLinePlot[{Transpose[{wintshrink4[[Range[1,69],1]],wintshrink4[[Range[1,69],2]]}],Transpose[{wintshrink4[[Range[1,69],1]],wintshrink4[[Range[1,69],3]]}]},Frame->True,PlotRange->{{0,1},{-0.1,2.1}},PlotStyle->{{Gray,Dashed,Thickness[0.005]},{Gray,Thickness[0.005]}},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},ImageSize->400],
ListLinePlot[{Transpose[{wintshrink4[[Range[70,101],1]],wintshrink4[[Range[70,101],2]]}],Transpose[{wintshrink4[[Range[70,101],1]],wintshrink4[[Range[70,101],3]]}]},Frame->True,PlotRange->{{0,1},{-0.1,2.1}},PlotStyle->{{Gray,Dashed,Thickness[0.005]},{Gray,Thickness[0.005]}},ImageSize->400],
Plot[{x1,x2}/.equ2/.paramP4/.paramLP4,{p,0,0.68},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate
Equilibrium"}],
Plot[{x1,x2}/.equ3/.paramP4/.paramLP4,{p,0.07,1},PlotStyle->{Red,Dotted}]
}]

Figure S2.1d: Counterintuitive Mean Density Change

Additional parameter combinations are used for more examples which are found in the supplementary material of “Coexistence in Periodic Environments”.

paramP5=Dispatch[{
αP21->0.4,
αP12->0.31,
αP11->0.63,
αP22->0.67,
rP1->1.70,
rP2->1.20
}];
paramLP5=Dispatch[{
αLP21->1.14,
αLP12->0.8,
αLP11->2,
αLP22->1.43,
rLP1->0.3,
rLP2->1
}];
Seaslength5=Dispatch[{
p ->0.5
}];

solSeas5[fr_,stepSize_,x10_,x20_] :=
Module[{solyear},
soltestP1=NDSolve[LVcompP[{x10,x20}]/.paramP5,{x1,x2},{t,0,fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
soltestLP1=NDSolve[LVcompLP[{x1[fr]/.soltestP1,x2[fr]/.soltestP1}]/.paramLP5,{x1,x2},{t,0,1.0-fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
sol1=Join[Flatten[First[x1/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x1/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
sol2=Join[Flatten[First[x2/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x2/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
solyear=Transpose[{sol1,sol2}]
]

StitchedSol5[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas5[fr,stpSz,0.25,0.25]];
If[i>0,holdadd=solSeas5[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

wintshrink5=Table[{frs,Mean[StitchedSol5[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol5[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

Show[{
ListLinePlot[{Transpose[{wintshrink5[[All,1]],wintshrink5[[All,2]]}],Transpose[{wintshrink5[[All,1]],wintshrink5[[All,3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},
FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP5/.paramLP5,{p,0,1},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]

Figure S2.1e: Mean Density Change with Competitive Exclusion in the Less Productive Season

paramP6=Dispatch[{
αP21->0.48,
αP12->0.25,
αP11->0.56,
αP22->0.67,
rP1->1.70,
rP2->1.20
}];
paramLP6=Dispatch[{
αLP21->0.59,
αLP12->1.25,
αLP11->1.11,
αLP22->0.77,
rLP1->0.3,
rLP2->1
}];
Seaslength6=Dispatch[{
p ->0.5
}];

solSeas6[fr_,stepSize_,x10_,x20_] :=
Module[{solyear},
soltestP1=NDSolve[LVcompP[{x10,x20}]/.paramP6,{x1,x2},{t,0,fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
soltestLP1=NDSolve[LVcompLP[{x1[fr]/.soltestP1,x2[fr]/.soltestP1}]/.paramLP6,{x1,x2},{t,0,1.0-fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
sol1=Join[Flatten[First[x1/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x1/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
sol2=Join[Flatten[First[x2/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x2/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
solyear=Transpose[{sol1,sol2}]
]

StitchedSol6[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas6[fr,stpSz,0.25,0.25]];
If[i>0,holdadd=solSeas6[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

wintshrink6=Table[{frs,Mean[StitchedSol6[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol6[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

Show[{
ListLinePlot[{Transpose[{wintshrink6[[All,1]],wintshrink6[[All,2]]}],Transpose[{wintshrink6[[All,1]],wintshrink6[[All,3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},
Epilog->{LightGray,Opacity[0.5],Rectangle[{0.0,0},{0.15,2.2}],Black,Line[{{0.15,0},{0.15,2.2}}]},
FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP6/.paramLP6,{p,0.15,1},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}],
Plot[{x1,x2}/.equ2/.paramP6/.paramLP6,{p,0,0.15},PlotStyle->{Red,Dotted}]}]

Figure S2.1f: Mean Density Change with Competitive Exclusion in the Productive Season

paramP7=Dispatch[{
αP21->0.53,
αP12->0.25,
αP11->0.5,
αP22->0.67,
rP1->1.70,
rP2->1.20
}];
paramLP7=Dispatch[{
αLP21->0.43,
αLP12->1.25,
αLP11->1.11,
αLP22->1.43,
rLP1->0.3,
rLP2->1
}];
Seaslength7=Dispatch[{
p ->0.5
}];

solSeas7[fr_,stepSize_,x10_,x20_] :=
Module[{solyear},
soltestP1=NDSolve[LVcompP[{x10,x20}]/.paramP7,{x1,x2},{t,0,fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
soltestLP1=NDSolve[LVcompLP[{x1[fr]/.soltestP1,x2[fr]/.soltestP1}]/.paramLP7,{x1,x2},{t,0,1.0-fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
sol1=Join[Flatten[First[x1/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x1/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
sol2=Join[Flatten[First[x2/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x2/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
solyear=Transpose[{sol1,sol2}]
]

StitchedSol7[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas7[fr,stpSz,0.25,0.25]];
If[i>0,holdadd=solSeas7[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

wintshrink7=Table[{frs,Mean[StitchedSol7[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol7[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

Show[{
ListLinePlot[{Transpose[{wintshrink7[[All,1]],wintshrink7[[All,2]]}],Transpose[{wintshrink7[[All,1]],wintshrink7[[All,3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},
Epilog->{LightGray,Opacity[0.5],Rectangle[{0.86,0},{1,2.1}],Black,Line[{{0.86,0},{0.86,2.2}}]},
FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP7/.paramLP7,{p,0,0.86},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}],
Plot[{x1,x2}/.equ3/.paramP7/.paramLP7,{p,0.86,1},PlotStyle->{Red,Dotted}]}]

S2: Time Series at p=0.5 for the Comparison of Isocline Approximation against Numerical Simulations

Time series plot for Figure S2.2a: Seasonally Mediated Coexistence

SolFinal1=StitchedSol1[0.5,100,.01];

Show[{ListLinePlot[{Transpose[{SolFinal1[[All,1]],SolFinal1[[All,2]]}],Transpose[{SolFinal1[[All,1]],SolFinal1[[All,3]]}]},PlotRange->{{0,100},{-0.02,2.1}},PlotStyle->{{Gray,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},Frame->True,FrameLabel->{Style["Time",FontSize->14,Bold,Black],Style["Density",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP2/.paramLP2/.p->0.5,{x1,0,100},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]

Time series plot for Figure S2.2b: Seasonally Mediated Competitive Exclusion

SolFinal2=StitchedSol2[0.5,100,.01];

Show[{ListLinePlot[{Transpose[{SolFinal2[[All,1]],SolFinal2[[All,2]]}],Transpose[{SolFinal2[[All,1]],SolFinal2[[All,3]]}]},PlotRange->{{0,100},{-0.02,2.1}},PlotStyle->{{Gray,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},Frame->True,FrameLabel->{Style["Time",FontSize->14,Bold,Black],Style["Density",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ3/.paramP3/.paramLP3/.p->0.5,{x1,0,100},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]

Time series plot for Figure S2.2ci: Seasonally Mediated Contingent Coexistence

SolFinal3=StitchedSol3[0.5,100,.01];

Show[{ListLinePlot[{Transpose[{SolFinal3[[All,1]],SolFinal3[[All,2]]}],Transpose[{SolFinal3[[All,1]],SolFinal3[[All,3]]}]},PlotRange->{{0,100},{-0.02,2.1}},PlotStyle->{{Gray,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},Frame->True,FrameLabel->{Style["Time",FontSize->14,Bold,Black],Style["Density",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ3/.paramP4/.paramLP4/.p->0.5,{x1,0,100},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]

Time series plot for Figure S2.2cii: Seasonally Mediated Contingent Coexistence

SolFinal4=StitchedSol4[0.5,100,.01];

Show[{ListLinePlot[{Transpose[{SolFinal4[[All,1]],SolFinal4[[All,2]]}],Transpose[{SolFinal4[[All,1]],SolFinal4[[All,3]]}]},PlotRange->{{0,100},{-0.02,2.1}},PlotStyle->{{Gray,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},Frame->True,FrameLabel->{Style["Time",FontSize->14,Bold,Black],Style["Density",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ2/.paramP4/.paramLP4/.p->0.5,{x1,0,100},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]

Time series plot for Figure S2.2d: Counterintuitive Mean Density Change

SolFinal5=StitchedSol5[0.5,100,.01];

Show[{ListLinePlot[{Transpose[{SolFinal5[[All,1]],SolFinal5[[All,2]]}],Transpose[{SolFinal5[[All,1]],SolFinal5[[All,3]]}]},PlotRange->{{0,100},{-0.02,2.1}},PlotStyle->{{Gray,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},Frame->True,FrameLabel->{Style["Time",FontSize->14,Bold,Black],Style["Density",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP5/.paramLP5/.p->0.5,{x1,0,100},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]

Time series plot for Figure S2.2e: Mean Density Change with Competitive Exclusion in the Less Productive Season

SolFinal6=StitchedSol6[0.5,100,.01];

Show[{ListLinePlot[{Transpose[{SolFinal6[[All,1]],SolFinal6[[All,2]]}],Transpose[{SolFinal6[[All,1]],SolFinal6[[All,3]]}]},PlotRange->{{0,100},{-0.02,2.1}},PlotStyle->{{Gray,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},Frame->True,FrameLabel->{Style["Time",FontSize->14,Bold,Black],Style["Density",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP6/.paramLP6/.p->0.5,{x1,0,100},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]

Time series plot for Figure S2.2e: Mean Density Change with Competitive Exclusion in the Productive Season

SolFinal7=StitchedSol7[0.5,100,.01];

Show[{ListLinePlot[{Transpose[{SolFinal7[[All,1]],SolFinal7[[All,2]]}],Transpose[{SolFinal7[[All,1]],SolFinal7[[All,3]]}]},PlotRange->{{0,100},{-0.02,2.1}},PlotStyle->{{Gray,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},Frame->True,FrameLabel->{Style["Time",FontSize->14,Bold,Black],Style["Density",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP7/.paramLP7/.p->0.5,{x1,0,100},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]

S3: Generality of Results for all Seasonally-Mediated Competitive Outcomes Using Isocline Approach

Seasonally Mediate Coexistence Bifurcation Diagrams across other Growth Rate Parameters

Bif1[param1_,Seasonality_]:=Table[If[((p rP1+(1-p)rLP1)/(p rP1 αP11 + (1-p)rLP1 αLP11 ))<((p rP2+(1-p)rLP2)/(p rP2 αP21 + (1-p)rLP2 αLP21 ))/.paramP2/.paramLP2/.Seaslength2,If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP2/.paramLP2/.Seaslength2,0,1],If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP2/.paramLP2/.Seaslength2,2,3]],{ param1,5,0,-0.01},{ Seasonality,0,1,0.01}];

Figure S3.1.1a Species 2’s productive growth rate over the length of the productive period

ArrayPlot[Bif1[rP2,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rP2",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.1.1b Species 1’s less productive growth rate over the length of the productive period

ArrayPlot[Bif1[rLP1,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rLP1",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.1.1c Species 2’s less productive growth rate over the length of the productive period

ArrayPlot[Bif1[rLP2,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rLP2",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Seasonally Mediate Coexistence Bifurcation Diagrams for Inter- and Intraspecific Competition

Figure S3.1.2a Species 1’s productive intraspecific competitive coefficient

ArrayPlot[Bif1[αP11,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP11",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.1.2b Species 2’s productive intraspecific competitive coefficient

ArrayPlot[Bif1[αP22,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP22",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.1.2c Productive interspecific competitive coefficient of species 2 on species 1

ArrayPlot[Bif1[αP12,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP12",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.1.2d Productive interspecific competitive coefficient of species 1 on species 2

ArrayPlot[Bif1[αP21,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP21",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.1.2e Species 1’s less productive intraspecific competitive coefficient

ArrayPlot[Bif1[αLP11,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP11",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.1.2f Species 2’s less productive intraspecific competitive coefficient

ArrayPlot[Bif1[αLP22,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP22",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.1.2g Less productive interspecific competitive coefficient of species 2 on species 1

ArrayPlot[Bif1[αLP12,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP12",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.1.2h Less productive interspecific competitive coefficient of species 1 on species 2

ArrayPlot[Bif1[αLP21,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP21",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Seasonally Mediate Competitive Exclusion Bifurcation Diagrams across other Growth Rate Parameters

Bif2[param1_,Seasonality_]:=Table[If[((p rP1+(1-p)rLP1)/(p rP1 αP11 + (1-p)rLP1 αLP11 ))<((p rP2+(1-p)rLP2)/(p rP2 αP21 + (1-p)rLP2 αLP21 ))/.paramP3/.paramLP3/.Seaslength3,If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP3/.paramLP3/.Seaslength3,0,1],If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP3/.paramLP3/.Seaslength3,2,3]],{ param1,5,0,-0.01},{ Seasonality,0,1,0.01}];

FigureS3.2.1a Species 2’s productive growth rate

ArrayPlot[Bif2[rP2,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rP2",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

FigureS3.2.1b Species 1’s less productive growth rate

ArrayPlot[Bif2[rLP1,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rLP1",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

FigureS3.2.1c Species 2’s less productive growth rate

ArrayPlot[Bif2[rLP2,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rLP2",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Seasonally Mediate Competitive Exclusion Bifurcation Diagrams for Inter- and Intraspecific Competition

Figure S3.2.2a Species 1’s productive intraspecific competitive coefficient

ArrayPlot[Bif2[αP11,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP11",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.2.2b Species 2’s productive intraspecific competitive coefficient

ArrayPlot[Bif2[αP22,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP22",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.2.2c Productive interspecific competitive coefficient of species 2 on species 1

ArrayPlot[Bif2[αP12,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP12",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.2.2d Productive interspecific competitive coefficient of species 1 on species 2

ArrayPlot[Bif2[αP21,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP21",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.2.2e Species 1’s less productive intraspecific competitive coefficient

ArrayPlot[Bif2[αLP11,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP11",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.2.2f Species 2’s less productive intraspecific competitive coefficient

ArrayPlot[Bif2[αLP22,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP22",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.2.2g Less productive interspecific competitive coefficient of species 2 on species 1

ArrayPlot[Bif2[αLP12,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP12",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.2.2h Less productive interspecific competitive coefficient of species 1 on species 2

ArrayPlot[Bif2[αLP21,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP21",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Seasonally Mediate Contingent Coexistence Bifurcation Diagrams across other Growth Rate Parameters

Bif3[param1_,Seasonality_]:=Table[If[((p rP1+(1-p)rLP1)/(p rP1 αP11 + (1-p)rLP1 αLP11 ))<((p rP2+(1-p)rLP2)/(p rP2 αP21 + (1-p)rLP2 αLP21 ))/.paramP4/.paramLP4/.Seaslength4,If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP4/.paramLP4/.Seaslength4,0,1],If[((p rP2+(1-p)rLP2)/(p rP2 αP22 + (1-p)rLP2 αLP22))<((p rP1+(1-p)rLP1)/(p rP1 αP12+ (1-p)rLP1 αLP12 ))/.paramP4/.paramLP4/.Seaslength4,2,3]],{ param1,5,0,-0.01},{ Seasonality,0,1,0.01}];

FigureS3.3.1a Species 1’s productive growth rate

ArrayPlot[Bif3[rP1,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rP1",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

FigureS3.3.1b Species 2’s productive growth rate

ArrayPlot[Bif3[rP2,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rP2",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

FigureS3.3.1c Species 1’s less productive growth rate

ArrayPlot[Bif3[rLP1,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["rLP1",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Seasonally Mediate Contingent Coexistence Bifurcation Diagrams for Inter- and Intraspecific Competition

Figure S3.3.2a Species 1’s productive intraspecific competitive coefficient

ArrayPlot[Bif3[αP11,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP11",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.3.2b Species 2’s productive intraspecific competitive coefficient

ArrayPlot[Bif3[αP22,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP22",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.3.2c Productive interspecific competitive coefficient of species 2 on species 1

ArrayPlot[Bif3[αP12,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP12",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.3.2d Productive interspecific competitive coefficient of species 1 on species 2

ArrayPlot[Bif3[αP21,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αP21",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.3.2e Species 1’s less productive intraspecific competitive coefficient

ArrayPlot[Bif3[αLP11,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP11",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.3.2f Species 2’s less productive intraspecific competitive coefficient

ArrayPlot[Bif3[αLP22,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP22",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.3.2g Less productive interspecific competitive coefficient of species 2 on species 1

ArrayPlot[Bif3[αLP12,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP12",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

Figure S3.3.2h Less productive interspecific competitive coefficient of species 1 on species 2

ArrayPlot[Bif3[αLP21,p],DataRange->{{0,1},{0,5}},AspectRatio->1/1,FrameTicks->{{Range[0,5,1],None},{Range[0,1,0.2],None}},ColorRules->{0->White,1->Lighter[LightGray],2->Lighter[LightGray],3->LightGray},FrameLabel->{Style["αLP21",FontSize->14,Black,Bold],Style["Length of Productive Period (p)",FontSize->14,Black,Bold]}]

S5: Approximation Breakdown

Figure S5ai and S5aii Inaccuracy of Seasonally Mediated Competitive Exclusion

paramP8=Dispatch[{
αP21->0.4,
αP12->0.3125,
αP11->0.625,
αP22->0.66,
rP1->40,
rP2->12
}];
paramLP8=Dispatch[{
αLP21->1.14,
αLP12->0.8,
αLP11->2,
αLP22->1.43,
rLP1->3,
rLP2->10
}];
Seaslength8=Dispatch[{
p ->0.5
}];

solSeas8[fr_,stepSize_,x10_,x20_] :=
Module[{solyear},
soltestP1=NDSolve[LVcompP[{x10,x20}]/.paramP8,{x1,x2},{t,0,fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
soltestLP1=NDSolve[LVcompLP[{x1[fr]/.soltestP1,x2[fr]/.soltestP1}]/.paramLP8,{x1,x2},{t,0,1.0-fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
sol1=Join[Flatten[First[x1/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x1/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
sol2=Join[Flatten[First[x2/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x2/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
solyear=Transpose[{sol1,sol2}]
]

StitchedSol8[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas8[fr,stpSz,0.25,0.25]];
If[i>0,holdadd=solSeas8[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

wintshrink8=Table[{frs,Mean[StitchedSol8[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol8[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

Show[{
ListLinePlot[{Transpose[{wintshrink8[[All,1]],wintshrink8[[All,2]]}],Transpose[{wintshrink8[[All,1]],wintshrink8[[All,3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},
FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP8/.paramLP8,{p,0,0.18},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}],
Plot[{x1,x2}/.equ3/.paramP8/.paramLP8,{p,0.18,0.54},PlotStyle->{Red,Dotted}],
Plot[{x1,x2}/.equ1/.paramP8/.paramLP8,{p,0.54,1},PlotStyle->{Red,Dotted}]}]

SolFinal8=StitchedSol8[0.5,100,.01];

meanval8=Table[{frs,Mean[StitchedSol8[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol8[frs,1100,.01][[90000;;100000,3]]]},{frs,0.5,0.5,.01}];

Show[{ListLinePlot[{Transpose[{SolFinal8[[All,1]],SolFinal8[[All,2]]}],Transpose[{SolFinal8[[All,1]],SolFinal8[[All,3]]}]},PlotRange->{{96,100},{-0.02,2.1}},PlotStyle->{{Gray,Dashed},Gray},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},Frame->True,Epilog->{LightGray,Opacity[0.4],Rectangle[{96.5,0},{97,2.2}],Rectangle[{97.5,0},{98,2.2}],Rectangle[{98.5,0},{99,2.2}],Rectangle[{99.5,0},{100,2.2}]}, FrameLabel->{Style["Time",FontSize->14,Bold,Black],Style["Density",FontSize->14,Bold,Black]},ImageSize->400],
Plot[x1=meanval8[[1]][[2]],{x1,95,105},PlotStyle->{Black,Dashed},PlotLegends->{"x1 Numerical
     Mean"}],Plot[x1=meanval8[[1]][[3]],{x1,95,105},PlotStyle->Black,PlotLegends->{"x2 Numerical 
     Mean"}],Plot[{x1,x2}/.equ3/.paramP8/.paramLP8/.p->0.5,{x1,0,100},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]

Figure S5bi and S5bii Seasonally Mediated Coexistence

paramP9=Dispatch[{
αP21->0.53,
αP12->0.25,
αP11->0.5,
αP22->0.66,
rP1->17,
rP2->12
}];
paramLP9=Dispatch[{
αLP21->0.583,
αLP12->1.22,
αLP11->1.11,
αLP22->0.83,
rLP1->3,
rLP2->10
}];
Seaslength9=Dispatch[{
p ->0.5
}];

solSeas9[fr_,stepSize_,x10_,x20_] :=
Module[{solyear},
soltestP1=NDSolve[LVcompP[{x10,x20}]/.paramP9,{x1,x2},{t,0,fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
soltestLP1=NDSolve[LVcompLP[{x1[fr]/.soltestP1,x2[fr]/.soltestP1}]/.paramLP9,{x1,x2},{t,0,1.0-fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
sol1=Join[Flatten[First[x1/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x1/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
sol2=Join[Flatten[First[x2/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x2/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
solyear=Transpose[{sol1,sol2}]
]

StitchedSol9[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas9[fr,stpSz,0.25,0.25]];
If[i>0,holdadd=solSeas9[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

wintshrink9=Table[{frs,Mean[StitchedSol9[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol9[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

Show[{
ListLinePlot[{Transpose[{wintshrink9[[All,1]],wintshrink9[[All,2]]}],Transpose[{wintshrink9[[All,1]],wintshrink9[[All,3]]}]},Frame->True,PlotRange->{{0,1},{-0.02,2.1}},PlotStyle->{{Black,Dashed},Black},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},
Epilog->{LightGray,Opacity[0.5],Rectangle[{0.0,0},{0.11,2.2}],Rectangle[{0.79,0},{1,2.1}],Black,Line[{{0.11,0},{0.11,2.2}}],Line[{{0.79,0},{0.79,2.2}}]},
FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400],
Plot[{x1,x2}/.equ1/.paramP9/.paramLP9,{p,0.116,0.692},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}],Plot[{x1,x2}/.equ2/.paramP9/.paramLP9,{p,0,0.116},PlotStyle->{Red,Dotted}],Plot[{x1,x2}/.equ3/.paramP2/.paramLP2,{p,0.692,1},PlotStyle->{Red,Dotted}]}]

SolFinal9=StitchedSol9[0.5,100,.01];

meanval9=Table[{frs,Mean[StitchedSol9[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol9[frs,1100,.01][[90000;;100000,3]]]},{frs,0.5,0.5,.01}];

Show[{ListLinePlot[{Transpose[{SolFinal9[[All,1]],SolFinal9[[All,2]]}],Transpose[{SolFinal9[[All,1]],SolFinal9[[All,3]]}]},PlotRange->{{96,100},{-0.02,2.1}},PlotStyle->{{Gray,Dashed},Gray},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},Frame->True,Epilog->{LightGray,Opacity[0.4],Rectangle[{96.5,0},{97,2.2}],Rectangle[{97.5,0},{98,2.2}],Rectangle[{98.5,0},{99,2.2}],Rectangle[{99.5,0},{100,2.2}]}, FrameLabel->{Style["Time",FontSize->14,Bold,Black],Style["Density",FontSize->14,Bold,Black]},ImageSize->400],
Plot[x1=meanval9[[1]][[2]],{x1,95,105},PlotStyle->{Black,Dashed},PlotLegends->{"x1 Numerical
     Mean"}],Plot[x1=meanval9[[1]][[3]],{x1,95,105},PlotStyle->Black,PlotLegends->{"x2 Numerical 
     Mean"}],Plot[{x1,x2}/.equ1/.paramP9/.paramLP9/.p->0.5,{x1,0,100},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]

Figure S5ci and S5cii Seasonally Mediated Contingent Coexistence

paramP10=Dispatch[{
αP21->0.83,
αP12->0.5,
αP11->0.5,
αP22->0.55,
rP1->17,
rP2->12
}];
paramLP10=Dispatch[{
αLP21->2,
αLP12->2.86,
αLP11->2.5,
αLP22->1,
rLP1->3,
rLP2->10
}];
Seaslength10=Dispatch[{
p ->0.5
}];

solSeas10[fr_,stepSize_,x10_,x20_] :=
Module[{solyear},
soltestP1=NDSolve[LVcompP[{x10,x20}]/.paramP10,{x1,x2},{t,0,fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
soltestLP1=NDSolve[LVcompLP[{x1[fr]/.soltestP1,x2[fr]/.soltestP1}]/.paramLP10,{x1,x2},{t,0,1.0-fr},AccuracyGoal->70,MaxStepSize->0.05,MaxSteps->500000];
sol1=Join[Flatten[First[x1/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x1/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
sol2=Join[Flatten[First[x2/.soltestP1][Range[0.0, fr, stepSize]]],Flatten[First[x2/.soltestLP1][Range[stepSize, 1.0-fr, stepSize]]]];
solyear=Transpose[{sol1,sol2}]
]

StitchedSol10[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas10[fr,stpSz,20,0.25]];
If[i>0,holdadd=solSeas10[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

StitchedSol11[fr_,tnd_,stpSz_]:=
Module[{StitSol},
Do[
If[i==0,holdadd=solSeas10[fr,stpSz,0.25,20]];
If[i>0,holdadd=solSeas10[fr,stpSz,hold[[-1,1]],hold[[-1,2]]]];
If[i==0,newhold=holdadd,newhold=Join[newhold,Delete[holdadd,1]]];
hold=holdadd,
{i,0,tnd}];
Do[
If[Length[hold]==100,solTime=Table[i,{i,0,tnd+(0.99-0.01(tnd)),.01}]];
If[Length[hold]==101,solTime=Table[i,{i,0,tnd+1,.01}]];
solTime,
{i,0,tnd}];
TimeSol=Transpose[Insert[Transpose[newhold],solTime,1]];
StitSol=TimeSol
]

wintshrink10=Table[{frs,Mean[StitchedSol10[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol10[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

wintshrink11=Table[{frs,Mean[StitchedSol11[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol11[frs,1100,.01][[90000;;100000,3]]]},{frs,0.0,1,.01}];

Show[{
ListLinePlot[{Transpose[{wintshrink10[[Range[1,8],1]],wintshrink10[[Range[1,8],2]]}],Transpose[{wintshrink10[[Range[1,8],1]],wintshrink10[[Range[1,8],3]]}]},Frame->True,
PlotRange->{{0,1},{-0.2,2.1}},
Epilog->{LightGray,Opacity[0.5],Rectangle[{0,0},{0.075,2.1}],Rectangle[{0.735,0},{1,2.1}],
Gray,Opacity[0.3],Rectangle[{0.075,0},{0.735,2.1}],
Black,Line[{{0.075,0},{0.075,2.2}}],Line[{{0.735,0},{0.735,2.2}}]},
PlotStyle->{{Black,Dashed,Thickness[0.008]},{Black,Thickness[0.008]}},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},FrameLabel->{Style["Length of Productive Period (p)",FontSize->14,Bold,Black],Style["Mean Densities",FontSize->14,Bold,Black]},ImageSize->400],
ListLinePlot[{Transpose[{wintshrink10[[Range[9,101],1]],wintshrink10[[Range[9,101],2]]}],Transpose[{wintshrink10[[Range[9,101],1]],wintshrink10[[Range[9,101],3]]}]},Frame->True,PlotRange->{{0,1},{-0.1,2.1}},PlotStyle->{{Black,Dashed,Thickness[0.008]},{Black,Thickness[0.008]}},ImageSize->400],
ListLinePlot[{Transpose[{wintshrink11[[Range[1,74],1]],wintshrink11[[Range[1,74],2]]}],Transpose[{wintshrink11[[Range[1,74],1]],wintshrink11[[Range[1,74],3]]}]},Frame->True,PlotRange->{{0,1},{-0.1,2.1}},PlotStyle->{{Gray,Dashed,Thickness[0.005]},{Gray,Thickness[0.005]}},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},ImageSize->400],
ListLinePlot[{Transpose[{wintshrink11[[Range[75,101],1]],wintshrink11[[Range[75,101],2]]}],Transpose[{wintshrink11[[Range[75,101],1]],wintshrink11[[Range[75,101],3]]}]},Frame->True,PlotRange->{{0,1},{-0.1,2.1}},PlotStyle->{{Gray,Dashed,Thickness[0.005]},{Gray,Thickness[0.005]}},ImageSize->400],
Plot[{x1,x2}/.equ2/.paramP10/.paramLP10,{p,0,0.735},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate
Equilibrium"}],
Plot[{x1,x2}/.equ3/.paramP10/.paramLP10,{p,0.075,1},PlotStyle->{Red,Dotted}]
}]

SolFinal11=StitchedSol11[0.7,100,.01];

meanval11=Table[{frs,Mean[StitchedSol11[frs,1100,.01][[90000;;100000,2]]],Mean[StitchedSol11[frs,1100,.01][[90000;;100000,3]]]},{frs,0.7,0.7,.01}];

Show[{ListLinePlot[{Transpose[{SolFinal11[[All,1]],SolFinal11[[All,2]]}],Transpose[{SolFinal11[[All,1]],SolFinal11[[All,3]]}]},PlotRange->{{96,100},{-0.02,2.1}},PlotStyle->{{Gray,Dashed},{Gray}},PlotLegends->{"x1 Numerical 
     Simulation","x2 Numerical 
     Simulation"},Frame->True,Epilog->{LightGray,Opacity[0.4],Rectangle[{96.7,0},{97,2.2}],Rectangle[{97.7,0},{98,2.2}],Rectangle[{98.7,0},{99,2.2}],Rectangle[{99.7,0},{100,2.2}]}, FrameLabel->{Style["Time",FontSize->14,Bold,Black],Style["Density",FontSize->14,Bold,Black]},ImageSize->400],
Plot[x1=meanval11[[1]][[2]],{x1,95,105},PlotStyle->{Black, Dashed},PlotLegends->{"x1 Numerical 
     Mean"}],Plot[x1=meanval11[[1]][[3]],{x1,95,105},PlotStyle->{Black},PlotLegends->{"x2 Numerical 
     Mean"}],Plot[{x1,x2}/.equ2/.paramP10/.paramLP10/.p->0.7,{x1,0,100},PlotStyle->{Red,Dotted},PlotLegends->{"Approximate 
Equilibrium"}]}]