Stochastic predator-prey 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)
 A: 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.
A: I'll try to clear up the confusion you have about the meaning of $P$ and $E(P)$. 
First, the meaning of $P$. $P$ is a stochastic variable. This means that $P$ is really a function from the space of all possible realizations of the stochastic process to a number (or, since $P$ is time-dependent here, a real function of time). 
How does this connect with the picture? Each curve corresponds to one realization of the process. Let's call those $\omega_i$ with $i$ and index to label the realizations. So what is plotted is really $P(\omega_i)$ for several different $\omega_i$.
They then talk about the percentage of infections. Apparently, infection is defined as a situation where the number of pathogens increases indefinitely. So, all you have to do to determine that percentage is to count how many $\omega_i$ will have $P(\omega_i)$ diverging. Since it's impossible to do for all possible realizations, they only use a sample of them to estimate the percentage. If the sample is large enough, the estimate will be good.
What about $E(P)$? $E(P)$ is not a function of the realizations $\omega_i$ like $P$ is. In fact, $E(P)$ is the mean of $P$ over all possible realizations. One could formally notate this as
$$E(P)=\sum_{\omega_i}P(\omega_i)\text{Prob}(\omega_i) \; .$$
But since by simulating we don't have access to all possible realizations and the probability distribution of those, we again just make an estimate by sampling. All you have to do is take the arithmetic mean over all sampled realizations.
But now you face a problem. Since some trajectories diverge, they will dominate the mean. A smarter thing to do would be to calculate maybe two separate means: one for diverging trajectories and one for converging trajectories. 
I hope this helps clearing up the confusion somewhat. If you're still stuck, I think you should look more seriously into stochastic processes. May I advise you the following book:
N. G. van Kampen: Stochastic processes in physics and chemistry, North Holland 1981, 3rd edn., 2007, ISBN 0-444-89349-0
I want to stress one point again though: $E(P)$ is not necessary to compute the percentage of infections at all!
EDIT: My slightly altered code in R:
    ai<-0.4; # parameters
gi<-0.01;
di<-0.005;
li<-0.05;
tp<-0.15;
dp<-0.01;
T<-1;
ap<-25/T;

dt<-0.01; # timestep

I<-c(2); # initial conditions
P<-c(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:10000) { 
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(log(Time)/log(10),log(I)/log(10),type="l",col="blue");
plot(log(Time)/log(10),log(P)/log(10),type="l",col="red");

