Gillespie's event driven algorithm for SIR model

Question

The SIR model is implemented as a event driven model.

The SIR model is:

$\frac{dX}{dt} =\mu N-\beta X{Y\over N}-\mu X$
$\frac{dY}{dt}= \beta X{Y\over N} - \gamma Y-\mu Y$
$\frac{dZ}{dt}= \gamma Y-\mu Z$

The events are

• Births occur at rate $\mu N$. Result: $X\to X+1.$
• Transmission occurs at rate $\beta\frac{XY}{N}$. Result: $Y\to Y + 1$ and $X \to X-1$.
• Recovery occurs at rate $\gamma Y$. Result: $Z\to Z+1$ and $Y\to Y - 1.$
• Deaths of $X, Y, $ or $Z$ (three independent events) occur at rate $\mu X, \mu Y, $ and $\mu Z.$ Result: $X\to X-1,Y\to Y-1$ or $Z\to Z-1$.

The algorithm is:

1. Label all possible events $E_1,\dots, E_n$.
2. For each event determine the rate at which it occurs, $R_1, \dots, R_n$.
3. The rate at which any event occurs is $R_\text{total} = \sum_{m=1}^n R_m$.
4. The time until the next event is $\delta t = \frac{-1}{R_\text{total}}\log(RAND_1)$.
5. Generate a new random number, $RAND_2$. Set $P = RAND_2\times R_\text{total}$.
6. Event $p$ occurs if \begin{align} \sum_{m=1}^{p-1}R_m < P \leq \sum_{m=1}^{p}R_m.\end{align}
7. The time is now updated $t\to t+\delta t$ and event $p$ is performed.
8. Return to step 2.

This is the MATLAB code that I wrote. In this I don't understand how the next event is determined (Step 6).

function GillespieSIR()
if nargin==0
   beta=1;
   gamma=1/10;
   mu=5e-4;
   N=5000;
   X0=floor(gamma*N/beta);
   Y0=ceil(mu*N/gamma);
   MaxTime=2*365;
end

t=0;
X=X0;
Y=Y0;
Z=N-X-Y;
susceptible=X0;
infected=Y0;
recovered=N-X0-Y0;
time=0;
while(t<MaxTime )
    rate_transmission=beta*X*Y/N;
    rate_births=mu*N;
    rate_recovery=gamma*Y;
    death_rate_S=mu*X;
    death_rate_Y=mu*Y;
    death_rate_Z=mu*Z;
    total_rate=rate_transmission+rate_births+rate_recovery+death_rate_S+death_rate_Y+death_rate_Z;
    time_random=rand(1);
    time_to_next_event=-1*log(time_random)/total_rate;
    cumulative_rate=rate_births;
    event_random=rand(1);
    if event_random < (cumulative_rate/total_rate)
        X=X+1;
        cumulative_rate=cumulative_rate+rate_transmission;
    elseif event_random < (cumulative_rate/total_rate)
        X=X-1;
        Y=Y+1;
       cumulative_rate=cumulative_rate+rate_recovery;
    elseif event_random < (cumulative_rate/total_rate) 
        Y=Y-1;
        Z=Z+1;
        cumulative_rate=cumulative_rate+death_rate_S;
    elseif event_random < (cumulative_rate/total_rate)
        X=X-1;
        cumulative_rate=cumulative_rate+death_rate_Y;
    elseif event_random < (cumulative_rate/total_rate)
        Y=Y-1;
        cumulative_rate=cumulative_rate+cumulative_rate_Z;
    else
        Z=Z-1;
    end
    t=t+time_to_next_event;
    susceptible=[susceptible;X];
    infected=[infected;Y];
    recovered=[recovered;Z];
    time=[time;t];
end

figure
subplot(3,1,1)
plot(time/365,susceptible)
xlabel('time')
ylabel('susceptible')
subplot(3,1,2)
plot(time/365,infected,'r')
xlabel('time')
ylabel('infected')
subplot(3,1,3)
plot(time/365,recovered,'k')
xlabel('time')
ylabel('recovered')  

How can I select the next event. Is the way I am selecting the next event wrong because I am not getting the plot as it should look.
To find the next event the code in the book has used m=min(find(cumsum(Rate)>=R1*sum(Rate))); where R1 is a random number.

When $p=1$, how is the first event selected?
Please let me know what I am doing wrong in the MATLAB code.