Clear["Global`*"]

Manipulate[
e1=0.5;
e2=0.5;
r1=0.2;
r2=R2;
k=500;
d=0.5;
theta1=-0.2;
theta2=0.2;
alpha=0.02;
tau=0.1;
sigma=0.2;
h=1;
tmax=1000;
s=NDSolve[{
p'[t]==e1*(((alpha*tau)/Sqrt[sigma^2+tau^2])*Exp[-((xbar[t]-theta1)^2)/(2*(sigma^2+tau^2))])*n1[t]*p[t]+e2*(((alpha*tau)/Sqrt[sigma^2+tau^2])*Exp[-((xbar[t]-theta2)^2)/(2*(sigma^2+tau^2))])*n2[t]*p[t]-d*p[t],
n1'[t]==r1*n1[t]*(1-(n1[t]/k))-(((alpha*tau)/Sqrt[sigma^2+tau^2])*Exp[-((xbar[t]-theta1)^2)/(2*(sigma^2+tau^2))])*n1[t]*p[t],
n2'[t]==r2*n2[t]*(1-(n2[t]/k))-(((alpha*tau)/Sqrt[sigma^2+tau^2])*Exp[-((xbar[t]-theta2)^2)/(2*(sigma^2+tau^2))])*n2[t]*p[t],
xbar'[t] == (Sqrt[h*(sigma^2)])^2*(((e1*n1[t]*tau*alpha*(theta1-xbar[t])/((tau^2+sigma^2)^(3/2)))*Exp[-(xbar[t]-theta1)^2/(2*(tau^2+sigma^2))])+((e2*n2[t]*tau*alpha*(theta2-xbar[t])/((tau^2+sigma^2)^(3/2)))*Exp[-(xbar[t]-theta2)^2/(2*(tau^2+sigma^2))])),
p[0]==5,n1[0]==1,n2[0]==120,xbar[0]==0.2},
{p,n1,n2,xbar},{t,tmax}];
GraphicsGrid[{{Plot[Evaluate[{p[t],n1[t],n2[t]}/.s],{t,0,tmax}],
Plot[Evaluate[{xbar[t]}/.s],{t,0,tmax},PlotRange->{-0.3,0.3}]}}],
{R2,0.001,1}]