You try to make a stochastic model corresponding to the deterministic model above. The Master equation of that model is
$\begin{eqnarray}\frac{d}{dt}\mathbb{P}(I,P) & = & \mathbb{T}(I,P|I-1,P)\cdot \mathbb{P}(I-1,P) - \mathbb{T}(I+1,P|I,P)\cdot \mathbb{P}(I,P) \\ & & + \mathbb{T}(I,P|I+1,P)\cdot \mathbb{P}(I+1,P) - \mathbb{T}(I-1,P|I,P)\cdot \mathbb{P}(I,P) \\ & & + \mathbb{T}(I,P|I,P-1)\cdot \mathbb{P}(I,P-1) - \mathbb{T}(I,P+1|I,P)\cdot \mathbb{P}(I,P) \\ & & + \mathbb{T}(I,P|I,P+1)\cdot \mathbb{P}(I,P+1) - \mathbb{T}(I,P-1|I,P)\cdot \mathbb{P}(I,P) \; . \end{eqnarray}$
Where the transition probabilities are
$\begin{eqnarray} \mathbb{T}(I+1,P|I,P) & = & \alpha_I+P\lambda_I\\ \mathbb{T}(I-1,P|I,P) & = & I\gamma_I-PI\delta_I\\ \mathbb{T}(I,P+1|I,P) & = & \alpha_P+P\theta_P\\ \mathbb{T}(I,P-1|I,P) & = & PI\delta_P \end{eqnarray}$
It takes a bit of work, but computing the expectation values $\mathbb{E}(I)=\sum_{I,P} I\cdot \mathbb{P}(I,P)$ and $\mathbb{E}(P)=\sum_{I,P} P \cdot\mathbb{P}(I,P)$, can be done by multiplying the master equation by either $I$ or $P$ and then summing over all values that $I$ and $P$ can get. This will give you the following equations:
$\begin{cases} \frac{\mathrm{d}\mathbb{E}(I)}{\mathrm{d}t}=& \alpha_I+\lambda_I\mathbb{E}(P)- \gamma_I\mathbb{E}(I) -\delta_I\mathbb{E}(PI) \\ \frac{\mathrm{d}\mathbb{E}(P)}{\mathrm{d}t}=&\alpha_P+\theta_P\mathbb{E}(P)-\delta_P\mathbb{E}(PI) \end{cases}$
Note however that you get the expectation of the product $PI$ and not the product of expectations.
EDIT:
Here's my code in R for a run of the stochastic model
ai<-0; # parameters gi<-0.1; di<-0.05; li<-0.05; tp<-0.15; dp<-0.01; ap<-0.5; T<-10; plusI<-c(1,-1,0,0); plusP<-c(0,0,1,-1); Time<-c(0); # initial conditions I<-c(20); P<-c(5); for (k in 1:800) { freq<-ai+li*P[k]+gi*I[k]+di*P[k]*I[k]+ap+tp*P[k]+dp*P[k]*I[k]; t<-rexp(1,freq); Time<-c(Time,Time[k]+t); sel<-sample(c(1,2,3,4),1,prob=c(ai+li*P[k],gi*I[k]+di*P[k]*I[k],ap*(Time[k]
To explain the most essential part, I draw a time from an exponential distribution with parameter the sum of all the rates of the processes. Then, I decide which process takes place by sampling with a probability vector equal to the 4 rates of the processes. This way, I manage to simulate the stochastic process. Here's a run of it:

I haven't tried to explore the state space in detail yet, but all choices of initial conditions I have tried lead to an exponential explosion of the rabbit population and a collapse of the fox population.