diff --git a/core/lha.jl b/core/lha.jl
index 777f19713e7bb11f2ec5a6e43b9bf7d1a5ed23c1..245d6079a1e646a8949d630ba2b23f76311f58c6 100644
--- a/core/lha.jl
+++ b/core/lha.jl
@@ -117,7 +117,7 @@ end
 
 function next_state!(Snplus1::StateLHA, A::LHA, 
                      xnplus1::Vector{Int}, tnplus1::Float64, tr_nplus1::Transition, 
-                     Sn::StateLHA; verbose::Bool = false)
+                     Sn::StateLHA, xn::Vector{Int}; verbose::Bool = false)
     # En fait d'apres observation de Cosmos, après qu'on ait lu la transition on devrait stop.
     edge_candidates = Vector{Edge}(undef, 2)
     first_round::Bool = true
@@ -139,14 +139,14 @@ function next_state!(Snplus1::StateLHA, A::LHA,
         current_loc = getfield(Snplus1, :loc)
         edges_from_current_loc = getfield(A, :map_edges)[current_loc]
         # Save all edges that satisfies transition predicate (asynchronous ones)
-        nbr_candidates = _find_edge_candidates!(edge_candidates, edges_from_current_loc, Snplus1, constants, xnplus1, true)
+        nbr_candidates = _find_edge_candidates!(edge_candidates, edges_from_current_loc, Snplus1, constants, xn, true)
         # Search the one we must chose, here the event is nothing because 
         # we're not processing yet the next event
         ind_edge, detected_event = _get_edge_index(edge_candidates, nbr_candidates, detected_event, nothing)
         # Update the state with the chosen one (if it exists)
         # Should be xn here
         if ind_edge > 0
-            getfield(edge_candidates[ind_edge], :update_state!)(Snplus1, constants, xnplus1)
+            getfield(edge_candidates[ind_edge], :update_state!)(Snplus1, constants, xn)
         end
         first_round = false
         if verbose
@@ -240,8 +240,8 @@ function read_trajectory(A::LHA, σ::Trajectory; verbose = false)
     Snplus1 = copy(Sn)
     if verbose println("Init: ") end
     if verbose @show Sn end
-    for n in 1:length_states(σ)
-        next_state!(Snplus1, A_new, σ[n], l_t[n], l_tr[n], Sn; verbose = verbose)
+    for n in 2:length_states(σ)
+        next_state!(Snplus1, A_new, σ[n], l_t[n], l_tr[n], Sn, σ[n-1]; verbose = verbose)
         copyto!(Sn, Snplus1)
         if Snplus1.loc in A_new.locations_final 
             break 
diff --git a/core/model.jl b/core/model.jl
index 7cc3b2dc75731f33d47393096e30504a283bcaff..f4a117b3e5d94563ed85be86df0ce763522a7262 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -166,7 +166,7 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
             _finish_bounded_trajectory!(values, times, transitions, time_bound)
         end
         if isabsorbing && !isacceptedLHA
-            next_state!(Snplus1, A, xn, time_bound, nothing, Sn; verbose = verbose)
+            next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn; verbose = verbose)
             copyto!(Sn, Snplus1)
         end
         return SynchronizedTrajectory(Sn, product, values, times, transitions)
@@ -174,15 +174,16 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
     # First we fill the allocated array
     for i = 2:estim_min_states
         getfield(m, :f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
-        tn = l_t[1]
-        if tn > time_bound || vec_x == xn
+        if l_t[1] > time_bound || vec_x == xn
+            tn = l_t[1]
             isabsorbing = (vec_x == xn)
             break
         end
         n += 1
+        next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn; verbose = verbose)
         copyto!(xn, vec_x)
+        tn = l_t[1]
         tr_n = l_tr[1]
-        next_state!(Snplus1, A, xn, tn, tr_n, Sn; verbose = verbose)
         _update_values!(full_values, times, transitions, xn, tn, tr_n, i)
         copyto!(Sn, Snplus1)
         isacceptedLHA = isaccepted(Snplus1)
@@ -198,7 +199,7 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
             _finish_bounded_trajectory!(values, times, transitions, time_bound)
         end
         if isabsorbing && !isacceptedLHA
-            next_state!(Snplus1, A, xn, time_bound, nothing, Sn; verbose = verbose)
+            next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn; verbose = verbose)
             copyto!(Sn, Snplus1)
         end
         return SynchronizedTrajectory(Sn, product, values, times, transitions)
@@ -212,8 +213,8 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
         while i < buffer_size
             i += 1
             getfield(m, :f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
-            tn = l_t[1]
-            if tn > time_bound
+            if l_t[1] > time_bound
+                tn = l_t[1]
                 i -= 1
                 break
             end
@@ -222,9 +223,10 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
                 i -= 1
                 break
             end
+            next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn; verbose = verbose)
             copyto!(xn, vec_x)
+            tn = l_t[1]
             tr_n = l_tr[1]
-            next_state!(Snplus1, A, xn, tn, tr_n, Sn; verbose = verbose)
             _update_values!(full_values, times, transitions, 
                             xn, tn, tr_n, estim_min_states+size_tmp+i)
             copyto!(Sn, Snplus1)
@@ -247,7 +249,7 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
         _finish_bounded_trajectory!(values, times, transitions, time_bound)
     end
     if isabsorbing && !isacceptedLHA
-        next_state!(Snplus1, A, xn, time_bound, nothing, Sn; verbose = verbose)
+        next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn; verbose = verbose)
         copyto!(Sn, Snplus1)
     end
     return SynchronizedTrajectory(Sn, product, values, times, transitions)
@@ -284,15 +286,15 @@ function volatile_simulate(product::SynchronizedModel;
     # If x0 is absorbing
     if isabsorbing || isacceptedLHA 
         if !isacceptedLHA
-            next_state!(Snplus1, A, xn, time_bound, nothing, Sn; verbose = verbose)
+            next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn; verbose = verbose)
             copyto!(Sn, Snplus1)
         end
         return Sn
     end
     while !isabsorbing && tn <= time_bound && !isacceptedLHA
         getfield(m, :f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
-        tn = l_t[1]
-        if tn > time_bound
+        if l_t[1] > time_bound
+            tn = l_t[1]
             i -= 1
             break
         end
@@ -300,15 +302,16 @@ function volatile_simulate(product::SynchronizedModel;
             isabsorbing = true
             break
         end
+        next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn; verbose = verbose)
         copyto!(xn, vec_x)
+        tn = l_t[1]
         tr_n = l_tr[1]
-        next_state!(Snplus1, A, xn, tn, tr_n, Sn; verbose = verbose)
         copyto!(Sn, Snplus1)
         isacceptedLHA = isaccepted(Snplus1)
         n += 1
     end
     if isabsorbing && !isacceptedLHA
-        next_state!(Snplus1, A, xn, time_bound, nothing, Sn; verbose = verbose)
+        next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn; verbose = verbose)
         copyto!(Sn, Snplus1)
     end
     return Sn