P[x_,y_]:=a x[t]y[t]

Q[x_,y_]:=((b x[t]+d y[t]-b x0) (x[t]-1)+(c y[t]^2))

x0:=0.1

b:=0.1

c:=0.01

d:=-0.01

a:=0.05

θ:=-0.00001

Manipulate[Module[{sol=NDSolve[{x'[t]==a x[t]y[t] Cos[θ]-((b x[t]+d y[t]-b x0) (x[t]-1)+(c y[t]^2)) Sin[θ],y'[t]==a x[t]y[t]  Sin[θ]+((b x[t]+d y[t]-b x0) (x[t]-1)+(c y[t]^2)) Cos[θ],x[0]==p[[1]],y[0]==p[[2]]},{x,y},{t,-5000,0}]},ParametricPlot[Evaluate[{x[t],y[t]}/.sol],{t,-5000,0},PlotRange->{{-1.2,1.2},{-1.2,1.2}}]],{{p,{1.01,1.01}},Locator},{{θ,-0.00001},-0.1,0.1},{{a,0.05},0,0.1},{{b,0.1},0,0.2},{{d,-0.01},-0.05,0}]