Stochastic predator-prey

Question

My system is a simple $P$ vs $I$ foxes- vs rabbits model given by: $$ \begin{cases} \frac{\mathrm{d}I}{\mathrm{d}t}=& \alpha_I+\lambda_IP- \gamma_II -\delta_IPI;\\ \frac{\mathrm{d}P}{\mathrm{d}t}=&\alpha_P+\theta_PP-\delta_PPI \end{cases} $$

with a parameter set: $$ \begin{cases} \theta_P &=0.15 \\ \delta_P&=0.01 \end{cases}\quad\text{ and }\quad \begin{cases} \alpha_I&=0.4 \\ \gamma_I&=0.1 \\ \delta_I&=0.05 \\ \lambda_I&=0.05 \end{cases} $$ but a condition on the initial introduction of $D$ rabbits ($P$) over a specific timeframe $T$.

$\alpha_P = \begin{cases} &\dfrac{D}{T} &\mbox{if } t\leq T \\ \\ &0 & \mbox{otherwise.} \end{cases}$

METHOD I'm using is: over any small timestep $\delta t<<1$, one of the following can occur:

I=I+1, with probability $\alpha_I+P\lambda_I$

I=I-1, with probability $I\gamma_I-PI\delta_I$

P=P+1, with probability $ P\theta_P$

P=P-1, with probability $ PI\delta_P$

QUESTION: But is this true when $t<T$?

EDIT 18/04/2013:

Consider that P is actually Pathogens, and I is actually immunity cells in the human body. A literature search finds: (where Innoculation time is $t<T$)

Pujol 2009 -The Effect of Ongoing Exposure Dynamics in Dose Response Relationships (free access)

Answer 1

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]<T)+tp*P[k],dp*P[k]*I[k]));
I<-c(I,I[k]+plusI[sel]);
P<-c(P,P[k]+plusP[sel]);
}

layout(matrix(c(1,2),1,2));

plot(Time,I,type="l",col="blue");
plot(Time,P,type="l",col="red");

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.