diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index 36b234fd3ff565fd87ccb85be62813ad12beb709..f34499975a778fbf38636c71f6455039d82037b1 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -40,11 +40,15 @@ export load_model, load_automaton, load_plots
 # Algorithms
 export automaton_abc, abc_smc
 
+# About biochemical networks
+export @biochemical_network
+
 include("common.jl")
 include("trajectory.jl")
 include("lha.jl")
 include("model.jl")
 include("utils.jl")
+include("biochemical_network.jl")
 include("../algorithms/abc_smc.jl")
 
 end
diff --git a/models/ER.jl b/models/ER.jl
index 971c53d8c88dc5e62cf3e27d885eebc851644ef4..f8a6fe7f54ca43ff52b17f8dc806b412e9ad0356 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -14,13 +14,13 @@ function ER_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{No
     a1 = p[1] * xn[1] * xn[2]
     a2 = p[2] * xn[3]
     a3 = p[3] * xn[3]
-    l_a = SVector(a1, a2, a3)
+    l_a = (a1, a2, a3)
     asum = sum(l_a)
-    nu_1 = SVector(-1, -1, 1, 0)
-    nu_2 = SVector(1, 1, -1, 0)
-    nu_3 = SVector(1, 0, -1, 1) 
-    l_nu = SVector(nu_1, nu_2, nu_3)
-    l_str_R = SVector("R1", "R2", "R3")
+    nu_1 = (-1, -1, 1, 0)
+    nu_2 = (1, 1, -1, 0)
+    nu_3 = (1, 0, -1, 1) 
+    l_nu = (nu_1, nu_2, nu_3)
+    l_str_R = ("R1", "R2", "R3")
 
     u1 = rand()
     u2 = rand()
diff --git a/models/SIR.jl b/models/SIR.jl
index 5da6d2c1e53f399e5e32b0d654aea7a30b479e5b..82ea26f7bbb0a4f4323232e338c92737548ff4b1 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -13,13 +13,13 @@ function SIR_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{N
                 xn::Vector{Int}, tn::Float64, p::Vector{Float64})
     a1 = p[1] * xn[1] * xn[2]
     a2 = p[2] * xn[2]
-    l_a = SVector(a1, a2)
+    l_a = (a1, a2)
     asum = sum(l_a)
     # column-major order
-    nu_1 = SVector(-1, 1, 0)
-    nu_2 = SVector(0, -1, 1)
-    l_nu = SVector(nu_1, nu_2)
-    l_str_R = SVector("R1","R2")
+    nu_1 = (-1, 1, 0)
+    nu_2 = (0, -1, 1)
+    l_nu = (nu_1, nu_2)
+    l_str_R = ("R1","R2")
 
     u1 = rand()
     u2 = rand()
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index 251f35c9ec153b2da246e9f2bf04aa47befd9941..f12a93078d9ad28bbcdf356c9d4bff4a55ae69b0 100644
--- a/tests/run_unit.jl
+++ b/tests/run_unit.jl
@@ -26,6 +26,7 @@ using Test
     @test include("unit/load_module.jl")
     @test include("unit/long_sim_er.jl")
     
+    @test include("unit/macro_biochemical_network.jl")
     @test include("unit/model_prior.jl")
     @test include("unit/models_exps_er_1d.jl")
     @test include("unit/observe_all.jl")
diff --git a/tests/unit/macro_biochemical_network.jl b/tests/unit/macro_biochemical_network.jl
new file mode 100644
index 0000000000000000000000000000000000000000..5fab6c55eddbb7f9d7a8b1f3c0d7b6f76772dbe0
--- /dev/null
+++ b/tests/unit/macro_biochemical_network.jl
@@ -0,0 +1,37 @@
+
+using MarkovProcesses 
+
+model_SIR = @biochemical_network "SIR" begin
+    R1: (S+I => 2I, ki*S*I)
+    R2: (I => R, kr*I)
+end
+set_x0!(model_SIR, [95,5,0])
+set_param!(model_SIR, [0.012, 0.05])
+
+model_ER = @biochemical_network "ER" begin
+    R1: (E+S => ES, k1*E*S)
+    R2: (ES => E+S, k2*ES)
+    R3: (ES => E+P, k3*ES)
+end
+set_x0!(model_ER, [100,100,0,0])
+set_param!(model_ER, [1.0,1.0,1.0])
+
+return true
+
+#=
+@biochemical_network "test1" begin
+    R1: (S+I => 2I+Z,  ki*S*I)
+    R2: (I => R, kr*I)
+end
+
+@biochemical_network "test2" begin
+    R1: (S+I => 2I,  ki*S*I)
+    R2: (I => R, kr*I)
+end
+
+@biochemical_network "test3" begin
+    R1: (S+I => I,  ki)
+    R2: (2I => R, kr*I*c)
+end
+=#
+