diff --git a/core/plots.jl b/core/plots.jl
index caf17e75b48d0825b1f68856a40a31df7e80d1e9..2f260ffc833dd50789d9cc6568b14b1e4795e472 100644
--- a/core/plots.jl
+++ b/core/plots.jl
@@ -1,6 +1,6 @@
 
 import Plots: plot, plot!, scatter!, hline!, Shape, text
-import Plots: current, palette, display, png, close
+import Plots: current, palette, display, png, close, savefig
 
 """
     `plot(σ, var...; plot_transitions=false)`
@@ -27,14 +27,14 @@ function plot(σ::AbstractTrajectory, vars::VariableModel...; plot_transitions::
         @assert var in get_obs_var(σ) "Variable $var is not observed." 
         plot!(p, times(σ), σ[var], 
               xlabel = "Time", ylabel = "Number of species",
-              label = var,
+              label = "$var",
               linetype=:steppost)
     end
     if plot_transitions
         for (i, var) in enumerate(to_plot)
             for tr in l_tr
                 idx_tr = findall(x->x==tr, transitions(σ))
-                label = (tr == nothing || i > 1) ? "" : tr
+                label = (tr == nothing || i > 1) ? "" : "$tr"
                 alpha = (tr == nothing) ? 0.0 : 0.5
                 scatter!(p, times(σ)[idx_tr], σ[var][idx_tr], label=label, 
                          markershape=:cross, markeralpha=alpha, 
@@ -49,7 +49,7 @@ function plot(σ::AbstractTrajectory, vars::VariableModel...; plot_transitions::
     if filename == ""
         display(p)
     else
-        png(p, filename)
+        savefig(p, filename)
     end
 end
 
@@ -124,7 +124,7 @@ function plot_periodic_trajectory(A::LHA, σ::SynchronizedTrajectory, sym_obs::S
     if filename == ""
         display(p)
     else
-        png(p, filename)
+        savefig(p, filename)
     end
 end
 
diff --git a/examples/scripts/doping_3way_oscillator.jl b/examples/scripts/doping_3way_oscillator.jl
index cbfba14af7fad35a5e19e77e5207969620e3c996..bfd1167a46bd9ba0a28457455a1796999bf91b63 100644
--- a/examples/scripts/doping_3way_oscillator.jl
+++ b/examples/scripts/doping_3way_oscillator.jl
@@ -10,7 +10,7 @@ sync_doping = doping_3way_oscillator * A_per
 set_time_bound!(sync_doping, 0.1)
 set_x0!(doping_3way_oscillator, [:A, :B, :C, :DA, :DB, :DC], [333, 333, 333, 10, 10, 10])
 σ = simulate(sync_doping)
-plot(σ; A = A_per, filename = "traj_full.png")
-plot(σ, :A; A = A_per, filename = "traj_A.png")
-plot_periodic_trajectory(A_per, σ, :A, filename = "traj_automaton.png")
+plot(σ; A = A_per, filename = "traj_full.svg")
+plot(σ, :A; A = A_per, filename = "traj_A.svg")
+plot_periodic_trajectory(A_per, σ, :A, filename = "traj_automaton.svg")
 
diff --git a/examples/scripts/enzym_1d.jl b/examples/scripts/enzym_1d.jl
index 6e385584ac5e698c04319430d3799b22a2b8bd19..bc6196efca080b69dc4f5311df94d0d7ccd23f8d 100644
--- a/examples/scripts/enzym_1d.jl
+++ b/examples/scripts/enzym_1d.jl
@@ -45,7 +45,7 @@ samples_weights = r.weights
 
 # Histogram
 histogram(samples_abc_post, weights = r.weights, normalize = :density)
-png(path_results * "histogram.png")
+savefig(path_results * "histogram.svg")
 
 ## Satisfaction function
 
@@ -92,5 +92,5 @@ y_MC = readdlm("/home/moud/plot_R1-3/estim_MC/$(exp)/satisfaction_func.csv", ','
 inf_x, sup_x = 0.0, 100.0
 x_MC = inf_x:((sup_x-inf_x)/(length(y_MC)-1)):sup_x
 plot!(x_MC, y_MC, label = "MC spf")
-png(path_results * "satisfaction_prob_function.png")
+savefig(path_results * "satisfaction_prob_function.svg")
 
diff --git a/examples/scripts/enzym_2d.jl b/examples/scripts/enzym_2d.jl
index dc74706f672c2f91df80744e4ee8fa10695a74a9..1ebb12dc51ad49081250291ea6eb83a637919827 100644
--- a/examples/scripts/enzym_2d.jl
+++ b/examples/scripts/enzym_2d.jl
@@ -47,7 +47,7 @@ samples_weights = r.weights
 
 # Histogram
 histogram2d(samples_abc_post[1,:], samples_abc_post[2,:], weights = samples_weights, normalize = :density)
-png(path_results * "histogram.png")
+savefig(path_results * "histogram.svg")
 =#
 ## Satisfaction function
 
@@ -88,11 +88,11 @@ xaxis = 0:0.1:1.5
 yaxis = 0:1.0:100.0
 p = plot(title = "Multivariate KDE", dpi = 480, background_color_legend = :transparent)
 plot!(p, xaxis, yaxis, prob_func, st = :surface, c = :coolwarm, camera = (30, 45))
-png(path_results * "estim_abc_satisfaction_prob_function.png")
+savefig(path_results * "estim_abc_satisfaction_prob_function.svg")
 x_MC = readdlm("/home/moud/results_last_automata/estim_satisfaction_func_MC/$(exp)/grid_X.csv", ',')
 y_MC = readdlm("/home/moud/results_last_automata/estim_satisfaction_func_MC/$(exp)/grid_Y.csv", ',')
 z_MC = readdlm("/home/moud/results_last_automata/estim_satisfaction_func_MC/$(exp)/satisfaction_func.csv", ',')
 p = plot(title = "ABC MC", dpi = 480, background_color_legend = :transparent)
 plot!(p, [x_MC...], [y_MC...], [z_MC...], st = :surface, c = :coolwarm, camera = (30, 45), label = "MC spf")
-png(path_results * "estim_MC_satisfaction_prob_function.png")
+savefig(path_results * "estim_MC_satisfaction_prob_function.svg")
 
diff --git a/tests/automaton_abc/R1.jl b/tests/automaton_abc/R1.jl
index cde4a9f19dfb9a9ae2dc84dcebbde5b579f558f0..ee76c2ef89641fc011f9b3dbb119b7367937749a 100644
--- a/tests/automaton_abc/R1.jl
+++ b/tests/automaton_abc/R1.jl
@@ -20,5 +20,5 @@ return test
 
 #using Plots
 #histogram(r.mat_p_end', weights = r.weights, normalize = :density)
-#png("R1_hist.png")
+#savefig("R1_hist.svg")
 
diff --git a/tests/automaton_abc/R2.jl b/tests/automaton_abc/R2.jl
index ac624b37f349561e8edc10dacac3165876652a6d..369050f63ff2f0946ba26fec97ce69c04599d900 100644
--- a/tests/automaton_abc/R2.jl
+++ b/tests/automaton_abc/R2.jl
@@ -13,5 +13,5 @@ pm_sync_ER = ParametricModel(sync_ER, (:k3, Uniform(0.0, 100.0)))
 @show r.nbr_sim
 
 histogram(r.mat_p_end', bins = :sqrt, weights = r.weights, normalize = :density, xlims = (0.0, 100.0))
-png("R2_hist.png")
+savefig("R2_hist.svg")
 
diff --git a/tests/automaton_abc/R3.jl b/tests/automaton_abc/R3.jl
index bc41c684a53242cd64077865669fd549e9ec6a95..b263f94bd9bb2e7621550831ee5faa0ce658e61d 100644
--- a/tests/automaton_abc/R3.jl
+++ b/tests/automaton_abc/R3.jl
@@ -11,5 +11,5 @@ pm_sync_ER = ParametricModel(sync_ER, (:k3, Uniform(0.0, 100.0)))
 r = automaton_abc(pm_sync_ER; nbr_particles = 1000)
 
 histogram(r.mat_p_end', weights = r.weights, normalize = :density, xlims = (0.0, 100.0))
-png("R3_hist.png")
+savefig("R3_hist.svg")
 
diff --git a/tests/automaton_abc/R6.jl b/tests/automaton_abc/R6.jl
new file mode 100644
index 0000000000000000000000000000000000000000..b8c8cef09d5c80267cf9bf5646a68ec8f7ba8338
--- /dev/null
+++ b/tests/automaton_abc/R6.jl
@@ -0,0 +1,16 @@
+
+using MarkovProcesses
+
+load_model("ER")
+observe_all!(ER)
+load_automaton("automaton_G_and_F")
+A_G_F_R6 = create_automaton_G_and_F(ER, 50.0, 100.0, 0.0, 0.8, :E, 30.0, 100.0, 0.8, 0.9, :P)
+sync_ER = A_G_F_R6 * ER
+pm_sync_ER = ParametricModel(sync_ER, (:k1, Uniform(0.0,100.0)), (:k2, Uniform(0.0,100.0)))
+
+r = automaton_abc(pm_sync_ER; alpha=0.2)
+
+using Plots
+histogram2d(r.mat_p_end[1,:], r.mat_p_end[2,:], bins=50)
+savefig("R6_hist.svg")
+
diff --git a/tests/automaton_abc/distributed_R1.jl b/tests/automaton_abc/distributed_R1.jl
index f10b7affc0cf981118958238bd8eeced531fd3d9..5fda0e80c3ac1a5068dd964100a06ea8fc6bb859 100644
--- a/tests/automaton_abc/distributed_R1.jl
+++ b/tests/automaton_abc/distributed_R1.jl
@@ -34,5 +34,5 @@ test = size(r.mat_p_end)[1] == pm_sync_ER.df &&
 return test
 
 #histogram(r.mat_p_end', weights = r.weights, normalize = :density)
-#png("R1_hist.png")
+#savefig("R1_hist.svg")
 
diff --git a/tests/simulation/sim_er_row_buffer_bounded.jl b/tests/simulation/sim_er_row_buffer_bounded.jl
index 940fdce7c5579a2ea1cdb4534a7bb4a09edb4907..c92dd9e82f2958db59cde9452419e5f67ef32263 100644
--- a/tests/simulation/sim_er_row_buffer_bounded.jl
+++ b/tests/simulation/sim_er_row_buffer_bounded.jl
@@ -9,7 +9,7 @@ ER_row_buffer.time_bound = 10.0
 σ = _simulate_row_buffer(ER_row_buffer)
 plt.figure()
 plt.step(times(σ), _get_values_row(σ,:P), "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_er_row_buffer_bounded.png")
+plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_er_row_buffer_bounded.svg")
 plt.close()
 
 return true
diff --git a/tests/simulation/sim_pm_er.jl b/tests/simulation/sim_pm_er.jl
index b5e3e5a075c75a2fe2962caff77407bf2171b1e1..492d46438f3a36e042329708996075bea6b9d62b 100644
--- a/tests/simulation/sim_pm_er.jl
+++ b/tests/simulation/sim_pm_er.jl
@@ -9,7 +9,7 @@ pm_ER = ParametricModel(ER, (:k1, Uniform(0.0,100.0)), (:k2, Uniform(0.0,100.0))
 prior_p = [0.2, 40.0]
 
 σ = simulate(pm_ER, prior_p)
-MarkovProcesses.plot(σ; filename = get_module_path() * "/tests/simulation/res_pics/sim_pm_er_long.png")
+MarkovProcesses.plot(σ; filename = get_module_path() * "/tests/simulation/res_pics/sim_pm_er_long.svg")
 
 return true
 
diff --git a/tests/simulation/sim_pm_sync_er.jl b/tests/simulation/sim_pm_sync_er.jl
index 56e92b64c777c1cb3652b3b80b77d7bb34e9a365..19559d7609459417d04b5f0ac69ce60eea56ade8 100644
--- a/tests/simulation/sim_pm_sync_er.jl
+++ b/tests/simulation/sim_pm_sync_er.jl
@@ -11,7 +11,7 @@ pm_sync_ER = ParametricModel(ER*A_F, (:k1, Uniform(0.0,100.0)), (:k2, Uniform(0.
 prior_p = [0.2, 40.0]
 
 σ = simulate(pm_sync_ER, prior_p)
-MarkovProcesses.plot(σ; filename = get_module_path() * "/tests/simulation/res_pics/sim_pm_sync_er_long.png")
+MarkovProcesses.plot(σ; filename = get_module_path() * "/tests/simulation/res_pics/sim_pm_sync_er_long.svg")
 
 return true
 
diff --git a/tests/simulation/sim_poisson.jl b/tests/simulation/sim_poisson.jl
index 40aa61dc92be289aa9eb1e6bf52b03ea942d89f8..ca050a56b874eee02a610ef67480c1ecd14c9995 100644
--- a/tests/simulation/sim_poisson.jl
+++ b/tests/simulation/sim_poisson.jl
@@ -4,7 +4,7 @@ load_plots()
 
 load_model("poisson")
 σ = simulate(poisson)
-MarkovProcesses.plot(σ; filename = get_module_path() * "/tests/simulation/res_pics/sim_poisson.png")
+MarkovProcesses.plot(σ; filename = get_module_path() * "/tests/simulation/res_pics/sim_poisson.svg")
 
 return true
 
diff --git a/tests/simulation/sim_sir.jl b/tests/simulation/sim_sir.jl
index 9ee9755b93e5bc32185ad46ff084df8dab8702b7..fe2740a63c5a057d7f832dc8c2aac027cbd2a6fe 100644
--- a/tests/simulation/sim_sir.jl
+++ b/tests/simulation/sim_sir.jl
@@ -7,7 +7,7 @@ load_model("SIR")
 σ = simulate(SIR)
 plt.figure()
 plt.step(times(σ), σ[:I], "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir.png")
+plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir.svg")
 plt.close()
 
 return true
diff --git a/tests/simulation/sim_sir_bounded.jl b/tests/simulation/sim_sir_bounded.jl
index 19a3a121e22b95ea04a0e126eefdad5dc5f7bcc7..fe3dcc76264b74fcce2b83bc51b11ed7610508fc 100644
--- a/tests/simulation/sim_sir_bounded.jl
+++ b/tests/simulation/sim_sir_bounded.jl
@@ -8,7 +8,7 @@ SIR.time_bound = 100.0
 σ = simulate(SIR)
 plt.figure()
 plt.step(times(σ), σ[:I], "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir_bounded.png")
+plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir_bounded.svg")
 plt.close()
 
 return true
diff --git a/tests/simulation/sim_sir_col_buffer_bounded.jl b/tests/simulation/sim_sir_col_buffer_bounded.jl
index 472195d2dffeaecbf16ba650407699bbff819a1c..d68f9901db947487057269c01e643af529e6fb6b 100644
--- a/tests/simulation/sim_sir_col_buffer_bounded.jl
+++ b/tests/simulation/sim_sir_col_buffer_bounded.jl
@@ -9,7 +9,7 @@ SIR_col_buffer.time_bound = 100.0
 σ = _simulate_col_buffer(SIR_col_buffer)
 plt.figure()
 plt.step(times(σ), _get_values_col(σ,:I), "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir_col_buffer_bounded.png")
+plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir_col_buffer_bounded.svg")
 plt.close()
 
 return true
diff --git a/tests/simulation/sim_sir_row_buffer_bounded.jl b/tests/simulation/sim_sir_row_buffer_bounded.jl
index 034862d0801dcc3d20239fcc178f48bcec1a1616..d10633dc4f8131fd1224e99e5e27ed6768a4c451 100644
--- a/tests/simulation/sim_sir_row_buffer_bounded.jl
+++ b/tests/simulation/sim_sir_row_buffer_bounded.jl
@@ -9,7 +9,7 @@ SIR_row_buffer.time_bound = 100.0
 σ = _simulate_row_buffer(SIR_row_buffer)
 plt.figure()
 plt.step(times(σ), _get_values_row(σ,:I), "ro--", marker="x", where="post", linewidth=1.0)
-plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir_row_buffer_bounded.png")
+plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir_row_buffer_bounded.svg")
 plt.close()
 
 return true
diff --git a/tests/simulation/sim_sir_tauleap.jl b/tests/simulation/sim_sir_tauleap.jl
index a1fa7cf4bf2a337d14d8d72a817f5b7974f9c53d..3b22d6db2404f5b3d68beb2096968b029fc02b47 100644
--- a/tests/simulation/sim_sir_tauleap.jl
+++ b/tests/simulation/sim_sir_tauleap.jl
@@ -4,7 +4,7 @@ load_plots()
 
 load_model("SIR_tauleap")
 σ = simulate(SIR_tauleap)
-MarkovProcesses.plot(σ; filename = get_module_path() * "/tests/simulation/res_pics/sim_sir_tauleap.png")
+MarkovProcesses.plot(σ; filename = get_module_path() * "/tests/simulation/res_pics/sim_sir_tauleap.svg")
 
 return true