Probability of the single molecule reaction to happen
When two species $A$ and $B$ react, $A + B \to product$, the rate of reaction is $k\mathrm{[A][B]}$ where $k$ is the rate constant and $\mathrm{[A]}$ and $k\mathrm{[B]}$ the concentrations at any time during the reaction, which clearly change with time. To work out what happens in each small time interval means calculating the change in the number of reactant molecules in a given volume. Doing this means calculating two things. The first is the time at which any of possibly several reactions will occur, and the second is to decide which reaction this will be.
In the reaction $A + B \to product$, there is only one type of reaction but there are very many schemes where there are two or more reactions; for example, the photo-dissociation of halocarbons with UV light and subsequent reactions to destroy ozone.

$$\mathrm{CFCl_3} \overset{h\nu}\to \mathrm{CHCl_2  +Cl}\qquad \mathrm{Cl+O_3} \to \mathrm{ClO+O} \qquad \mathrm{ClO+O_2} \to \mathrm{Cl+O_2} $$ 

Given such a reaction scheme, the first step is to decide what reaction will happen next given the current state of the reaction at time $t$. At the level of considering individual molecules, once the reaction has started, it does not follow sequentially as shown in the scheme. This only happens when a whole ensemble of molecules, typically as large as Avogadro's number is examined, and then normal kinetic equations apply because of the averaging over many events that occurs. In any given time interval when considering molecules one at a time, it is not possible to predict exactly which reaction will happen next. It is quite possible, for instance, that several $\mathrm{Cl+O_3}$ reactions may occur before either of the other two. The best that can be done is to guess the outcome, and repeat the calculation many times and this approach is now described.

Suppose that the reactions are labelled with a number $R$ that represents reaction 1 or 2 or 3 etc. The chance of a reaction of type $R$ occurring will be defined as

$$prob(T + dT, R) \equiv a_RdT$$

Specifically, this is the probability that a reaction of type $R$ occurs in time interval $T$ to $T + dT$ in the reaction volume, given the current state of all reactions at an earlier time $t$. The values of $a_R$ can be easily worked out, and this is shown later. What is now needed is the probability $P_0(T)$ that no reaction of any type occurs up to time $T$ past the present time $t$, multiplied by the probability that only a reaction of type $R$ occurs in the small time interval $T \to T + dT$, which is defined as $a_RdT$. The probability of reaction type $R$ happening is therefore 

$$ P(R,T)= P_0(T)a_RdT$$

and shown schematically in the figure.



Calculation of reaction probability.



To calculate the probability of no reaction occurring, consider any small interval $d\tau$ after some arbitrary time $\tau$, making an interval $\tau \to \tau + d\tau$. The probability of no reaction is then

$$P_0(\tau+d\tau) =P_0(\tau)\left(1-\sum_{R=1}^n a_Rd\tau  \right)$$ 

where the summation is over all $n$ reaction types $R$. The summation is made because this gives the chance that any one of the $R$ reactions will occur, and therefore one minus this number is the chance that no reaction occurs. Rearranged, this equation gives

$$ \frac{P_0(\tau+d\tau)-P_0(\tau)}{d\tau}=-\sum_{R=1}^n a_R$$

and in the limit of a small time interval $d\tau$, it becomes the true derivative $\displaystyle \frac{dP_0}{d\tau}=-\sum_{R=1}^na_R$.

For clarity only, a notational change is now made so that

$$a_0=\sum_{R=1}^n a_R \tag{8}$$

and therefore $\displaystyle \frac{dP_0}{d\tau}=-a_0$.The chance of no reaction up to time $T$ is found by integrating, and is $\displaystyle P_0(T)=e^{-a_0T}$.

The final result for the probability of reaction R occurring at time $T$ is the chance of no reaction up to $T $ multiplied by the chance of reacting at $T$, which is $a_R$, giving

$$P(T,R) = a_Re^{-a_0T} \tag{9}$$

an expression that shows that these times are Poisson distributed. Equation 9 is a joint probability density function, as it is the chance that reaction occurs in the time interval and that it will be reaction $R$.

Two terms now remain undefined. The first is the individual $a_R$ and their sum $a_0$, and the second determines which reaction happens next. Calculating $a_R$ is very easy, as each $a_R$ is the product of the number of molecules present in a reaction of type $R$ multiplied with the rate constant per unit volume. 

If there is just a single reaction $\mathrm{A} \overset{k_1}\to \mathrm{B}$ then $a_1=k_1N_A$ where $N_A$ is the number of $A$ molecules present at time $T$. If the reaction is $\mathrm{A+B} \overset{k_2} \to \mathrm{C}$ then $a_2=k_2N_aN_B$ and if both reactions are present in the overall reaction scheme which would then be

$$\mathrm{A} \overset{k_1}\to \mathrm{B} \qquad \mathrm{A+B} \overset{k_2} \to \mathrm{C} \qquad a_0=\sum a_R=a_1+a_2 \tag{10}$$

In the special case where a reaction combines two identical molecules $\mathrm{A+A} \overset{k_2} \to $ the number of indistinguishable combinations, i.e the number of distinct A-A pairs of the two species must be calculated and this makes $\displaystyle a = k_2\frac{N_A(N_A - 1)}{2!}$. Normally, in dealing with chemical kinetics, $N_A$ is so vast that $\displaystyle a = k_2\frac{N_A^2}{2}$, but in the Monte Carlo simulations, $N_A$ is small because individual molecules are dealt with, and the exact formula must be used. Thus, macroscopic, (kinetic), chemical equations are only valid when deviations from the Poisson distribution are negligible.



An aside

The term $\displaystyle a = k_2\frac{N_A(N_A - 1)}{2!}$ implies that these numbers are averages in the reaction volume at time $t$ and should be written as the average of the product, i.e. as $\displaystyle \frac{\langle N_a(t)(N_A(t)-1)\rangle }{2!}$. The question then arises as how to deal with this. Two approximations can be made; the first is to assume that $N_A$ is large compared to 1 and so $\langle N_A(t)(N_A(t)-1)\rangle \rightarrow \langle N_A^2(t)\rangle $, which seems to be a reasonably sensible approximation, the next step is more problematic and is $\langle N_A^2(t)\rangle \sim \langle N_A(t)\rangle^2$. This is a drastic step and assumes that the variance in the number of molecules is zero, i.e. that $N_A(t)=\langle N_A(t)\rangle$ which means that $N_A$ is a deterministic process, so that the random nature of the numbers of $N_A$ entering into the reaction volume is always the same, which is not what was assumed at the outset. (To see that the variance is zero use the definition $\sigma^2_{N_A} = \langle N_A^2 \rangle -\langle N_A \rangle^2$ where $\sigma$ is the standard deviation) 

In any normal kinetic experiment the number of molecules is vast, $\langle N_A(t)\rangle \sim 10^{18}$ would not be untypical, so that the fluctuations in the  number of molecules is tiny and the approximation made above is valid. 

$$ \frac{\sigma_{N_A}}{\langle N_A \rangle} \sim \frac{\sqrt{\langle N_A \rangle}}{\langle N_A \rangle} \sim \frac{1}{\sqrt{\langle N_A \rangle}}\sim 10^{-9}$$

Even if there are as few as 1000 molecules the relative fluctuation  is not too drastic but in experiments using small numbers of molecules, as in the case of many 'single molecule' experiments then clearly the rate equation approach needs to be modified. 



The final part of the calculation determines which particular reaction step is going to occur and at what time this is done by choosing pairs of numbers that have the distribution $P(T, R)$ equation 9. Gillespie (1997, 2007) argues that $P(T, R)$ is equivalent to two distributions, $P_1$ and $P_2$;

$$P(T,R)= P_1\cdot P_2 = \left(a_0e^{-a_0 T} \right)\left( \frac{a_R}{a_0}\right)$$

and which can be determined by choosing two uniformly distributed random numbers, $r$, each in the range 0 to 1. The equation

$$T=-\frac{1}{a_0}\ln(r) \tag{11}$$

generates a random number $T$ according to the probability density function $\displaystyle P_1=a_0e^{-a_0 T}$. 

Equation 11 produces times that are Poisson distributed, and the time $T$ is added to the current time of the reaction. The second distribution, $P_2$, is calculated by choosing a uniform random number according to the probability density function $a_R/a_0$. For example, the chance that reaction 1, (equation 10), occurs is $a_1/(a_1 + a_2)$ and for reaction 2 this chance is $1 - a_1/(a_1 + a_2) = a_2/(a_1 + a_2)$. Which reaction occurs is determined by choosing a random number between 0 and 1 and seeing which of the two regions it falls into. If it falls between 0 and $a_1/(a_1 + a_2)$, reaction 1 occurs; otherwise, reaction 2 occurs. If reaction 1 occurs, $\mathrm{A} \to \mathrm{B}$, then 1 is subtracted from the original number of $A$ molecules in each step of the reaction and 1 is added to the number of $B$ molecules.

If reaction 2 happens, $\mathrm{A+B} \overset{k_2} \to \mathrm{C}$ then 1 molecule is removed from each of the numbers for $A$ and $B$. How many molecules are added and subtracted, of course, depends on the stoicheiometry of each particular reaction.

In the case of a termolecular reaction $\mathrm{ A + B + C} \rightarrow $ when the probability of three molecules simultaneously colliding is calculated the probability the chance that this occurs in  an infinitesimal time interval $d t$ is proportional to $dt^2$ so becomes vanishingly small and thus does not occur with any physical reality. It is often convenient, however, to consider a termolecular scheme as a simplified approximation to a scheme such as $N_A+N_B \leftrightharpoons X^*; \; X^* + N_C \rightarrow$ 

The outline plan for the Gillespie algorithm is 

(1) Define the initial amounts of each species. 


  Define the arrays of data points to hold results.   Set the maximum time. Set the time to zero; t = 0.


(2) Start a loop around the calculation until the time is up or no molecules are left. 


  Calculate $a_1,\;a_2,\;a_3$ etc. and $a_0 =a_1 +a_2 +a_3 +\cdots$ Calculate time $\displaystyle T = -\frac{1}{a_0}\ln(r)$
  Increment time $t = t + T$
  Calculate which reaction occurs: if a random number falls between 0 and $a_1/a_0$, reaction 1 occurs, and so forth for other reactions.
  Add or subtract numbers of molecules from the totals, depending on which reaction has occurred and its stoicheiometry.
  Store the results.


(3) Continue the loop until finished.



Refs

'Markov Processes. An Introduction for Physical Scientists' D. T. Gillespie, publ Academic Press 1992.

D. Gillespie  J. Phys. Chem. 81, p2340, 1977 eqn 29

The python code for a Lotka-Volterra predator-prey reaction is given below. I'm assuming that matplotlib and numpy, as np.blabla() etc are used.  

def lotka():                # Lotka reaction Gillespie method      
                            # x0 + Y = 2Y  (k1) ,  Y +  Z = 2Z (k2) , Z = Q(k3)

    #k1 =10.0   k2 = 0.01   k3 = 10.0
    #y0 =2000.0  z0 =2000.0
    ycount = np.zeros(bins,dtype=int)
    zcount = np.zeros(bins,dtype=int)
    dtime =  np.zeros(bins,dtype=float)
    for i in range(bins):
        dtime[i] = maxt*(i+1/2)/bins
        pass
    ny = y0
    nz = z0
    ycount[0] = y0
    zcount[0] = z0
    t = 0.0
    indx = 0
    while indx &lt; bins :
        r1 = k1*ny
        r2 = k2*nz*ny
        r3 = k3*nz
        a0 = r1 + r2 + r3+1e-20
        t0 = -np.log( np.random.rand() )/a0
        t = t + t0
        indx = int(np.trunc( t*bins/maxt ))
        if indx &lt; bins :
            r = np.random.rand()
            if  r &gt; ( r1 + r2 )/a0  :
                nz = nz - 1
                zcount[indx] = nz          # react 3
            elif r &gt; r1/a0 :
                nz = nz + 1
                ny = ny - 1
                zcount[indx] = nz
                ycount[indx] = ny          # react  2
            else:
                ny = ny + 1
                ycount[indx] = ny          # react 1
                pass
            pass
        pass         # events

    return zcount,ycount,dtime      # end lotka

    bins = 1000
    maxt = 200.0
    k1 = 0.5
    k2 = 0.001
    k3 = 0.5
    y0 = 400.0
    z0 = 100.0

    zcount,ycount,dtime = lotka()

    fig1= plt.figure(figsize=(10, 6))
    plt.rcParams.update({'font.size' : 16})  # set font size for plots

    plt.plot(dtime,zcount,color='red')
    plt.plot(dtime,ycount,color='blue')#'ro',markersize=2)

    plt.xlabel('time')
    plt.ylabel('population')
    plt.xlim([0,maxt])

    plt.tight_layout()

    plt.show()