diff --git a/bench/pkg/catalyst.jl b/bench/pkg/catalyst.jl
index 830c58b054643de5139af3b2bc4c21a3e4709fbd..332c5bb60b91a7c2adc4b1f5a6a5296ba0a783c7 100644
--- a/bench/pkg/catalyst.jl
+++ b/bench/pkg/catalyst.jl
@@ -4,14 +4,17 @@ using MarkovProcesses
 using Catalyst
 using DiffEqJump
 
+path_latex = "./"
+
+# Bench 1: A long simulation of ER
+
 load_model("ER")
 set_param!(ER, :k1, 0.2)
 set_param!(ER, :k2, 40.0)
 ER.buffer_size = 100
 ER.estim_min_states = 8000
 
-b_pkg = @benchmark simulate(ER)
-@show minimum(b_pkg), mean(b_pkg), maximum(b_pkg)
+bench1_pkg = @benchmark simulate(ER)
 
 rs = @reaction_network begin
   c1, S + E --> SE
@@ -21,14 +24,78 @@ end c1 c2 c3
 p = (0.2,40.0,1.0)   # [c1,c2,c3]
 tspan = (0., 100.)
 u0 = [100., 100., 0., 0.]  # [S,E,SE,P]
-
 # solve JumpProblem
 dprob = DiscreteProblem(rs, u0, tspan, p)
 jprob = JumpProblem(rs, dprob, Direct())
 jsol = solve(jprob, SSAStepper())
 
-b_catalyst = @benchmark solve(jprob, SSAStepper())
-@show minimum(b_catalyst), mean(b_catalyst), maximum(b_catalyst)
+bench1_catalyst = @benchmark solve(jprob, SSAStepper())
+
+str_latex_bench1 = "
+\\begin{tabular}{|c|c|c|c|c|}
+    \\hline
+    Bench 1 & Mean time (ms) & Max. time (ms) &
+Min. time (ms) & \\begin{tabular}[c]{@{}c@{}}Mean\\\\Memory (KB)\\end{tabular} \\\\
+    \\hline
+    Package & $(round(time(mean(bench1_pkg))/1E6, digits=2))     & $(round(time(maximum(bench1_pkg))/1E6, digits=2)) &
+    $(round(time(minimum(bench1_pkg))/1E6, digits=2))  & $(round(memory(mean(bench1_pkg))/1024, digits=2)) \\\\
+    \\hline
+    Catalyst.jl  & $(round(time(mean(bench1_catalyst))/1E6, digits=2))  & $(round(time(maximum(bench1_catalyst))/1E6, digits=2)) &
+    $(round(time(minimum(bench1_catalyst))/1E6, digits=2)) & $(round(memory(mean(bench1_catalyst))/1024, digits=2)) \\\\
+    \\hline
+\\end{tabular}"
+
+# Bench 2: Creation of model + a long simulation of ER
+
+bench2_pkg = @benchmark begin
+    ER = @network_model begin
+        R1: (E+S => ES, k1*E*S)
+        R2: (ES => E+S, k2*ES)
+        R3: (ES => E+P, k3*ES)
+    end "ER"
+    set_param!(ER, :k1, 0.2)
+    set_param!(ER, :k2, 40.0)
+    ER.buffer_size = 100
+    ER.estim_min_states = 8000
+    simulate(ER)
+end
+
+bench2_catalyst = @benchmark begin
+    rs = @reaction_network begin
+        c1, S + E --> SE
+        c2, SE --> S + E
+        c3, SE --> P + E
+    end c1 c2 c3
+    p = (0.2,40.0,1.0)   # [c1,c2,c3]
+    tspan = (0., 100.)
+    u0 = [100., 100., 0., 0.]  # [S,E,SE,P]
+    # solve JumpProblem
+    dprob = DiscreteProblem(rs, u0, tspan, p)
+    jprob = JumpProblem(rs, dprob, Direct())
+    jsol = solve(jprob, SSAStepper())
+end
+
+str_latex_bench2 = "
+\\begin{tabular}{|c|c|c|c|c|}
+    \\hline
+    Bench 1 & Mean time (ms) & Max. time (ms) &
+Min. time (ms) & \\begin{tabular}[c]{@{}c@{}}Mean\\\\Memory (KB)\\end{tabular} \\\\
+    \\hline
+    Package & $(round(time(mean(bench2_pkg))/1E6, digits=2))     & $(round(time(maximum(bench2_pkg))/1E6, digits=2)) &
+    $(round(time(minimum(bench2_pkg))/1E6, digits=2))  & $(round(memory(mean(bench2_pkg))/1024, digits=2)) \\\\
+    \\hline
+    Catalyst.jl  & $(round(time(mean(bench2_catalyst))/1E6, digits=2))  & $(round(time(maximum(bench2_catalyst))/1E6, digits=2)) &
+    $(round(time(minimum(bench2_catalyst))/1E6, digits=2)) & $(round(memory(mean(bench2_catalyst))/1024, digits=2)) \\\\
+    \\hline
+\\end{tabular}"
+
+@show minimum(bench1_pkg), mean(bench1_pkg), maximum(bench1_pkg)
+@show minimum(bench1_catalyst), mean(bench1_catalyst), maximum(bench1_catalyst)
 
-#plot(jsol,lw=2,title="Gillespie: Michaelis-Menten Enzyme Kinetics")
+open(path_latex * "bench1.tex", "w+") do io
+    write(io, str_latex_bench1)
+end;
+open(path_latex * "bench2.tex", "w+") do io
+    write(io, str_latex_bench2)
+end;
 
diff --git a/models/ER.jl b/models/ER.jl
index b2d0d0c88dfabdbef74a606d91c1e0c9ce98e33f..350707d2307d30b2aaca49fe336cdadaff242442 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -52,7 +52,7 @@ isabsorbing_ER(p::Vector{Float64},xn::Vector{Int}) =
     @inbounds(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3] === 0.0)
 g_ER = [:P]
 
-ER = ContinuousTimeModel(d,k,dict_var_ER,dict_p_ER,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,isabsorbing_ER; g=g_ER,name="ER pkg")
+ER = ContinuousTimeModel(d,k,dict_var_ER,dict_p_ER,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,isabsorbing_ER; g=g_ER,name="ER SSA pkg")
 
 function create_ER(new_p::Vector{Float64})
     ER_new = deepcopy(ER)
diff --git a/models/SIR.jl b/models/SIR.jl
index 7161169bd7b103eac1a133fc9c0676a63744c1f1..e2bfadf7cb857e131931356b1d5315060ed4a53d 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -50,7 +50,7 @@ end
 isabsorbing_SIR(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
 g_SIR = [:I]
 
-SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR, name="SIR pkg")
+SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR, name="SIR SSA pkg")
 
 function create_SIR(new_p::Vector{Float64})
     SIR_new = deepcopy(SIR)