Supplementary Mathematica Notebook for “Habitat distribution affects the connectivity and population size in migratory networks”

Swati Patel (Mathematics Department, Oregon State University, patelswa@oregonstate.edu) and Caz Taylor (Biology, Tulane University)

This notebook provides the symbolic computational analysis used in identifying the invasion criteria for all relevant equilibrium points of our migratory species population model without dispersal (see main text).  This notebook follows along with the Supplementary Information and the Supplementary Figure 1.

Model:

f1=(1-c11*(x12+x11)/Kb1-c11*(x21+x11)/Kw1);
f2=(1-c12*(x12+x11)/Kb1-c12*(x22+x12)/Kw2);
f3=(1-c21*(x21+x22)/Kb2-c21*(x21+x11)/Kw1);
f4=(1-c22*(x21+x22)/Kb2-c22*(x22+x12)/Kw2);

x11dot= x11*r11*f1;
x12dot= x12*r12*f2;
x21dot= x21*r21*f3;
x22dot=x22*r22 * f4;

nondimsub = {Kb1-> γ1*s*K, Kb2-> γ2*(1-s)*K, Kw1-> (1-γ1)*s*K, Kw2-> (1-γ2)*(1-s)*K}

Equilibria:

SOLVE FOR SINGLE POPULATION EQUILIBRIA:

EQA=Solve[{f1/.x21->0/.x12->0}==0, {x11}][[1]]

EQB=Solve[{f4/.x21->0/.x12->0}==0, {x22}][[1]]

EQC=Solve[{f2/.x22->0/.x11->0}==0, {x12}][[1]]

EQD=Solve[{f3/.x22->0/.x11->0}==0, {x21}][[1]]

SOLVE FOR TWO POPULATION EQUILIBRIA:

EQE=Solve[{f1==0, f4==0}/.x12->0/.x21->0, {x11,x22}][[1]]

EQF=Solve[{f1==0, f2==0}/.x22->0/.x21->0, {x11,x12}][[1]]

EQG=Solve[{f1==0, f3==0}/.x12->0/.x22->0, {x11,x21}][[1]]

EQH=Solve[{f3==0, f4==0}/.x12->0/.x11->0, {x21,x22}][[1]]

EQI=Solve[{f2==0, f4==0}/.x11->0/.x21->0, {x12,x22}][[1]]

EQJ=Solve[{f2== 0, f3==0}/.x11->0/.x22->0, {x21,x12}] [[1]]

SOLVE FOR THREE POPULATION EQUILIBRIA:

EQK=Solve[{f1==0,f4 == 0, f3==0}/.x12->0, {x11,x21,x22}] [[1]]

EQL=Solve[{f1==0,f2== 0, f4==0}/.x21->0, {x11,x12,x22}] [[1]]

EQM=Solve[{f1==0,f2== 0, f3==0}/.x22->0, {x11,x21,x12}] [[1]]

EQN=Solve[{f4==0,f2== 0, f3==0}/.x11->0, {x22,x21,x12}] [[1]]

SOLVE FOR FOUR POPULATIONS:

Solve[{f4==0,f2== 0, f3==0, f1==0}, {x22,x21,x12, x11}]

Equilibria:

CONDITIONS FOR INVASION of second population into equilibrium with one population:
	notation: INV2A is the condition for INVasion of population 2 into equilibrium given by A.

INV2A = f2/.EQA/.x12->0/.x22->0

Does INV2A imply EQF is positive and vice versa? The following shows that if INV2A is positive then there exists a positive equilibrium and vice versa, assuming c11=c22=1.

FindInstance[{INV2A>0,x11>0, x12>0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQF/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

FindInstance[{INV2A<0,x11>0, x12>0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQF/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

FindInstance[{INV2A>0,x11<0, x12<0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQF/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

INV1C = f1/.EQC/.x11->0/.x21->0/.b1->1/Kb1/.w1->1/Kw1/.b2->1/Kb2/.w2->1/Kw2

CONDITIONS FOR INVASION of third population into equilibrium with two populations:

Invasion of cross path into two direct paths:

INV2E=f2/.EQE/.x12->0

INV3E = f3/.EQE/.x21->0

The following shows that these two conditions are mutually exclusive when c11=c22=1, c12>1 and c21>1.

FindInstance[{INV3E>0,INV2E>0, c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

RegionPlot[{INV2E>0/.c11->1/.c12->1.01/.c22->1/.Kb1->1/.022/.Kb2->1/.022,INV3E>0/.c11->1/.c21->1.01/.c22->1/.Kb1->1/.022/.Kb2->1/.022},{Kw1,0,1000},{Kw2,0,1000}]

Does INV2E imply EQL is positive and vice versa? The following shows that if INV2E is positive then there exists a positive equilibrium and vice versa, assuming c11=c22=1.

FindInstance[{INV2E>0,x11>0, x12>0,x22>0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQL/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]
FindInstance[{INV2E<0,x11>0, x12>0,x22>0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQL/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]
FindInstance[{INV2E>0, x12<0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQL/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

Invasion of direct into two cross paths:

INV1J=f1/.EQJ/.x11->0

INV4J=f4/.EQJ/.x22->0

The following shows that these two conditions have overlap when c11=c22=1, c12>1 and c21>1.

FindInstance[{INV1J>0,INV4J>0, c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

The following shows that at least one of the direct populations will always be able to invade the equilibrium with only cross populations.

FindInstance[{INV1J<0,INV4J<0, c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

RegionPlot[{INV1J>0/.c11->1/.c12->1.01/.c21->1.01/.c22->1/.Kb1->1/.022/.Kb2->1/.022,INV4J>0/.c11->1/.c12->1.01/.c21->1.01/.c22->1/.Kb1->1/.022/.Kb2->1/.022},{Kw1,0,1000},{Kw2,0,1000}]

Invasion of direct into community of one direct and one cross:

INV4F=f4/.EQF/.x22->0/.x21->0

INV4G=f4/.EQG/.x22->0/.x12->0

The following shows that the direct path can always invade a community of one direct and one cross path when costs are higher in cross paths:

FindInstance[{INV4G<0, c21>1,Kb1>0, Kb2>0, Kw1>0, Kw2>0}/.c11->1/.c22->1,{c21, Kb1, Kb2, Kw1, Kw2}]

Invasion of cross into community of one direct and one cross:

INV3F=f3/.EQF/.x22->0/.x21->0

The following shows that the cross population sometimes invades a community of one direct and one cross:

FindInstance[{INV3F>0, c21>1, c12>1, Kb1>0, Kb2>0, Kw1>0, Kw2>0}/.c11->1/.c22->1,{c21, c12,Kb1, Kb2, Kw1, Kw2}]

FindInstance[{INV3F<0, c21>1,c12>1, Kb1>0, Kb2>0, Kw1>0, Kw2>0}/.c11->1/.c22->1,{c21, c12,Kb1, Kb2, Kw1, Kw2}]

Does INV3F>0 imply that EQM is positive? The following shows that positive EQM implies INV3F but not the other way around.

FindInstance[{INV3F>0,x11>0, x12>0,x21>0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQM/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

FindInstance[{INV3F<0,x11>0, x12>0,x21>0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQM/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

FindInstance[{INV3F>0,x11<0, x12>0,x21>0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQM/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

Does INV1J>0 imply that EQM is positive? The following shows that positive EQM implies positive INV1J.

FindInstance[{INV1J>0,x11>0, x12>0,x21>0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQM/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

FindInstance[{INV1J<0,x11>0, x12>0,x21>0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQM/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

FindInstance[{INV1J>0,x11>0, x12<0,x21>0,c12>1, Kb1>0, Kb2>0,Kw1>0,Kw2>0, c21>1}/.EQM/.c11->1/.c22->1, {c21,Kb1, c12, Kw1, Kb2,Kw2}]

CONDITIONS FOR INVASION of fourth population into equilibrium with three populations:

Assuming c11=1, c22=1, c12>1, c21>1, the following shows a missing direct population always invades:

INV4M=Simplify[f4/.EQM/.x22->0]

INV4M/.c11->1/.c22->100/.c12->1.01/.c21->1.01

RegionPlot[{INV4M>0/.c11->1/.c22->1},{c12,1,2},{c21,1,2}]

The following shows that a missing cross pop never invades three populations at equilibrium

INV2K=Simplify[f2/.EQK/.x12->0]

FindInstance[{INV3K>0, c11>0, c22>0, c12>c11, c12>c22, c21>c11, c21>c22}, {c12, c21, c11, c22}]