diff --git a/core/biochemical_network.jl b/core/biochemical_network.jl
index c2592c4a9ef031506cab2fd2aaa6868d8774e670..0866b412ee7e70530cfae8615de7030abfbdda6c 100644
--- a/core/biochemical_network.jl
+++ b/core/biochemical_network.jl
@@ -35,7 +35,8 @@ function get_str_propensity(propensity::Symbol, dict_species::Dict, dict_params:
     return str_propensity
 end
 
-macro biochemical_network(expr_name,expr_network)
+macro biochemical_network(expr_network,expr_name...)
+    model_name = isempty(expr_name) ? "Unnamed macro generated" : expr_name[1]
     transitions = String[]
     dict_species = Dict{String,Int}()
     dict_params = Dict{String,Int}()
@@ -103,7 +104,7 @@ macro biochemical_network(expr_name,expr_network)
     # Let's write some lines that creates the function f! (step of a simulation) for this biochemical network
     nbr_rand = rand(1:1000)
     nbr_reactions = length(list_expr_reactions)
-    basename_func = "$(replace(expr_name, ' '=>'_'))_$(nbr_rand)"
+    basename_func = "$(replace(model_name, ' '=>'_'))_$(nbr_rand)"
     expr_model_f! = "function $(basename_func)_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, xn::Vector{Int}, tn::Float64, p::Vector{Float64})\n\t"
     # Computation of nu and propensity functions in f!
     str_l_a = "l_a = ("
@@ -178,6 +179,7 @@ macro biochemical_network(expr_name,expr_network)
     model_f! = eval(Meta.parse(expr_model_f!))
     model_isabsorbing = eval(Meta.parse(expr_model_isabsorbing))
     return :(ContinuousTimeModel($dim_state, $dim_params, $dict_species, $dict_params, $transitions, 
-                                 $(zeros(dim_params)), $(zeros(Int, dim_state)), 0.0, $model_f!, $model_isabsorbing; g = $list_species))
+                                 $(zeros(dim_params)), $(zeros(Int, dim_state)), 0.0, $model_f!, $model_isabsorbing; 
+                                 g = $list_species, name=$model_name))
 end
 
diff --git a/core/common.jl b/core/common.jl
index e283f8330e438d2c63827436463205b358c51183..718d60936ccaf328acfa63455c343efb1ba679cf 100644
--- a/core/common.jl
+++ b/core/common.jl
@@ -10,6 +10,7 @@ const Location = String
 const VariableAutomaton = String
 
 mutable struct ContinuousTimeModel <: Model
+    name::String
     dim_state::Int # state space dim
     dim_params::Int # parameter space dim
     map_var_idx::Dict{String,Int} # maps variable str to index in the state space
@@ -85,7 +86,8 @@ end
 function ContinuousTimeModel(dim_state::Int, dim_params::Int, map_var_idx::Dict, map_param_idx::Dict, transitions::Vector{String}, 
               p::Vector{Float64}, x0::Vector{Int}, t0::Float64, 
               f!::Function, isabsorbing::Function; 
-              g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10, estim_min_states::Int = 50)
+              g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, 
+              buffer_size::Int = 10, estim_min_states::Int = 50, name::String = "Unnamed")
     dim_obs_state = length(g)
     _map_obs_var_idx = Dict()
     _g_idx = Vector{Int}(undef, dim_obs_state)
@@ -100,7 +102,8 @@ function ContinuousTimeModel(dim_state::Int, dim_params::Int, map_var_idx::Dict,
     if length(methods(isabsorbing)) >= 2
         @warn "You have possibly redefined a function Model.isabsorbing used in a previously instantiated model."
     end
-    new_model = ContinuousTimeModel(dim_state, dim_params, map_var_idx, _map_obs_var_idx, map_param_idx, transitions, p, x0, t0, f!, g, _g_idx, isabsorbing, time_bound, buffer_size, estim_min_states)
+    new_model = ContinuousTimeModel(name, dim_state, dim_params, map_var_idx, _map_obs_var_idx, map_param_idx, transitions, 
+                                    p, x0, t0, f!, g, _g_idx, isabsorbing, time_bound, buffer_size, estim_min_states)
     @assert check_consistency(new_model)
     return new_model
 end
diff --git a/core/model.jl b/core/model.jl
index 26136ba705e92e7eee7c74187f893d18cdaecab8..13809b739633b01e032150bf402dc3656142e653 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -364,6 +364,23 @@ function distribute_prob_accept_lha(sm::SynchronizedModel, nbr_sim::Int)
     return sum_val / nbr_sim
 end
 
+function Base.show(io::IO, m::ContinuousTimeModel)
+    print(io, "$(m.name) model\n")
+    print(io, "- variables :\n")
+    for (var, idx) in m.map_var_idx
+        print(io, "* $var (idx = $idx in state space)\n")
+    end
+    print(io, "- parameters :\n")
+    for (param, idx) in m.map_param_idx
+        print(io, "* $param (idx = $idx in parameter space)\n")
+    end
+    print(io, "- transitions : $(m.transitions)\n")
+    print(io, "- observed variables :\n")
+    for i in eachindex(m.g)
+        print(io, "* $(m.g[i]) (idx = $i in observed state space)\n")
+    end
+end
+
 isbounded(m::Model) = get_proba_model(m).time_bound < Inf
 function check_consistency(m::ContinuousTimeModel) 
     @assert m.dim_state == length(m.map_var_idx) 
diff --git a/models/ER.jl b/models/ER.jl
index 11e703b2b6f81c79fc4c81b7ac27fd190ce586b5..2af5b5f409ec0155b6aa095a34d530e105335b4c 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -52,7 +52,7 @@ isabsorbing_ER(p::Vector{Float64},xn::Vector{Int}) =
     (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)
+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")
 
 function create_ER(new_p::Vector{Float64})
     ER_new = deepcopy(ER)
diff --git a/models/SIR.jl b/models/SIR.jl
index 93fbbb920ef26eedc60d3abcf03c1d3dd7dedf10..592835f1d44553fa8b858ba686ee59fbbdfe6876 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)
+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")
 
 function create_SIR(new_p::Vector{Float64})
     SIR_new = deepcopy(SIR)
diff --git a/models/SIR_tauleap.jl b/models/SIR_tauleap.jl
index 6dac452d92feb7e4f0478010da969f0d6787abf8..3d9cb6e7e832be5b3a92f59424d6f0c01539593c 100644
--- a/models/SIR_tauleap.jl
+++ b/models/SIR_tauleap.jl
@@ -37,7 +37,7 @@ g_SIR_tauleap = ["I"]
 
 SIR_tauleap = ContinuousTimeModel(d,k,dict_var_SIR_tauleap,dict_p_SIR_tauleap,l_tr_SIR_tauleap,
                                   p_SIR_tauleap,x0_SIR_tauleap,t0_SIR_tauleap,
-                                  SIR_tauleap_f!,isabsorbing_SIR_tauleap; g=g_SIR_tauleap)
+                                  SIR_tauleap_f!,isabsorbing_SIR_tauleap; g=g_SIR_tauleap, name="SIR tauleap pkg")
 function create_SIR_tauleap(new_p::Vector{Float64})
     SIR_tauleap_new = deepcopy(SIR_tauleap)
     @assert length(SIR_tauleap_new.p) == length(new_p)
diff --git a/models/poisson.jl b/models/poisson.jl
index f5991aecc075af058d83a5c09492153f69761870..db080272085cd79e5f9648ecaf8bcdc4d23d8cf3 100644
--- a/models/poisson.jl
+++ b/models/poisson.jl
@@ -24,7 +24,8 @@ 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; g=g_poisson, time_bound=1.0)
+                                  poisson_f!,isabsorbing_poisson; 
+                                  g=g_poisson, time_bound=1.0, name="Poisson process pkg")
 function create_poisson(new_p::Vector{Float64})
     poisson_new = deepcopy(poisson)
     @assert length(poisson_new.p) == length(new_p)
diff --git a/tests/unit/macro_biochemical_network.jl b/tests/unit/macro_biochemical_network.jl
index 5fab6c55eddbb7f9d7a8b1f3c0d7b6f76772dbe0..13295975103277cbb929fe7322bd0167ccaf6d55 100644
--- a/tests/unit/macro_biochemical_network.jl
+++ b/tests/unit/macro_biochemical_network.jl
@@ -1,18 +1,25 @@
 
 using MarkovProcesses 
 
-model_SIR = @biochemical_network "SIR" begin
+model_SIR = @biochemical_network begin
     R1: (S+I => 2I, ki*S*I)
     R2: (I => R, kr*I)
-end
+end "SIR"
 set_x0!(model_SIR, [95,5,0])
 set_param!(model_SIR, [0.012, 0.05])
 
-model_ER = @biochemical_network "ER" begin
+model_unnamed_SIR = @biochemical_network begin
+    R1: (S+I => 2I, ki*S*I)
+    R2: (I => R, kr*I)
+end
+set_x0!(model_unnamed_SIR, [95,5,0])
+set_param!(model_unnamed_SIR, [0.012, 0.05])
+
+model_ER = @biochemical_network begin
     R1: (E+S => ES, k1*E*S)
     R2: (ES => E+S, k2*ES)
     R3: (ES => E+P, k3*ES)
-end
+end "ER"
 set_x0!(model_ER, [100,100,0,0])
 set_param!(model_ER, [1.0,1.0,1.0])