diff --git a/core/model.jl b/core/model.jl
index 4e2a797dcdcf52a97a8afbd368aaebe47d3b4247..83b53133ce7ac2c1a7e7b2e1e5e16323e93ba481 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -133,7 +133,10 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
     if p != nothing
         p_sim = p
     end
+    x0 = getfield(m, :x0)
+    t0 = getfield(m, :t0)
     time_bound = getfield(m, :time_bound)
+    S0 = init_state(A, x0, t0)
     buffer_size = getfield(m, :buffer_size)
     estim_min_states = getfield(m, :estim_min_states)
     # First alloc
@@ -142,15 +145,15 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
     times = zeros(Float64, estim_min_states)
     transitions = Vector{Transition}(undef, estim_min_states)
     # Initial values
-    for i = eachindex(full_values) full_values[i][1] = getfield(m, :x0)[i] end
-    times[1] = getfield(m, :t0)
+    for i = eachindex(full_values) full_values[i][1] = x0[i] end
+    times[1] = t0
     transitions[1] = nothing
-    S0 = init_state(A, getfield(m, :x0), getfield(m, :t0))
     # Values at time n
     n = 1
-    xn = copy(getfield(m, :x0))
-    tn = getfield(m, :t0) 
+    xn = copy(x0)
+    tn = copy(t0) 
     Sn = copy(S0)
+    next_state!(Sn, A, x0, t0, nothing, S0, x0, p_sim; verbose = verbose)
     isabsorbing::Bool = getfield(m, :isabsorbing)(p_sim,xn)
     isacceptedLHA::Bool = isaccepted(Sn)
     # Alloc of vectors where we stock n+1 values
@@ -269,13 +272,16 @@ function volatile_simulate(product::SynchronizedModel;
     if p != nothing
         p_sim = p
     end
+    x0 = getfield(m, :x0)
+    t0 = getfield(m, :t0)
     time_bound = getfield(m, :time_bound)
-    S0 = init_state(A, getfield(m, :x0), getfield(m, :t0))
+    S0 = init_state(A, x0, t0)
     # Values at time n
     n = 1
-    xn = copy(getfield(m, :x0))
-    tn = getfield(m, :t0) 
+    xn = copy(x0)
+    tn = copy(t0) 
     Sn = copy(S0)
+    next_state!(Sn, A, x0, t0, nothing, S0, x0, p_sim; verbose = verbose)
     isabsorbing::Bool = getfield(m, :isabsorbing)(p_sim,xn)
     isacceptedLHA::Bool = isaccepted(Sn)
     # Alloc of vectors where we stock n+1 values
diff --git a/core/plots.jl b/core/plots.jl
index 2f260ffc833dd50789d9cc6568b14b1e4795e472..7e1a15f940648e5428ef9306e5cd21e07157ef57 100644
--- a/core/plots.jl
+++ b/core/plots.jl
@@ -53,25 +53,27 @@ function plot(σ::AbstractTrajectory, vars::VariableModel...; plot_transitions::
     end
 end
 
-function plot!(A::LHA)
+function plot!(A::LHA; label::String = "")
+    label_region = (label == "") ? "Region LHA" : label
     if A.name in ["F property", "G property"]
         if haskey(A.constants, :x1) && haskey(A.constants, :x2) && haskey(A.constants, :t1) && haskey(A.constants, :t2) 
             x1, x2, t1, t2 = A.constants[:x1], A.constants[:x2], A.constants[:t1], A.constants[:t2] 
-            plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = "Region LHA")
+            plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = label_region)
         end
     end
     if A.name in ["G and F property"]    
         if haskey(A.constants, :x1) && haskey(A.constants, :x2) && haskey(A.constants, :t1) && haskey(A.constants, :t2) 
             x1, x2, t1, t2 = A.constants[:x1], A.constants[:x2], A.constants[:t1], A.constants[:t2] 
-            plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = "Region LHA")
+            plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = label_region)
         end
         if haskey(A.constants, :x3) && haskey(A.constants, :x4) && haskey(A.constants, :t3) && haskey(A.constants, :t4) 
             x3, x4, t3, t4 = A.constants[:x3], A.constants[:x4], A.constants[:t3], A.constants[:t4] 
-            plot!(Shape([(t3,x3), (t3,x4), (t4,x4), (t4,x3), (t3,x3)]), opacity = 0.5, label = "Region LHA")
+            plot!(Shape([(t3,x3), (t3,x4), (t4,x4), (t4,x3), (t3,x3)]), opacity = 0.5, label = label_region)
         end
     end
+    label_region = (label == "") ? "L, H" : label
     if A.name == "Period"
-        hline!([A.constants[:L], A.constants[:H]], label = "L, H", linestyle = :dot)
+        hline!([A.constants[:L], A.constants[:H]], label = label_region, linestyle = :dot)
     end
 end
 
diff --git a/models/poisson.jl b/models/poisson.jl
index c0133aeed16d3497158b1a695910e532239cc320..02de0954abd42098742319c99e43bc3ba97a227b 100644
--- a/models/poisson.jl
+++ b/models/poisson.jl
@@ -10,8 +10,8 @@ l_tr_poisson = [:R]
 p_poisson = [5.0]
 x0_poisson = [0]
 t0_poisson = 0.0
-function poisson_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
-                    xn::Vector{Int}, tn::Float64, p::Vector{Float64})
+@everywhere function poisson_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
+                                xn::Vector{Int}, tn::Float64, p::Vector{Float64})
     
     u1 = rand()
     tau = (-log(u1)/p[1])
@@ -19,12 +19,12 @@ function poisson_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Tra
     l_t[1] = tn + tau
     l_tr[1] = :R
 end
-isabsorbing_poisson(p::Vector{Float64}, xn::Vector{Int}) = p[1] === 0.0
+@everywhere isabsorbing_poisson(p::Vector{Float64}, xn::Vector{Int}) = p[1] === 0.0
 g_poisson = [:N]
 
 poisson = ContinuousTimeModel(d,k,dict_var_poisson,dict_p_poisson,l_tr_poisson,
                                   p_poisson,x0_poisson,t0_poisson,
-                                  poisson_f!,isabsorbing_poisson; 
+                                  getfield(Main, :poisson_f!), getfield(Main, :isabsorbing_poisson); 
                                   g=g_poisson, time_bound=1.0, name="Poisson process pkg")
 function create_poisson(new_p::Vector{Float64})
     poisson_new = deepcopy(poisson)
diff --git a/tests/automata/absorbing_x0_p.jl b/tests/automata/absorbing_x0_p.jl
index 2fd9843c40f9a9e2822d983dca26814a4b058d60..14c62197a36c366ef9f5c78a98cd16198e2b3b6f 100644
--- a/tests/automata/absorbing_x0_p.jl
+++ b/tests/automata/absorbing_x0_p.jl
@@ -1,6 +1,8 @@
 
 using MarkovProcesses
 
+test_all = true
+
 load_model("ER")
 observe_all!(ER)
 load_automaton("automaton_F")
@@ -14,7 +16,6 @@ A_G_F_R6 = create_automaton_G_and_F(ER, x1, x2, t1, t2, :E,
                                     x3, x4, t3, t4, :P) 
 
 set_param!(ER, [0.0,0.0,0.0])
-test_all = true
 for A in [A_F_R1, A_G_R5, A_G_F_R6]
     local sync = A*ER
     local σ, σ2 = simulate(sync), simulate(sync)
@@ -23,12 +24,18 @@ end
 
 set_param!(ER, [1.0,1.0,1.0])
 set_x0!(ER, [0, 0, 0, 5])
-test_all = true
 for A in [A_F_R1, A_G_R5, A_G_F_R6]
     local sync = A*ER
     local σ, σ2 = simulate(sync), simulate(sync)
     global test_all = test_all && isaccepted(σ) && isaccepted(σ2)
 end
 
+load_model("poisson")
+set_param!(poisson, [0.0])
+aut = create_automaton_F(poisson, 5.0, 7.0, 0.75, 1.0, :N)
+sync_poisson = aut * poisson 
+test_all = test_all && simulate(sync_poisson).state_lha_end[:d] > 0
+test_all = test_all && volatile_simulate(sync_poisson)[:d] > 0
+
 return true