Code optimization

Question

I designed this Gillespie's algorithm by the First Reaction method, and I would like this runs faster. I mean, to reduce the absolute timing of the Do bucle. So, I was thinking what kind of Mathematica's functions will make it faster. Could be with a Compile function?

pF[abundancias_] := Through[reactionsPF @@ abundancias];

listData = {}; 
times = {}; 
t = 0; 
AppendTo[times, t];
abundancia = {1000, 100, 1, 0, 100}(*S,E,I,R,V, 5 species*);
AppendTo[listData, abundancia];

Do[tis = 
  ParallelMap[
   If[# > 0, 
     RandomVariate[ExponentialDistribution[#]], \[Infinity]] &, 
   pF[abundancia](*Funciones de propensión para cada reacción*)]; 
 reaccionYt = {FirstPosition[#, Min[#]], Min[#]} &[\[Tau]is]; 
 abundancia = Through[{reaccionYt[[1]] /. changesByRXN}[[1, 1]] @@ abundancia][[
   1]];
 AppendTo[listData, abundancia]; 
 t = t + reaccionYt[[2]]; 
 AppendTo[times, t]; 
 Print["t: ", t], 40000]

Where reactionPF and changesbyRXN are list and associations array of pure functions:

reactionsPF = 
 Function /@ {.., \[Mu] #1,..};

changesByRXN = 
  Function /@ <|1 -> {#1 + 1, #2, #3, #4, #5}, 
    2 -> {#1 - 1, #2, #3, #4, #5}, ....|>;