diff --git a/bench/pygmalion/access_val_er.jl b/bench/pygmalion/access_val_er.jl
new file mode 100644
index 0000000000000000000000000000000000000000..05c444c5d43a69ceec90a45d03b3898824be0c10
--- /dev/null
+++ b/bench/pygmalion/access_val_er.jl
@@ -0,0 +1,41 @@
+
+using BenchmarkTools
+using pygmalion
+using MarkovProcesses
+
+println("Pygmalion:")
+
+str_m = "enzymatic_reaction"
+str_d = "abc_er"
+pygmalion.load_model(str_m)
+str_oml = "P,R,time"
+ll_om = split(str_oml, ",")
+x0 = State(95.0, 95.0, 0.0, 0.0, 0.0, 0.0)
+p_true = Parameters(1.0, 1.0, 1.0)
+u = Control(10.0)
+tml = 1:400
+g_all = create_observation_function([ObserverModel(str_oml, tml)]) 
+so = pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true)
+function random_trajectory_value_pyg(so::SystemObservation)
+    n_states = get_number_of_observations(so, "P")
+    return to_vec(so, "P", rand(1:n_states))
+end
+# Bench
+@timev random_trajectory_value_pyg(so)
+b1_pyg = @benchmark random_trajectory_value_pyg($so)
+@show minimum(b1_pyg), mean(b1_pyg), maximum(b1_pyg)
+
+println("MarkovProcesses:")
+
+MarkovProcesses.load_model("ER")
+ER.time_bound = 10.0
+σ = MarkovProcesses.simulate(ER)
+function random_trajectory_value(σ::AbstractTrajectory)
+    n_states = get_states_number(σ)
+    return σ["P"][rand(1:n_states)]
+end
+# Bench
+@timev random_trajectory_value(σ) 
+b1 = @benchmark random_trajectory_value($σ)
+@show minimum(b1), mean(b1), maximum(b1)
+
diff --git a/bench/pygmalion/access_var_values.jl b/bench/pygmalion/access_var_values.jl
new file mode 100644
index 0000000000000000000000000000000000000000..a108da24f8418c53c4a3cfb78db3b3c58b532641
--- /dev/null
+++ b/bench/pygmalion/access_var_values.jl
@@ -0,0 +1,62 @@
+
+using BenchmarkTools
+using pygmalion
+using MarkovProcesses
+
+println("Pygmalion:")
+
+str_m = "enzymatic_reaction"
+str_d = "abc_er"
+pygmalion.load_model(str_m)
+str_oml = "P,R,time"
+ll_om = split(str_oml, ",")
+x0 = State(95.0, 95.0, 0.0, 0.0, 0.0, 0.0)
+p_true = Parameters(1.0, 1.0, 1.0)
+u = Control(10.0)
+tml = 1:400
+g_all = create_observation_function([ObserverModel(str_oml, tml)]) 
+so = pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true)
+function read_trajectory_v1(so::SystemObservation)
+    vals = to_vec(so, "P")
+    n_states = length(vals)
+    res = 0.0
+    for i = 1:n_states
+        res += vals[i]
+    end
+    return res
+end
+function read_trajectory_v2(so::SystemObservation)
+    idx_P = om_findfirst("P", so.oml)
+    n_states = length(so.otll[idx_P]) 
+    res = 0.0
+    for i = 1:n_states
+        res += so.otll[idx_P][i][2][1]
+    end
+    return res
+end
+# Bench
+@timev read_trajectory_v1(so)
+b1_pyg = @benchmark read_trajectory_v1($so) 
+@show minimum(b1_pyg), mean(b1_pyg), maximum(b1_pyg)
+@timev read_trajectory_v2(so)
+b2_pyg = @benchmark read_trajectory_v2($so) 
+@show minimum(b2_pyg), mean(b2_pyg), maximum(b2_pyg)
+
+println("MarkovProcesses:")
+
+MarkovProcesses.load_model("ER")
+ER.time_bound = 10.0
+σ = MarkovProcesses.simulate(ER)
+function read_trajectory(σ::AbstractTrajectory)
+    n_states = get_states_number(σ)
+    res = 0
+    for i = 1:n_states
+        res += (σ["P"])[i]
+    end
+    return res
+end
+# Bench
+@timev read_trajectory(σ) 
+b1 = @benchmark read_trajectory($σ)
+@show minimum(b1), mean(b1), maximum(b1)
+
diff --git a/bench/pygmalion/long_sim_er.jl b/bench/pygmalion/long_sim_er.jl
new file mode 100644
index 0000000000000000000000000000000000000000..6a066ba129051f41bfe5cbf4b6303d5095f87855
--- /dev/null
+++ b/bench/pygmalion/long_sim_er.jl
@@ -0,0 +1,35 @@
+
+using BenchmarkTools
+
+println("Pygmalion:")
+
+using pygmalion
+str_m = "enzymatic_reaction"
+str_d = "abc_er"
+pygmalion.load_model(str_m)
+str_oml = "P,R,time"
+ll_om = split(str_oml, ",")
+x0 = State(95.0, 95.0, 0.0, 0.0, 0.0, 0.0)
+p_true = Parameters(1.0, 1.0, 1.0)
+u = Control(10.0)
+tml = 1:400
+g_all = create_observation_function([ObserverModel(str_oml, tml)]) 
+b1_pyg = @benchmark pygmalion.simulate($f, $g_all, $x0, $u, $p_true; on = nothing, full_timeline = true)
+b2_pyg = @benchmark pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true)
+
+@timev pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true)
+@show minimum(b1_pyg), mean(b1_pyg), maximum(b1_pyg)
+@show minimum(b2_pyg), mean(b2_pyg), maximum(b2_pyg)
+
+println("MarkovProcesses:")
+
+using MarkovProcesses
+MarkovProcesses.load_model("ER")
+ER.time_bound = 10.0
+b1 = @benchmark MarkovProcesses.simulate($ER)
+b2 = @benchmark MarkovProcesses.simulate(ER)
+
+@timev MarkovProcesses.simulate(ER)
+@show minimum(b1), mean(b1), maximum(b1)
+@show minimum(b2), mean(b2), maximum(b2)
+
diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index eec79d03650ad8209f3d18f9f88901bb9692d71c..04ac4c7ad6f5bc108b4d180f4de93d0a5598df57 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -8,7 +8,7 @@ export load_model
 include("model.jl")
 
 export Observations, AbstractTrajectory
-export +,-,δ,get_obs_variables
+export +,-,δ,get_obs_variables,get_states_number
 include("observations.jl")
 
 end
diff --git a/core/model.jl b/core/model.jl
index ed108270f69e28dcc1db055e5c97fda659acda8b..5e87907d3f768fa7d084aa40453da2ab663f781b 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -38,7 +38,7 @@ end
 function simulate(m::ContinuousTimeModel)
     full_values = zeros(m.d, 0)
     times = zeros(0)
-    transitions = Vector{String}(undef, 0)
+    transitions = Vector{Union{String,Nothing}}(undef, 0)
     xn = m.x0
     tn = m.t0 
     n = 0
@@ -51,6 +51,14 @@ function simulate(m::ContinuousTimeModel)
         n += 1
     end
     values = full_values[m._g_idx,:] 
+    if is_bounded(m)
+        if times[end] > m.time_bound
+            values[:, end] = values[:,end-1]
+            times[end] = m.time_bound
+            transitions[end] = nothing
+        end
+    end
+    
     return Trajectory(m, values, times, transitions)
 end
 
@@ -62,6 +70,7 @@ function simulate(m::ContinuousTimeModel, n::Int)
     return obs
 end
 
+is_bounded(m::Model) = m.time_bound < Inf
 function check_consistency(m::Model) end
 function simulate(m::Model, n::Int; bound::Real = Inf)::AbstractObservations end
 function set_param!(m::Model, p::AbstractVector{Real})::Nothing end
diff --git a/core/observations.jl b/core/observations.jl
index 86fca14a7de6e8f6382d3d263f5cc6710a899170..7a35b81907c40ed965f52b8ae1310353feb92d75 100644
--- a/core/observations.jl
+++ b/core/observations.jl
@@ -6,7 +6,7 @@ struct Trajectory <: AbstractTrajectory
     m::ContinuousTimeModel
     values::AbstractMatrix{Real}
     times::AbstractVector{Real}
-    transitions::AbstractVector{Union{String,Missing}}
+    transitions::AbstractVector{Union{String,Nothing}}
 end
 
 function +(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
@@ -15,7 +15,10 @@ function δ(σ1::AbstractTrajectory,t::Real) end
 function get_obs_variables(σ::AbstractTrajectory) end
 
 get_values(σ::AbstractTrajectory, var::String) = 
-    σ.values[(σ.m)._map_obs_var_idx[var],:] 
+σ.values[(σ.m)._map_obs_var_idx[var],:] 
+
+get_states_number(σ::AbstractTrajectory) =
+length(σ.times)
 
 function getindex(σ::AbstractTrajectory, idx::String)
     if idx  == "times"
diff --git a/models/ER.jl b/models/ER.jl
new file mode 100644
index 0000000000000000000000000000000000000000..876b129a7510a3721b9119f0d1a7ab7e1859e680
--- /dev/null
+++ b/models/ER.jl
@@ -0,0 +1,54 @@
+
+
+import StaticArrays: SVector, SMatrix, @SMatrix
+
+State = SVector{4, Int}
+Parameters = SVector{3, Real}
+
+d=4
+k=3
+dict_var = Dict("E" => 1, "S" => 2, "ES" => 3, "P" => 4)
+dict_p = Dict("k1" => 1, "k2" => 2, "k3" => 3)
+l_tr = ["R1","R2","R3"]
+p = Parameters(1.0, 1.0, 1.0)
+x0 = State(100, 100, 0, 0)
+t0 = 0.0
+
+function f(xn::State, tn::Real, p::Parameters)
+    a1 = p[1] * xn[1] * xn[2]
+    a2 = p[2] * xn[3]
+    a3 = p[3] * xn[3]
+    l_a = SVector(a1, a2, a3)
+    asum = sum(l_a)
+    l_nu = @SMatrix [-1 1 1;
+                     -1 1 0;
+                     1 -1 -1;
+                     0 0 1]
+    u1, u2 = rand(), rand()
+    tau = - log(u1) / asum
+    b_inf = 0.0
+    b_sup = a1
+    reaction = 0
+    for i = 1:3
+        if b_inf < asum*u2 < b_sup
+            reaction = i
+            break
+        end
+        b_inf += l_a[i]
+        b_sup += l_a[i+1]
+    end
+ 
+    nu = l_nu[:,reaction]
+    xnplus1 = State(xn[1]+nu[1], xn[2]+nu[2], xn[3]+nu[3], xn[4]+nu[4])
+    tnplus1 = tn + tau
+    transition = "R$(reaction)"
+    
+    return xnplus1, tnplus1, transition
+end
+is_absorbing_sir(p::Parameters,xn::State) = 
+    (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) == 0.0
+g = SVector("P")
+
+ER = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,f,is_absorbing_sir; g=g)
+export ER
+
diff --git a/models/SIR.jl b/models/SIR.jl
index de72297c19c29f392692fcee202c303e6cff0dbd..122889bd1310e1cc0b37d72bf2b1a741ca1cd4a5 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -18,9 +18,9 @@ function f(xn::State, tn::Real, p::Parameters)
     l_a = SVector(a1, a2)
     asum = sum(l_a)
     # column-major order
-    l_nu = @SMatrix [-1.0 0.0;
-                     1.0 -1.0;
-                     0.0 1.0]
+    l_nu = @SMatrix [-1 0;
+                     1 -1;
+                     0 1]
     
     u1, u2 = rand(), rand()
     tau = - log(u1) / asum
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index 9aee0767b25e17afcdba9fa691dbb46470330a2a..e85df99ff42f211b4bcd059ccfda906e64a923d5 100644
--- a/tests/run_unit.jl
+++ b/tests/run_unit.jl
@@ -5,5 +5,7 @@ using Test
     @test include("unit/load_model.jl")
     @test include("unit/load_module.jl")
     @test include("unit/simulate_sir.jl")
+    @test include("unit/simulate_sir_bounded.jl")
+    @test include("unit/simulate_er.jl")
 end
 
diff --git a/tests/simulation/sim_er.jl b/tests/simulation/sim_er.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f46aaacbdf2bf477e026638f4834b02eb118387a
--- /dev/null
+++ b/tests/simulation/sim_er.jl
@@ -0,0 +1,12 @@
+
+using MarkovProcesses 
+using PyPlot
+
+load_model("ER")
+
+σ = simulate(ER)
+plt.figure()
+plt.step(σ["times"], σ["P"], "ro--", marker="x", where="post", linewidth=1.0)
+plt.savefig("sim_er.png")
+plt.close()
+
diff --git a/tests/simulation/sim_sir_bounded.jl b/tests/simulation/sim_sir_bounded.jl
new file mode 100644
index 0000000000000000000000000000000000000000..bbca0535d2090101c912ff99167c1d499a7e4f99
--- /dev/null
+++ b/tests/simulation/sim_sir_bounded.jl
@@ -0,0 +1,13 @@
+
+using MarkovProcesses 
+using PyPlot
+
+load_model("SIR")
+SIR.time_bound = 100.0
+
+σ = simulate(SIR)
+plt.figure()
+plt.step(σ["times"], σ["I"], "ro--", marker="x", where="post", linewidth=1.0)
+plt.savefig("sim_sir_bounded.png")
+plt.close()
+
diff --git a/tests/unit/simulate_er.jl b/tests/unit/simulate_er.jl
new file mode 100644
index 0000000000000000000000000000000000000000..32340d1d9fa954c7acc3983dffa70de5fe9265f7
--- /dev/null
+++ b/tests/unit/simulate_er.jl
@@ -0,0 +1,16 @@
+
+using MarkovProcesses
+
+load_model("ER")
+
+σ = simulate(ER)
+
+d1 = Dict("E" => 1, "S"=> 2, "ES" => 3, "P" => 4)
+d2 = Dict("P" => 1)
+
+bool_test = ER.g == ["P"] && ER._g_idx == [4] && 
+            ER.map_var_idx == d1 && 
+            ER._map_obs_var_idx == d2
+
+return bool_test
+
diff --git a/tests/unit/simulate_sir_bounded.jl b/tests/unit/simulate_sir_bounded.jl
new file mode 100644
index 0000000000000000000000000000000000000000..c0b68e6c5b394bec62f3d31b481a7a92e6a86aaf
--- /dev/null
+++ b/tests/unit/simulate_sir_bounded.jl
@@ -0,0 +1,17 @@
+
+using MarkovProcesses
+
+load_model("SIR")
+SIR.time_bound = 100.0
+
+σ = simulate(SIR)
+
+d1 = Dict("S" => 1, "I"=> 2, "R" => 3)
+d2 = Dict("I" => 1)
+
+bool_test = SIR.g == ["I"] && SIR._g_idx == [2] && 
+            SIR.map_var_idx == d1 && 
+            SIR._map_obs_var_idx == d2
+
+return bool_test
+