Here we produce, one by one, the nine graphs that appear in Figure 3 as a 3x3 grid of images.  Each shows equilibrium values of N and P, according to our model.  The following parameters stay the same throughout each of the nine images:

NInput = IN (in the paper) = 0.2
s = 0.01
f = 0.88
k = 30.0
m = 10.0
α = 0.2
q = 2.0

The following parameters vary throughout the images:

PInput = IP = 0.004, 0.008, 0.012
r = 0.15, 0.25, 0.35

Each possible combination of PInput and r appears in one image.

Clear[NInput,PInput,sink,f,r,m,k,q,N,P,x,y];
NInput=.2;
PInput=.004;
sink=.01;
f=0.88;
r=.15;
m=10;
k=30.0;
q=2.0;
α=0.2;

cp1=ContourPlot[{NInput-sink*Nit+f*P*((E^(-α*(Nit/P -k)))/(1+(E^(-α*(Nit/P -k)))))==0,PInput-sink*P+((r*P^q)/(m^q+P^q))==0},{Nit,0,1200},{P,0,35},ContourStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red}},FrameLabel->{"N","P"},LabelStyle->Directive[Black,FontFamily->"Helvetica",FontSize->14]];

plot1=Show[cp1]

Export["figure1A pt004.eps",%]

Clear[NInput,PInput,sink,f,r,m,k,q,N,P,x,y];
NInput=.2;
PInput=.004;
sink=.01;
f=0.88;
r=.25;
m=10;
k=30.0;
q=2.0;
α=0.2;

cp2=ContourPlot[{NInput-sink*Nit+f*P*((E^(-α*(Nit/P -k)))/(1+(E^(-α*(Nit/P -k)))))==0,PInput-sink*P+((r*P^q)/(m^q+P^q))==0},{Nit,0,1200},{P,0,35},ContourStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red}},FrameLabel->{"N","P"},LabelStyle->Directive[Black,FontFamily->"Helvetica",FontSize->14]];

plot2=Show[cp2]

Export["figure1B pt004.eps",%]

Clear[NInput,PInput,sink,f,r,m,k,q,N,P,x,y];
NInput=.2;
PInput=.004;
sink=.01;
f=0.88;
r=.35;
m=10;
k=30.0;
q=2.0;
α=0.2;

cp3=ContourPlot[{NInput-sink*Nit+f*P*((E^(-α*(Nit/P -k)))/(1+(E^(-α*(Nit/P -k)))))==0,PInput-sink*P+((r*P^q)/(m^q+P^q))==0},{Nit,0,1200},{P,0,35},ContourStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red}},FrameLabel->{"N","P"},LabelStyle->Directive[Black,FontFamily->"Helvetica",FontSize->14]];

plot3=Show[cp3]

Export["figure1C pt004.eps",%]

Clear[NInput,PInput,sink,f,r,m,k,q,N,P,x,y];
NInput=.2;
PInput=.008;
sink=.01;
f=0.88;
r=.15;
m=10;
k=30.0;
q=2.0;
α=0.2;

cp4=ContourPlot[{NInput-sink*Nit+f*P*((E^(-α*(Nit/P -k)))/(1+(E^(-α*(Nit/P -k)))))==0,PInput-sink*P+((r*P^q)/(m^q+P^q))==0},{Nit,0,1200},{P,0,35},ContourStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red}},FrameLabel->{"N","P"},LabelStyle->Directive[Black,FontFamily->"Helvetica",FontSize->14]];

plot4=Show[cp4]

Export["figure2A pt008.eps",%]

Clear[NInput,PInput,sink,f,r,m,k,q,α,Nit,P,x,y];
NInput=.2;
PInput=.008;
sink=.01;
f=0.88;
r=.25;
m=10;
k=30.0;
q=2.0;
α=0.2;

cp5=ContourPlot[{NInput-sink*Nit+f*P*((E^(-α*(Nit/P -k)))/(1+(E^(-α*(Nit/P -k)))))==0,PInput-sink*P+((r*P^q)/(m^q+P^q))==0},{Nit,0,1200},{P,0,35},ContourStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red}},FrameLabel->{"N","P"},LabelStyle->Directive[Black,FontFamily->"Helvetica",FontSize->14]];
plot5=Show[cp5]

Export["figure2B pt 008.eps",%]

Clear[NInput,PInput,sink,f,r,m,k,q,N,P,x,y];
NInput=.2;
PInput=.008;
sink=.01;
f=0.88;
r=.35;
m=10;
k=30.0;
q=2.0;
α=0.2;

cp6=ContourPlot[{NInput-sink*Nit+f*P*((E^(-α*(Nit/P -k)))/(1+(E^(-α*(Nit/P -k)))))==0,PInput-sink*P+((r*P^q)/(m^q+P^q))==0},{Nit,0,1200},{P,0,35},ContourStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red}},FrameLabel->{"N","P"},LabelStyle->Directive[Black,FontFamily->"Helvetica",FontSize->14]];

plot6=Show[cp6]

Export["figure2C pt008.eps",%]

Clear[NInput,PInput,sink,f,r,m,k,q,N,P,x,y];
NInput=.2;
PInput=.012;
sink=.01;
f=0.88;
r=.15;
m=10;
k=30.0;
q=2.0;
α=0.2;

cp7=ContourPlot[{NInput-sink*Nit+f*P*((E^(-α*(Nit/P -k)))/(1+(E^(-α*(Nit/P -k)))))==0,PInput-sink*P+((r*P^q)/(m^q+P^q))==0},{Nit,0,1200},{P,0,35},ContourStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red}},FrameLabel->{"N","P"},LabelStyle->Directive[Black,FontFamily->"Helvetica",FontSize->14]];

plot7=Show[cp7]

Export["figure3A pt012.eps",%]

Clear[NInput,PInput,sink,f,r,m,k,q,N,P,x,y];
NInput=.2;
PInput=.012;
sink=.01;
f=0.88;
r=.25;
m=10;
k=30.0;
q=2.0;
α=0.2;

cp8=ContourPlot[{NInput-sink*Nit+f*P*((E^(-α*(Nit/P -k)))/(1+(E^(-α*(Nit/P -k)))))==0,PInput-sink*P+((r*P^q)/(m^q+P^q))==0},{Nit,0,1200},{P,0,35},ContourStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red}},FrameLabel->{"N","P"},LabelStyle->Directive[Black,FontFamily->"Helvetica",FontSize->14]];

plot8=Show[cp8]

Export["figure3B pt012.eps",%]

Clear[NInput,PInput,sink,f,r,m,k,q,N,P,x,y];
NInput=.2;
PInput=.012;
sink=.01;
f=0.88;
r=.35;
m=10;
k=30.0;
q=2.0;
α=0.2;

cp9=ContourPlot[{NInput-sink*Nit+f*P*((E^(-α*(Nit/P -k)))/(1+(E^(-α*(Nit/P -k)))))==0,PInput-sink*P+((r*P^q)/(m^q+P^q))==0},{Nit,0,1200},{P,0,35},ContourStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red}},FrameLabel->{"N","P"},LabelStyle->Directive[Black,FontFamily->"Helvetica",FontSize->14]];

plot9=Show[cp9]

Export["figure3C pt012.eps",%]

Below we combine the nine images into a single grid for inclusion in the article.

GraphicsGrid[{{plot1,plot2,plot3},{plot4,plot5,plot6},{plot7,plot8,plot9}}]