diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml
index 4ba1f8e..1224ddc 100644
--- a/.github/workflows/python-package-conda.yml
+++ b/.github/workflows/python-package-conda.yml
@@ -12,19 +12,20 @@ jobs:
       max-parallel: 5
 
     steps:
-    - uses: actions/checkout@v2
-    - name: Set up Python 3.8
-      uses: actions/setup-python@v2
+    - uses: actions/checkout@v4
+    - name: Set up Python 3.11
+      uses: actions/setup-python@v5
       with:
-        python-version: 3.8
+        python-version: 3.11
     - name: Add conda to system path
       run: |
         # $CONDA is an environment variable pointing to the root of the miniconda directory
         echo $CONDA/bin >> $GITHUB_PATH
     - name: Install dependencies
       run: |
+        conda install -y python=3.11
         conda env update --file conda/environment.yml --name base
     - name: Test with pytest
       run: |
         conda install pytest pytest-cov
-        pytest --cov --cov-report term --cov-report xml --junitxml=xunit-result.xml --ignore='Testing/test_formal_methods.py'
\ No newline at end of file
+        pytest --cov --cov-report term --cov-report xml --junitxml=xunit-result.xml --ignore='Testing/test_formal_methods.py' --ignore='Testing/parsing/test_remote_parsing.py'
diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index 65025f3..09fedd1 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -5,12 +5,17 @@
 # Required
 version: 2
 
+# Set the version of Python and other tools you might need
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.8"
+
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
   builder: html
   configuration: docs/source/conf.py
 
 python:
-  version: "3.8"
   install:
     - requirements: docs/requirements.txt
diff --git a/Testing/models/get_model_str.py b/Testing/models/get_model_str.py
new file mode 100644
index 0000000..236c7cd
--- /dev/null
+++ b/Testing/models/get_model_str.py
@@ -0,0 +1,15 @@
+import os
+
+def get_model_str(model_name: str) -> str:
+    if not os.path.splitext(model_name)[1]:
+        model_name += ".txt"
+    script_dir = os.path.dirname(os.path.abspath(__file__))
+    file_path = os.path.join(script_dir, model_name)
+
+    try:
+        with open(file_path, "r") as file:
+            model_str = file.read()
+        return model_str
+    except FileNotFoundError:
+        print(f"File {file_path} not found")
+        return None
diff --git a/Testing/models/model1.txt b/Testing/models/model1.txt
new file mode 100644
index 0000000..8e9c574
--- /dev/null
+++ b/Testing/models/model1.txt
@@ -0,0 +1,12 @@
+#! rules
+X()::rep => @ k1*[X()::rep]
+Z()::rep => X()::rep
+=> Y()::rep @ 1/(1+([X()::rep])**4)
+
+#! inits
+2 X()::rep
+Y()::rep
+
+#! definitions
+k1 = 0.05
+k2 = 0.12 
diff --git a/Testing/models/model2.txt b/Testing/models/model2.txt
new file mode 100644
index 0000000..9b5d79b
--- /dev/null
+++ b/Testing/models/model2.txt
@@ -0,0 +1,12 @@
+#! rules
+X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep]
+X(T{a})::rep => X(T{o})::rep @ k2*[Z()::rep]
+=> Y(P{f})::rep @ 1/(1+([X()::rep])**4)
+
+#! inits
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+Y(P{g}, N{l})::rep
+
+#! definitions
+k1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/model3.txt b/Testing/models/model3.txt
new file mode 100644
index 0000000..1f9ec6d
--- /dev/null
+++ b/Testing/models/model3.txt
@@ -0,0 +1,12 @@
+#! rules
+Y{i}::rep => Y{a}::rep @ p*[Y{i}::rep]
+Y{i}::rep => Y{-}::rep @ (1-p)*[Y{i}::rep]
+X()::rep + Y{a}::rep => X().Y{a}::rep @ q*[X()::rep]*[Y{a}::rep]
+X(K{i}).Y{_}::rep => X(K{p}).Y{_}::rep @ p*[X(K{i}).Y{_}::rep] // also here
+
+#! inits
+2 X(K{i})::rep
+1 Y{i}::rep
+
+#! definitions
+p = 0.3
diff --git a/Testing/models/model4.txt b/Testing/models/model4.txt
new file mode 100644
index 0000000..9eadbe7
--- /dev/null
+++ b/Testing/models/model4.txt
@@ -0,0 +1,11 @@
+#! rules
+X()::rep => @ k1*[X()::rep]
+Z()::rep => X()::rep @ k2
+=> Y()::rep @ 1/(1+([X()::rep])**4)
+
+#! inits
+2 X()::rep
+Y()::rep
+
+#! definitions
+k2 = 5
diff --git a/Testing/models/model5.txt b/Testing/models/model5.txt
new file mode 100644
index 0000000..9d04053
--- /dev/null
+++ b/Testing/models/model5.txt
@@ -0,0 +1,12 @@
+#! rules
+// commenting
+X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep] // also here
+X(T{a})::rep => X(T{o})::rep @ k2*[Z()::rep]
+=> Y(P{f})::rep @ 1/(1+([X()::rep])**4) // ** means power (^)
+
+#! inits
+2 X(K{c}, T{e})::rep
+Y(P{g}, N{l})::rep
+
+#! definitions
+k2 = 0.05 // also comment
diff --git a/Testing/models/model6.txt b/Testing/models/model6.txt
new file mode 100644
index 0000000..98b3363
--- /dev/null
+++ b/Testing/models/model6.txt
@@ -0,0 +1,12 @@
+#! rules
+X()::rep => @ k1*[X()::rep]
+Z()::rep => X()::rep @ k2
+=> Y()::rep @ 1/(1+([X()::rep])**4)
+
+#! inits
+2 X()::rep
+Y()::rep
+
+#! definitions
+k2 = 5
+k1 = 2
diff --git a/Testing/models/model7.txt b/Testing/models/model7.txt
new file mode 100644
index 0000000..884c75d
--- /dev/null
+++ b/Testing/models/model7.txt
@@ -0,0 +1,13 @@
+#! rules
+// commenting
+X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep] // also here
+X(T{a})::rep => X(T{o})::rep @ k2*[Z()::rep]
+=> Y(P{f})::rep @ 1/(1+([X()::rep])**4) // ** means power (^)
+
+#! inits
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+Y(P{g}, N{l})::rep
+
+#! definitions
+k2 = 0.05 // also comment
+k1 = 2
diff --git a/Testing/models/model_TS.txt b/Testing/models/model_TS.txt
new file mode 100644
index 0000000..ca0b309
--- /dev/null
+++ b/Testing/models/model_TS.txt
@@ -0,0 +1,15 @@
+#! rules
+=> K(S{u},T{i})::cyt @ omega
+K(S{u})::cyt => K(S{p})::cyt @ alpha*[K(S{u})::cyt]
+K(S{p})::cyt + B{a}::cyt => K(S{p}).B{a}::cyt @ beta*[K(S{p})::cyt]*[B{a}::cyt]
+B{_}::cyt => @ gamma*[B{_}::cyt]
+K(S{u},T{i}).B{a}::cyt => @ 5
+
+#! inits
+1 B{a}::cyt
+
+#! definitions
+alpha = 10
+beta = 5
+gamma = 2
+omega = 3
diff --git a/Testing/models/model_abstract.txt b/Testing/models/model_abstract.txt
new file mode 100644
index 0000000..323643c
--- /dev/null
+++ b/Testing/models/model_abstract.txt
@@ -0,0 +1,10 @@
+#! rules
+=> X()::rep @ k2*[T{_}::rep]
+T{a}::rep => T{i}::rep @ k1*[T{_}::rep]
+
+#! inits
+10 T{a}::rep
+
+#! definitions
+k1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/model_bigger_TS.txt b/Testing/models/model_bigger_TS.txt
new file mode 100644
index 0000000..c6272e7
--- /dev/null
+++ b/Testing/models/model_bigger_TS.txt
@@ -0,0 +1,15 @@
+#! rules
+=> 2 K(S{u},T{i})::cyt @ omega
+K(S{u})::cyt => K(S{p})::cyt @ alpha*[K(S{u})::cyt]
+K(S{p})::cyt + B{a}::cyt => K(S{p}).B{a}::cyt @ beta*[K(S{p})::cyt]*[B{a}::cyt]
+B{_}::cyt => @ gamma*[B{_}::cyt]
+K(S{u},T{i}).B{a}::cyt => @ 5
+
+#! inits
+6 B{a}::cyt
+
+#! definitions
+alpha = 10
+beta = 5
+gamma = 2
+omega = 3
diff --git a/Testing/models/model_cmplx_in_abstr_seq1.txt b/Testing/models/model_cmplx_in_abstr_seq1.txt
new file mode 100644
index 0000000..35fc5d3
--- /dev/null
+++ b/Testing/models/model_cmplx_in_abstr_seq1.txt
@@ -0,0 +1,5 @@
+#! rules
+S{i}::A()::A2::cell  => A()::cell
+
+#! complexes
+A2 = A().A()
diff --git a/Testing/models/model_cmplx_in_abstr_seq2.txt b/Testing/models/model_cmplx_in_abstr_seq2.txt
new file mode 100644
index 0000000..a8aeb53
--- /dev/null
+++ b/Testing/models/model_cmplx_in_abstr_seq2.txt
@@ -0,0 +1,2 @@
+#! rules
+S{i}::A()::A().A()::cell => A()::cell
diff --git a/Testing/models/model_even_bigger_TS.txt b/Testing/models/model_even_bigger_TS.txt
new file mode 100644
index 0000000..9fb4726
--- /dev/null
+++ b/Testing/models/model_even_bigger_TS.txt
@@ -0,0 +1,15 @@
+#! rules
+=> K(S{u},T{i})::cyt @ omega
+K(S{u})::cyt => K(S{p})::cyt @ alpha*[K(S{u})::cyt]
+K(S{p})::cyt + B{a}::cyt => K(S{p}).B{a}::cyt @ beta*[K(S{p})::cyt]*[B{a}::cyt]
+B{_}::cyt => @ gamma*[B{_}::cyt]
+K(S{u},T{i}).B{a}::cyt => @ 5
+
+#! inits
+10 B{a}::cyt
+
+#! definitions
+alpha = 10
+beta = 5
+gamma = 2
+omega = 3
diff --git a/Testing/models/model_for_bound.txt b/Testing/models/model_for_bound.txt
new file mode 100644
index 0000000..4af40cc
--- /dev/null
+++ b/Testing/models/model_for_bound.txt
@@ -0,0 +1,9 @@
+#! rules
+=> A{i}::cyt @ k
+=> C{i}::cyt @ k
+
+#! inits
+2 B{i}::cyt
+
+#! definitions
+k = 1
diff --git a/Testing/models/model_for_matching.txt b/Testing/models/model_for_matching.txt
new file mode 100644
index 0000000..7d36857
--- /dev/null
+++ b/Testing/models/model_for_matching.txt
@@ -0,0 +1,16 @@
+#! rules
+K(S{i}).B()::cyt + C{_}::cell => K(S{i})::cyt + B()::cyt + C{_}::cell @ 3*[K(S{i}).B()::cyt]/2*v_1
+A{p}.K(S{i})::cyt => A{i}::cyt + K(S{a})::cyt + C{a}::cyt @ 0.3*[A{p}.K(S{i})::cyt]
+K(S{i},T{i})::cyt + C{_}::cell => K(S{a},T{i})::cyt @ k2*[K(S{i},T{i})::cyt]
+C{_}::cell + C{_}::cell => C{_}.C{_}::cell @ v_1*[C{_}::cell]**2
+C{_}::cell + K()::cell => C{_}.K()::cell @ v_1*[C{_}::cell]**2
+
+#! inits
+2 K(S{i},T{i}).B(T{a})::cyt
+1 A{p}.K(S{i},T{i})::cyt
+2 C{i}::cell
+1 C{a}::cell
+
+#! definitions
+v_1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/model_miyoshi.txt b/Testing/models/model_miyoshi.txt
new file mode 100644
index 0000000..ce014a7
--- /dev/null
+++ b/Testing/models/model_miyoshi.txt
@@ -0,0 +1,28 @@
+#! rules
+S{u}::KaiC()::KaiC6::cyt => S{p}::KaiC()::KaiC6::cyt @ (kcat1*[KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
+S{p}::KaiC()::KaiC6::cyt => S{u}::KaiC()::KaiC6::cyt @ (kcat2*[KaiB4{a}.KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
+T{u}::KaiC()::KaiC6::cyt => T{p}::KaiC()::KaiC6::cyt @ (kcat3*[KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
+T{p}::KaiC()::KaiC6::cyt => T{u}::KaiC()::KaiC6::cyt @ (kcat4*[KaiB4{a}.KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
+KaiB4{i}::cyt => KaiB4{a}::cyt @ (kcatb2*[KaiB4{i}::cyt])/(Kmb2 + [KaiB4{i}::cyt])
+KaiB4{a}::cyt => KaiB4{i}::cyt @ (kcatb1*[KaiB4{a}::cyt])/(Kmb1 + [KaiB4{a}::cyt])
+KaiB4{a}.KaiA2()::cyt => KaiB4{a}::cyt + KaiA2()::cyt @ k12*[KaiB4{a}.KaiA2()::cyt]
+KaiC6::cyt => 6 KaiC()::cyt @ kdimer*[KaiC6::cyt]
+6 KaiC()::cyt => KaiC6::cyt @ kdimer*[KaiC()::cyt]*([KaiC()::cyt] - 1)*([KaiC()::cyt] - 2)*([KaiC()::cyt] - 3)*([KaiC()::cyt] - 4)*([KaiC()::cyt] - 5)
+
+#! inits
+6 KaiC(S{p},T{p})::cyt
+1 KaiB4{a}.KaiA2()::cyt
+
+#! definitions
+kcat1 = 0.539
+kcat3 = 0.89
+Km = 0.602
+kcatb2 = 0.346
+kcatb1 = 0.602
+Kmb2 = 66.75
+Kmb1 = 2.423
+k12 = 0.0008756
+kdimer = 1.77
+
+#! complexes
+KaiC6 = KaiC().KaiC().KaiC().KaiC().KaiC().KaiC()
diff --git a/Testing/models/model_miyoshi_non_param.txt b/Testing/models/model_miyoshi_non_param.txt
new file mode 100644
index 0000000..d47c22c
--- /dev/null
+++ b/Testing/models/model_miyoshi_non_param.txt
@@ -0,0 +1,33 @@
+#! rules
+S{u}:KaiC():KaiC6::cyt => S{p}:KaiC():KaiC6::cyt @ (kcat1*[KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
+S{p}:KaiC():KaiC6::cyt => S{u}:KaiC():KaiC6::cyt @ (kcat2*[KaiB4{a}.KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
+T{u}:KaiC():KaiC6::cyt => T{p}:KaiC():KaiC6::cyt @ (kcat3*[KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
+T{p}:KaiC():KaiC6::cyt => T{u}:KaiC():KaiC6::cyt @ (kcat4*[KaiB4{a}.KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
+KaiB4{i}::cyt => KaiB4{a}::cyt @ (kcatb2*[KaiB4{i}::cyt])/(Kmb2 + [KaiB4{i}::cyt])
+KaiB4{a}::cyt => KaiB4{i}::cyt @ (kcatb1*[KaiB4{a}::cyt])/(Kmb1 + [KaiB4{a}::cyt])
+KaiB4{a}.KaiA2()::cyt => KaiB4{a}::cyt + KaiA2()::cyt @ k12*[KaiB4{a}.KaiA2()::cyt]
+KaiC6::cyt => 6 KaiC()::cyt @ kdimer*[KaiC6::cyt]
+6 KaiC()::cyt => KaiC6::cyt @ kdimer*[KaiC()::cyt]
+KaiC().KaiC()::cyt => KaiC()::cyt + KaiC()::cyt @ kcat1
+
+#! inits
+6 KaiC(S{p},T{p})::cyt
+1 KaiB4{a}.KaiA2()::cyt
+1 KaiC(S{u},T{p})::cyt
+1 KaiC(S{u},T{p}).KaiC(S{p},T{p})::cyt
+
+#! definitions
+kcat2 = 0.539
+kcat4 = 0.89
+kcat1 = 0.539
+kcat3 = 0.89
+Km = 0.602
+kcatb2 = 0.346
+kcatb1 = 0.602
+Kmb2 = 66.75
+Kmb1 = 2.423
+k12 = 0.0008756
+kdimer = 1.77
+
+#! complexes
+KaiC6 = KaiC().KaiC().KaiC().KaiC().KaiC().KaiC()
diff --git a/Testing/models/model_nonreachable.txt b/Testing/models/model_nonreachable.txt
new file mode 100644
index 0000000..6552fad
--- /dev/null
+++ b/Testing/models/model_nonreachable.txt
@@ -0,0 +1,11 @@
+#! rules
+K(S{i}).B()::cyt => K(S{a})::cyt + B()::cyt @ 3*[K(S{i}).B()::cyt]/2*v_1
+K(S{a})::cyt + A{i}::cyt => K(S{a}).A{i}::cyt
+
+#! inits
+2 K(S{i}).B()::cyt
+1 A{i}::cyt
+
+#! definitions
+v_1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/model_parametrised.txt b/Testing/models/model_parametrised.txt
new file mode 100644
index 0000000..5daea04
--- /dev/null
+++ b/Testing/models/model_parametrised.txt
@@ -0,0 +1,12 @@
+#! rules
+// commenting
+X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep] // also here
+X(T{a})::rep => X(T{o})::rep @ k2*[Z()::rep]
+=> Y(P{f})::rep @ 1/(v_3+([X()::rep])**4) // ** means power (^)
+
+#! inits
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+Y(P{g}, N{l})::rep // comment just 1 item
+
+#! definitions
+k1 = 0.05
diff --git a/Testing/models/model_parametrised2.txt b/Testing/models/model_parametrised2.txt
new file mode 100644
index 0000000..8a9d81a
--- /dev/null
+++ b/Testing/models/model_parametrised2.txt
@@ -0,0 +1,15 @@
+#! rules
+=> K(S{u},T{i})::cyt @ omega
+K(S{u})::cyt => K(S{p})::cyt @ alpha*[K(S{u})::cyt]
+K(S{p})::cyt + B{a}::cyt => K(S{p}).B{a}::cyt @ beta*[K(S{p})::cyt]*[B{a}::cyt]
+B{_}::cyt => @ gamma*[B{_}::cyt]
+K(S{u},T{i}).B{a}::cyt => @ 5
+
+#! inits
+1 B{a}::cyt
+
+#! definitions
+alpha = 10
+beta = 5
+//gamma = 2
+omega = 3
diff --git a/Testing/models/model_reachable.txt b/Testing/models/model_reachable.txt
new file mode 100644
index 0000000..fedfbb9
--- /dev/null
+++ b/Testing/models/model_reachable.txt
@@ -0,0 +1,12 @@
+#! rules
+K(S{i}).B()::cyt => K(S{a})::cyt + B()::cyt @ 3*[K(S{i}).B()::cyt]/2*v_1
+K(S{a})::cyt + A{i}::cyt => K(S{a}).A{i}::cyt
+K().A{i}::cyt => K().A{a}::cyt
+
+#! inits
+2 K(S{i}).B()::cyt
+1 A{i}::cyt
+
+#! definitions
+v_1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/model_with_comments.txt b/Testing/models/model_with_comments.txt
new file mode 100644
index 0000000..2022186
--- /dev/null
+++ b/Testing/models/model_with_comments.txt
@@ -0,0 +1,15 @@
+#! rules
+// commenting
+X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep] // also here
+X(T{a})::rep => X(T{o})::rep @ k2*[Z()::rep]
+=> Y(P{f})::rep @ 1/(1+([X()::rep])**4) // ** means power (^)
+
+#! inits
+// here
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+Y(P{g}, N{l})::rep // comment just 1 item
+
+#! definitions
+// and
+k1 = 0.05 // also
+k2 = 0.12
diff --git a/Testing/models/model_with_complexes.txt b/Testing/models/model_with_complexes.txt
new file mode 100644
index 0000000..ed08060
--- /dev/null
+++ b/Testing/models/model_with_complexes.txt
@@ -0,0 +1,20 @@
+#! rules
+// commenting
+X(T{a})::XX::rep => X(T{o})::XX::rep @ k2*[X().X()::rep]
+K{i}::X()::XYZ::rep => K{p}::X()::XYZ::rep @ k1*[X().Y().Z()::rep] // also here
+=> P{f}::XP::rep @ 1/(1+([X().P{_}::rep])**4) // ** means power (^)
+
+#! inits
+// here
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+Y(P{g}, N{l})::rep // comment just 1 item
+
+#! definitions
+// and
+k1 = 0.05 // also
+k2 = 0.12
+
+#! complexes
+XYZ = X().Y().Z() // a big complex
+XX = X().X()
+XP = X().P{_}
diff --git a/Testing/models/model_with_context.txt b/Testing/models/model_with_context.txt
new file mode 100644
index 0000000..80da740
--- /dev/null
+++ b/Testing/models/model_with_context.txt
@@ -0,0 +1,12 @@
+#! rules
+K(S{i}).B(T{a})::cyt => K(S{i})::cyt + B(T{a})::cyt @ 3*[K(S{i}).B(T{a})::cyt]/2*v_1
+A{p}.K(S{i},T{i})::cyt => A{i}::cyt + K(S{a},T{a})::cyt
+K(S{i},T{i})::cyt => K(S{a},T{i})::cyt
+
+#! inits
+2 K(S{i}).B(T{a})::cyt
+1 A{p}.K(S{i},T{i})::cyt
+
+#! definitions
+v_1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/model_with_labels.txt b/Testing/models/model_with_labels.txt
new file mode 100644
index 0000000..dc67192
--- /dev/null
+++ b/Testing/models/model_with_labels.txt
@@ -0,0 +1,12 @@
+#! rules
+r1_S ~ A(S{i})::cell => A(S{a})::cell @ k1*[A(S{i})::cell]
+r1_T ~ A(T{i})::cell => A(T{a})::cell @ k2*[A(T{i})::cell]
+r2 ~ A()::cell => A()::out @ k3*[A()::cell]
+
+#! inits
+1 A(S{i},T{i})::cell
+
+#! definitions
+k1 = 0.3
+k2 = 0.5
+k3 = 0.1
diff --git a/Testing/models/model_with_redundant.txt b/Testing/models/model_with_redundant.txt
new file mode 100644
index 0000000..b0beb8a
--- /dev/null
+++ b/Testing/models/model_with_redundant.txt
@@ -0,0 +1,13 @@
+#! rules
+K(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt + D(A{_})::cell @ 3*[K().B()::cyt]/2*v_1
+K().B()::cyt => K()::cyt + B()::cyt + D(A{_})::cell @ 3*[K().B()::cyt]/2*v_1
+K().K()::cyt => K()::cyt + K()::cyt
+K(S{i}).K()::cyt => K(S{a})::cyt + K()::cyt
+K(S{i}, T{p}).K()::cyt => K(S{a}, T{p})::cyt + K()::cyt
+
+#! inits
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+
+#! definitions
+v_1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/model_with_sinks.txt b/Testing/models/model_with_sinks.txt
new file mode 100644
index 0000000..986dc5b
--- /dev/null
+++ b/Testing/models/model_with_sinks.txt
@@ -0,0 +1,12 @@
+#! rules
+K(S{u})::cyt => K(S{p})::cyt @ alpha*[K(S{u})::cyt]
+K(S{p})::cyt + B{a}::cyt => K(S{p}).B{a}::cyt @ beta*[K(S{p})::cyt]*[B{a}::cyt]
+B{a}::cyt => B{i}::cyt @ alpha*[B{_}::cyt]
+
+#! inits
+1 B{a}::cyt
+1 K(S{u})::cyt
+
+#! definitions
+alpha = 10
+beta = 5
diff --git a/Testing/models/model_with_variable.txt b/Testing/models/model_with_variable.txt
new file mode 100644
index 0000000..04fa9a0
--- /dev/null
+++ b/Testing/models/model_with_variable.txt
@@ -0,0 +1,17 @@
+#! rules
+// commenting
+T{a}::X()::?::rep => T{o}::X()::?::rep @ k2*[X().X()::rep] ; ? = { XX, XY }
+K{i}::X()::XY::rep => K{p}::X()::XY::rep @ k1*[X().Y().Z().X()::rep] // also here
+
+#! inits
+// here
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+
+#! definitions
+// and
+k1 = 0.05 // also
+k2 = 0.12
+
+#! complexes
+XX = X().X()
+XY = X().Y()
diff --git a/Testing/models/model_without_complexes.txt b/Testing/models/model_without_complexes.txt
new file mode 100644
index 0000000..4fd3420
--- /dev/null
+++ b/Testing/models/model_without_complexes.txt
@@ -0,0 +1,15 @@
+#! rules
+// commenting
+X(T{a}).X()::rep => X(T{o}).X()::rep @ k2*[X().X()::rep]
+X(K{i}).Y().Z()::rep => X(K{p}).Y().Z()::rep @ k1*[X().Y().Z()::rep] // also here
+=> X().P{f}::rep @ 1/(1+([X().P{_}::rep])**4) // ** means power (^)
+
+#! inits
+// here
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+Y(P{g}, N{l})::rep // comment just 1 item
+
+#! definitions
+// and
+k1 = 0.05 // also
+k2 = 0.12
diff --git a/Testing/models/model_without_context.txt b/Testing/models/model_without_context.txt
new file mode 100644
index 0000000..3e9fede
--- /dev/null
+++ b/Testing/models/model_without_context.txt
@@ -0,0 +1,11 @@
+#! rules
+K().B()::cyt => K()::cyt + B()::cyt @ 3*[K().B()::cyt]/2*v_1
+A{_}.K()::cyt => A{_}::cyt + K()::cyt
+
+#! inits
+2 K().B()::cyt
+1 A{_}.K()::cyt
+
+#! definitions
+v_1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/model_without_redundant.txt b/Testing/models/model_without_redundant.txt
new file mode 100644
index 0000000..c23f6e8
--- /dev/null
+++ b/Testing/models/model_without_redundant.txt
@@ -0,0 +1,10 @@
+#! rules
+K().B()::cyt => K()::cyt + B()::cyt + D(A{_})::cell @ 3*[K().B()::cyt]/2*v_1
+K().K()::cyt => K()::cyt + K()::cyt
+
+#! inits
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+
+#! definitions
+v_1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/model_without_variable.txt b/Testing/models/model_without_variable.txt
new file mode 100644
index 0000000..e5cbe46
--- /dev/null
+++ b/Testing/models/model_without_variable.txt
@@ -0,0 +1,14 @@
+#! rules
+// commenting
+X(K{i}).Y()::rep => X(K{p}).Y()::rep @ k1*[X().Y().Z().X()::rep]
+X(T{a}).X()::rep => X(T{o}).X()::rep @ k2*[X().X()::rep]
+X(T{a}).Y()::rep => X(T{o}).Y()::rep @ k2*[X().X()::rep]
+
+#! inits
+// here
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+
+#! definitions
+// and
+k1 = 0.05 // also
+k2 = 0.12
diff --git a/Testing/models/model_wrong1.txt b/Testing/models/model_wrong1.txt
new file mode 100644
index 0000000..876a931
--- /dev/null
+++ b/Testing/models/model_wrong1.txt
@@ -0,0 +1,12 @@
+#! rules
+X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep]
+X(T{a})::rep => X(T{o}):;rep @ k2*[Z()::rep]
+=> Y(P{f})::rep @ 1/(1+([X()::rep])**4)
+
+#! inits
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+Y(P{g}, N{l})::rep
+
+#! definitions
+k1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/model_wrong2.txt b/Testing/models/model_wrong2.txt
new file mode 100644
index 0000000..1c73561
--- /dev/null
+++ b/Testing/models/model_wrong2.txt
@@ -0,0 +1,12 @@
+#! rules
+X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep]
+X(T{a})::rep = X(T{o})::rep @ k2*[Z()::rep]
+=> Y(P{f})::rep @ 1/(1+([X()::rep])**4)
+
+#! inits
+2 X(K{c}, T{e}).X(K{c}, T{j})::rep
+Y(P{g}, N{l})::rep
+
+#! definitions
+k1 = 0.05
+k2 = 0.12
diff --git a/Testing/models/regulation1.txt b/Testing/models/regulation1.txt
new file mode 100644
index 0000000..7edbece
--- /dev/null
+++ b/Testing/models/regulation1.txt
@@ -0,0 +1,4 @@
+#! regulation
+type programmed
+r1_S: {r1_T, r2}
+r1_T: {r1_S}
diff --git a/Testing/models/regulation2.txt b/Testing/models/regulation2.txt
new file mode 100644
index 0000000..3ed8fbc
--- /dev/null
+++ b/Testing/models/regulation2.txt
@@ -0,0 +1,3 @@
+#! regulation
+type ordered
+(r1_S, r2), (r1_T, r2)
diff --git a/Testing/models/regulation3.txt b/Testing/models/regulation3.txt
new file mode 100644
index 0000000..5339364
--- /dev/null
+++ b/Testing/models/regulation3.txt
@@ -0,0 +1,3 @@
+#! regulation
+type conditional
+r2: {A(S{a},T{i})::cell}
diff --git a/Testing/models/regulation4.txt b/Testing/models/regulation4.txt
new file mode 100644
index 0000000..81e1a1e
--- /dev/null
+++ b/Testing/models/regulation4.txt
@@ -0,0 +1,3 @@
+#! regulation
+type concurrent-free
+(r1_S, r2), (r1_T, r2)
diff --git a/Testing/models/regulation5.txt b/Testing/models/regulation5.txt
new file mode 100644
index 0000000..ae05ba3
--- /dev/null
+++ b/Testing/models/regulation5.txt
@@ -0,0 +1,3 @@
+#! regulation
+type regular
+(r1_S;r1_T;r2|r1_T;r1_S;r2)
diff --git a/Testing/models/regulation6.txt b/Testing/models/regulation6.txt
new file mode 100644
index 0000000..7edbece
--- /dev/null
+++ b/Testing/models/regulation6.txt
@@ -0,0 +1,4 @@
+#! regulation
+type programmed
+r1_S: {r1_T, r2}
+r1_T: {r1_S}
diff --git a/Testing/models/tumor_model.txt b/Testing/models/tumor_model.txt
new file mode 100644
index 0000000..2849009
--- /dev/null
+++ b/Testing/models/tumor_model.txt
@@ -0,0 +1,15 @@
+#! rules
+T(P{i})::x => T(P{m})::x @ a1*[T(P{i})::x]
+T(P{m})::x => T(P{i})::x + T(P{i})::x @ a2*[T(P{m})::x]
+T(P{i})::x => @ d1*[T(P{i})::x]
+T(P{m})::x => @ d2*[T(P{m})::x]
+
+#! inits
+2 T(P{m})::x
+1 T(P{i})::x
+
+#! definitions
+a1 = 1.2
+a2 = 2
+d1 = 0.8
+d2 = 0.5
diff --git a/Testing/models/tumor_parametric_model.txt b/Testing/models/tumor_parametric_model.txt
new file mode 100644
index 0000000..60c2fc1
--- /dev/null
+++ b/Testing/models/tumor_parametric_model.txt
@@ -0,0 +1,14 @@
+#! rules
+T(P{i})::x => T(P{m})::x @ a1*[T(P{i})::x]
+T(P{m})::x => T(P{i})::x + T(P{i})::x @ a2*[T(P{m})::x]
+T(P{i})::x => @ d1*[T(P{i})::x]
+T(P{m})::x => @ d2*[T(P{m})::x]
+
+#! inits
+2 T(P{m})::x
+1 T(P{i})::x
+
+#! definitions
+a1 = 1.2
+a2 = 2
+d1 = 0.8
diff --git a/Testing/objects_testing.py b/Testing/objects_testing.py
new file mode 100644
index 0000000..171ba58
--- /dev/null
+++ b/Testing/objects_testing.py
@@ -0,0 +1,440 @@
+import numpy as np
+import collections
+
+from eBCSgen.Core.Atomic import AtomicAgent
+from eBCSgen.Core.Rule import Rule
+from eBCSgen.Core.Structure import StructureAgent
+from eBCSgen.Core.Complex import Complex
+from eBCSgen.Core.Side import Side
+from eBCSgen.Core.Rate import Rate
+from eBCSgen.TS.State import Vector, State, Memory
+from eBCSgen.Core.Reaction import Reaction
+
+from eBCSgen.Parsing.ParseBCSL import Parser
+
+rate_parser = Parser("rate")
+atomic_parser = Parser("atomic")
+structure_parser = Parser("structure")
+complex_parser = Parser("complex")
+state_parser = Parser("state")
+side_parser = Parser("side")
+rate_complex_parser = Parser("rate_complex")
+rule_parser = Parser("rule")
+rules_parser = Parser("rules")
+model_parser = Parser("model")
+observables_parser = Parser("observables")
+observable_parser = Parser("observable")
+
+# atomic
+a1 = AtomicAgent("T", "s")
+a2 = AtomicAgent("S", "a")
+a3 = AtomicAgent("S", "s")
+a4 = AtomicAgent("T", "_")
+a5 = AtomicAgent("S", "_")
+a6 = AtomicAgent("T", "p")
+a7 = AtomicAgent("T", "u")
+
+a8 = AtomicAgent("T", "s")
+a9 = AtomicAgent("T", "i")
+a10 = AtomicAgent("S", "i")
+a11 = AtomicAgent("U", "a")
+a12 = AtomicAgent("U", "_")
+a13 = AtomicAgent("U", "b")
+a14 = AtomicAgent("T", "a")
+
+a15 = AtomicAgent("S", "u")
+a16 = AtomicAgent("S", "p")
+a17 = AtomicAgent("B", "_")
+a18 = AtomicAgent("B", "-")
+a19 = AtomicAgent("B", "+")
+
+a20 = AtomicAgent("B", "a")
+
+a4_p = AtomicAgent("C", "p")
+a4_u = AtomicAgent("C", "u")
+
+u2_c1_p = AtomicAgent("U", "p")
+u2_c1_u = AtomicAgent("U", "u")
+
+
+# structure
+s1 = StructureAgent("B", {a1})
+s2 = StructureAgent("D", set())
+s3 = StructureAgent("K", {a1, a3, a11})
+s4 = StructureAgent("B", {a4})
+s5 = StructureAgent("D", {a5, a6})
+s6 = StructureAgent("K", set())
+s7 = StructureAgent("K", {a2, a4})
+
+s8 = StructureAgent("X", {a1})
+s9 = StructureAgent("A", {a10, a11})
+s10 = StructureAgent("X", {a4})
+s11 = StructureAgent("A", {a10, a12})
+s12 = StructureAgent("A", {a5, a11})
+s13 = StructureAgent("A", set())
+s14 = StructureAgent("A", {a10, a13})
+
+s15 = StructureAgent("strA", {a1, a2})
+s16 = StructureAgent("strA", {a2, a4})
+s18 = StructureAgent("strD", set())
+s19 = StructureAgent("strA", {a2, a1})
+s20 = StructureAgent("strA", {a1})
+s21 = StructureAgent("strA", {a6})
+s23 = StructureAgent("strA", {a6, a2})
+s24 = StructureAgent("strA", {a1, a10})
+s25 = StructureAgent("strA", {a6, a10})
+s26 = StructureAgent("strA", set())
+
+s27 = StructureAgent("X", {a2})
+s28 = StructureAgent("X", {a10})
+s29 = StructureAgent("K", {a14})
+s30 = StructureAgent("K", {a9})
+s31 = StructureAgent("X", set())
+
+s32 = StructureAgent("S", set())
+
+s34 = StructureAgent("K", {a15})
+s35 = StructureAgent("B", set())
+s36 = StructureAgent("K", {a16})
+s37 = StructureAgent("B", set())
+s38 = StructureAgent("D", {a17})
+s39 = StructureAgent("K", {a18})
+s40 = StructureAgent("K", {a19})
+
+s41 = StructureAgent("Y", set())
+s42 = StructureAgent("Z", set())
+s43 = StructureAgent("W", set())
+s44 = StructureAgent("D", {a4})
+s45 = StructureAgent("B", {a14, a16})
+
+s46 = StructureAgent("K", {a16, a9})
+s47 = StructureAgent("K", {a15, a9})
+
+s6_c1_p = StructureAgent("D", {a4_p})
+s6_c1_u = StructureAgent("D", {a4_u})
+
+s2_c1_p = StructureAgent("B", {u2_c1_p})
+s2_c1_u = StructureAgent("B", {u2_c1_u})
+
+s1_c1_a = StructureAgent("K", {a15, a14})
+
+s3_c1_a = StructureAgent("K", {a16, a14})
+
+# complex
+c1 = Complex([s1], "cell")
+c2 = Complex([s1, s2, s3], "cell")
+c3 = Complex([s1, s3, s5], "cyt")
+c4 = Complex([s1, a1], "cyt")
+c5 = Complex([s2] * 3 + [s4] * 3, "cell")
+c6 = Complex([s6] * 2 + [a6], "cyt")
+c7 = Complex([a7, a4], "cell")
+
+c8 = Complex([s8, s9, s9], "cyt")
+c9 = Complex([s10, s11, s12], "cyt")
+c10 = Complex([s9, s9, s8], "cyt")
+c11 = Complex([s9, s9, s8], "cell")
+
+c12 = Complex([s9, s11, a1], "cell")
+c13 = Complex([s13, s13, a4], "cell")
+
+large_c1 = Complex([s11] * 6 + [s10] * 5, "cell")
+large_c2 = Complex([s12] * 7 + [s10] * 6, "cell")
+
+c14 = Complex([s6], "cyt")  # K()::cyt
+c15 = Complex([s29], "cyt")  # K(T{a})::cyt
+c16 = Complex([s30], "cyt")  # K(T{i})::cyt
+c17 = Complex([s30, s27], "cyt")  # K(T{i}).X(S{a})::cyt
+c18 = Complex([s30, s28], "cyt")  # K(T{i}).X(S{i})::cyt
+c19 = Complex([s29, s27], "cyt")  # K(T{a}).X(S{a})::cyt
+c20 = Complex([s29, s28], "cyt")  # K(T{a}).X(S{i})::cyt
+
+c21 = Complex([s32], "rep")
+c22 = Complex([s2], "rep")
+
+c23 = Complex([s34, s35], "cyt")
+c24 = Complex([s36], "cyt")
+c25 = Complex([s35], "cyt")
+c26 = Complex([s38], "cell")
+
+c27 = Complex([s31], "rep")
+c28 = Complex([s41], "rep")
+c29 = Complex([s42], "rep")
+c30 = Complex([s43], "rep")
+
+c1_c1 = Complex([s2_c1_u], "cyt")  # B(U{u})::cyt
+c1_c2 = Complex([s2_c1_p], "cyt")  # B(U{p})::cyt
+c1_c3 = Complex([s1_c1_a], "cyt")  # K(S{u},T{a})::cyt
+c1_c5 = Complex([s3_c1_a, s2_c1_u], "cyt")  # K(S{p},T{a}).B(U{u})::c
+c1_c6 = Complex([s46, s2_c1_u], "cyt")  # K(S{p},T{i}).B(U{u})::c
+c1_c7 = Complex([s46, s2_c1_p], "cyt")  # K(S{p},T{i}).B(U{p})::c
+c1_c8 = Complex([s3_c1_a, s2_c1_p], "cyt")  # K(S{p},T{a}).B(U{p})::c
+c1_c9 = Complex([s6_c1_p], "cyt")  # D(C{p})::cyt
+c1_c10 = Complex([s6_c1_u], "cyt")  # D(C{u})::cyt
+sequence_no_change = (s1_c1_a, s2_c1_u, s3_c1_a, s2_c1_u, s6_c1_p)
+
+counter_c1 = Complex(collections.Counter({s35: 2}), "cyt")
+counter_c2 = Complex(collections.Counter({s36: 1}), "cyt")
+counter_c3 = Complex(collections.Counter({s35: 1}), "cyt")
+counter_c4 = Complex(collections.Counter({s44: 1}), "cell")
+counter_c5 = Complex(collections.Counter({s37: 1}), "cell")
+counter_c6 = Complex(collections.Counter({s45: 1}), "cell")
+
+cx1 = Complex([a20], "cyt")
+cx2 = Complex([s46], "cyt")
+cx3 = Complex([s47], "cyt")  # K(S{u},T{i})::cyt
+cx4 = Complex([s46, a20], "cyt")
+cx5 = Complex([s47, a20], "cyt")
+
+# side
+side1 = Side([c1])
+side2 = Side([c4, c5])
+side3 = Side([c1, c2, c2])
+side4 = Side([c7, c6, c5])
+side5 = Side([])
+
+lhs = Side([c23])
+rhs = Side([c24, c25, c26])
+
+side6 = Side([counter_c1, counter_c2, counter_c4])
+side7 = Side([counter_c2, counter_c3, counter_c4])
+side8 = Side([counter_c3, counter_c4, counter_c2])
+side9 = Side([counter_c6, counter_c1])
+side10 = Side([counter_c5, counter_c1])
+side11 = Side([counter_c5, counter_c1, counter_c2])
+side12 = Side([counter_c6, counter_c1, counter_c3, counter_c4])
+
+# states
+state1 = State(Vector(np.array([2, 3])), Memory(0))
+state2 = State(Vector(np.array([2, 0, 3, 1, 6, 2])), Memory(0))
+state3 = State(Vector(np.array((0, 0))), Memory(0))
+state4 = State(Vector(np.array((1, 0))), Memory(0))
+state5 = State(Vector(np.array((2, 0))), Memory(0))
+state6 = State(Vector(np.array((3, 0))), Memory(0))
+state7 = State(Vector(np.array((4, 0))), Memory(0))
+state8 = State(Vector(np.array((5, 0))), Memory(0))
+state9 = State(Vector(np.array((6, 0))), Memory(0))
+state10 = State(Vector(np.array((7, 1))), Memory(0))
+state11 = State(Vector(np.array((7, 2))), Memory(0))
+state12 = State(Vector(np.array((7, 3))), Memory(0))
+state13 = State(Vector(np.array((7, 4))), Memory(0))
+state14 = State(Vector(np.array((7, 5))), Memory(0))
+state15 = State(Vector(np.array((7, 6))), Memory(0))
+
+# rules
+
+sequence_1 = (s31,)
+mid_1 = 1
+compartments_1 = ["rep"]
+complexes_1 = [(0, 0)]
+pairs_1 = [(0, None)]
+rate_1 = Rate("k1*[X()::rep]")
+
+r1 = Rule(sequence_1, mid_1, compartments_1, complexes_1, pairs_1, rate_1)
+
+sequence_2 = (s42, s31)
+mid_2 = 1
+compartments_2 = ["rep"] * 2
+complexes_2 = [(0, 0), (1, 1)]
+pairs_2 = [(0, 1)]
+
+r2 = Rule(sequence_2, mid_2, compartments_2, complexes_2, pairs_2, None)
+
+sequence_3 = (s41,)
+mid_3 = 0
+compartments_3 = ["rep"]
+complexes_3 = [(0, 0)]
+pairs_3 = [(None, 0)]
+rate_3 = Rate("1.0/(1.0+([X()::rep])**4.0)")
+
+r3 = Rule(sequence_3, mid_3, compartments_3, complexes_3, pairs_3, rate_3)
+
+sequence_4 = (s34, s35, s36, s37)
+reversed_sequence_4 = (s36, s37, s34, s35)
+mid_4 = 2
+compartments_4 = ["cyt"] * 4
+complexes_4 = [(0, 1), (2, 2), (3, 3)]
+reversed_complexes_4 = [(0, 0), (1, 1), (2, 3)]
+pairs_4 = [(0, 2), (1, 3)]
+rate_4 = Rate("3.0*[K()::cyt]/2.0*v_1")
+reversed_rate_4 = Rate("2.0*[K()::cyt]/3.0*v_1")
+
+r4 = Rule(sequence_4, mid_4, compartments_4, complexes_4, pairs_4, rate_4)
+reversed_r4a = Rule(
+    reversed_sequence_4, mid_4, compartments_4, reversed_complexes_4, pairs_4, rate_4
+)
+reversed_r4b = Rule(
+    reversed_sequence_4,
+    mid_4,
+    compartments_4,
+    reversed_complexes_4,
+    pairs_4,
+    reversed_rate_4,
+)
+sequence_one_side_bidirectional = (s36, s37)
+mid_one_side_bidirectional_a = 2
+mid_one_side_bidirectional_b = 0
+compartments_one_side_bidirectional = ["cyt"] * 2
+complexes_one_side_bidirectional = [(0, 0), (1, 1)]
+pairs_one_side_bidirectional_a = [(0, None), (1, None)]
+pairs_one_side_bidirectional_b = [(None, 0), (None, 1)]
+one_side_bidirectional_a = Rule(
+    sequence_one_side_bidirectional,
+    mid_one_side_bidirectional_a,
+    compartments_one_side_bidirectional,
+    complexes_one_side_bidirectional,
+    pairs_one_side_bidirectional_a,
+    rate_4,
+)
+one_side_bidirectional_b = Rule(
+    sequence_one_side_bidirectional,
+    mid_one_side_bidirectional_b,
+    compartments_one_side_bidirectional,
+    complexes_one_side_bidirectional,
+    pairs_one_side_bidirectional_b,
+    rate_4,
+)
+one_side_bidirectional_b_reversed_rate = Rule(
+    sequence_one_side_bidirectional,
+    mid_one_side_bidirectional_b,
+    compartments_one_side_bidirectional,
+    complexes_one_side_bidirectional,
+    pairs_one_side_bidirectional_b,
+    reversed_rate_4,
+)
+one_side_bidirectional_a_no_rate = Rule(
+    sequence_one_side_bidirectional,
+    mid_one_side_bidirectional_a,
+    compartments_one_side_bidirectional,
+    complexes_one_side_bidirectional,
+    pairs_one_side_bidirectional_a,
+    None,
+)
+one_side_bidirectional_b_no_rate = Rule(
+    sequence_one_side_bidirectional,
+    mid_one_side_bidirectional_b,
+    compartments_one_side_bidirectional,
+    complexes_one_side_bidirectional,
+    pairs_one_side_bidirectional_b,
+    None,
+)
+
+
+sequence_5 = (s34, s35, s36, s37, s38)
+mid_5 = 2
+compartments_5 = ["cyt"] * 4 + ["cell"]
+complexes_5 = [(0, 1), (2, 2), (3, 3), (4, 4)]
+pairs_5 = [(0, 2), (1, 3), (None, 4)]
+rate_5 = Rate("3.0*[K()::cyt]/2.0*v_1")
+
+r5 = Rule(sequence_5, mid_5, compartments_5, complexes_5, pairs_5, rate_5)
+
+sequence_6 = (s39, s35, s38, s40, s37)
+mid_6 = 3
+compartments_6 = ["cyt"] * 2 + ["cell"] + ["cyt"] * 2
+complexes_6 = [(0, 1), (2, 2), (3, 3), (4, 4)]
+pairs_6 = [(0, 3), (1, 4), (2, None)]
+rate_6 = Rate("3.0*[K(T{3+})::cyt]/2.0*v_1")
+
+r6 = Rule(sequence_6, mid_6, compartments_6, complexes_6, pairs_6, rate_6)
+rule_no_rate = Rule(sequence_4, mid_4, compartments_4, complexes_4, pairs_4, None)
+reversed_no_rate = Rule(
+    reversed_sequence_4, mid_4, compartments_4, reversed_complexes_4, pairs_4, None
+)
+
+sequence_c1 = (s34, s35, s36, s37, s2)
+mid_c1 = 2
+compartments_c1 = ["cyt"] * 5
+complexes_c1 = [(0, 0), (1, 1), (2, 3), (4, 4)]
+pairs_c1 = [(0, 2), (1, 3), (None, 4)]
+rate_c1 = Rate("3*[K()::cyt]/2*v_1")
+
+rule_c1 = Rule(sequence_c1, mid_c1, compartments_c1, complexes_c1, pairs_c1, rate_c1)
+
+rule_no_change = Rule(
+    sequence_no_change, mid_c1, compartments_c1, complexes_c1, pairs_c1, rate_c1
+)
+
+sequence_repl1 = (s31, s31, s31)
+mid_repl1 = 1
+compartments_repl1 = ["rep"] * 3
+complexes_repl1 = [(0, 0), (1, 1), (2, 2)]
+pairs_repl1 = [(0, 1), (0, 2)]
+rate_repl1 = Rate("3.0*[X()::rep]/2.0*v_1")
+
+rule_repl1 = Rule(
+    sequence_repl1, mid_repl1, compartments_repl1, complexes_repl1, pairs_repl1, None
+)
+rule_repl1_rate = Rule(
+    sequence_repl1,
+    mid_repl1,
+    compartments_repl1,
+    complexes_repl1,
+    pairs_repl1,
+    rate_repl1,
+)
+
+repl_sequence2 = (s31, s31, s31, s31)
+mid_repl2 = 1
+compartments_repl2 = ["rep"] * 4
+complexes_repl2 = [(0, 0), (1, 1), (2, 2), (3, 3)]
+pairs_repl2 = [(0, 1), (0, 2), (0, 3)]
+
+rule_repl2 = Rule(
+    repl_sequence2, mid_repl2, compartments_repl2, complexes_repl2, pairs_repl2, None
+)
+
+# reactions
+
+reaction1 = Reaction(lhs, rhs, rate_5)
+
+reaction_c1_1 = Reaction(
+    Side([c1_c1, c1_c3]),
+    Side([c1_c5, c1_c9]),
+    rate_c1,
+)
+reaction_c1_2 = Reaction(
+    Side([c1_c1, c1_c3]),
+    Side([c1_c5, c1_c10]),
+    rate_c1,
+)
+reaction_c1_3 = Reaction(
+    Side([c1_c2, cx3]),
+    Side([c1_c7, c1_c10]),
+    rate_c1,
+)
+reaction_c1_4 = Reaction(
+    Side([c1_c1, cx3]),
+    Side([c1_c6, c1_c9]),
+    rate_c1,
+)
+reaction_c1_5 = Reaction(
+    Side([c1_c2, c1_c3]),
+    Side([c1_c8, c1_c9]),
+    rate_c1,
+)
+reaction_c1_6 = Reaction(
+    Side([c1_c2, c1_c3]),
+    Side([c1_c8, c1_c10]),
+    rate_c1,
+)
+reaction_c1_7 = Reaction(
+    Side([c1_c1, cx3]),
+    Side([c1_c6, c1_c10]),
+    rate_c1,
+)
+reaction_c1_8 = Reaction(
+    Side([c1_c2, cx3]),
+    Side([c1_c7, c1_c9]),
+    rate_c1,
+)
+
+reactions_c1 = {
+    reaction_c1_1,
+    reaction_c1_2,
+    reaction_c1_3,
+    reaction_c1_4,
+    reaction_c1_5,
+    reaction_c1_6,
+    reaction_c1_7,
+    reaction_c1_8,
+}
diff --git a/Testing/parsing/__init__.py b/Testing/parsing/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/Testing/parsing/test_PCTL.py b/Testing/parsing/test_PCTL.py
new file mode 100644
index 0000000..0db9b4a
--- /dev/null
+++ b/Testing/parsing/test_PCTL.py
@@ -0,0 +1,17 @@
+import pytest
+
+from eBCSgen.Parsing.ParsePCTLformula import PCTLparser
+
+parser = PCTLparser()
+
+def test_parse():
+    formula = "P =? [F T(P{m})::x >= 2 & T(P{i})::x = 0]"
+    assert parser.parse(formula).success
+    formula = "P >= 0.5 [G T(P{m})::x = 2 & T(P{i})::x <= 0]"
+    assert parser.parse(formula).success
+    formula = "P <= 0.5 [T(P{m})::x < 2 U T(P{m})::x > 0]"
+    assert parser.parse(formula).success
+    formula = "P > 0.5 [F T(P{m})::x < 2 & ( T(P{m})::x = 0 | T()::x >= 7)]"
+    assert parser.parse(formula).success
+    formula = "P < 0.5 [G (T(P{m})::x = 2 & T(P{m})::x = 0) | (T(P{m})::x <= 2 & T(P{m})::x >= 0) ]"
+    assert parser.parse(formula).success
\ No newline at end of file
diff --git a/Testing/parsing/test_atomic.py b/Testing/parsing/test_atomic.py
new file mode 100644
index 0000000..0f519af
--- /dev/null
+++ b/Testing/parsing/test_atomic.py
@@ -0,0 +1,18 @@
+import Testing.objects_testing as objects
+
+
+def test_parser():
+    assert objects.atomic_parser.parse("T{s}").data == objects.a1
+    assert objects.atomic_parser.parse("S{a}").data == objects.a2
+    assert objects.atomic_parser.parse("S{s}").data == objects.a3
+    assert objects.atomic_parser.parse("T{_}").data == objects.a4
+    assert objects.atomic_parser.parse("S{_}").data == objects.a5
+    assert objects.atomic_parser.parse("T{p}").data == objects.a6
+    assert objects.atomic_parser.parse("T{u}").data == objects.a7
+
+    assert not objects.atomic_parser.parse("T{s").success
+    assert not objects.atomic_parser.parse("T{}").success
+    assert not objects.atomic_parser.parse("Ts}").success
+    assert not objects.atomic_parser.parse("{s}").success
+    assert not objects.atomic_parser.parse("x").success
+    assert not objects.atomic_parser.parse("").success
diff --git a/Testing/parsing/test_cmplx_in_abstr_seq.py b/Testing/parsing/test_cmplx_in_abstr_seq.py
new file mode 100644
index 0000000..068f47d
--- /dev/null
+++ b/Testing/parsing/test_cmplx_in_abstr_seq.py
@@ -0,0 +1,17 @@
+import pytest
+
+from Testing.models.get_model_str import get_model_str
+import Testing.objects_testing as objects
+
+
+def test_complexes_in_abstract_sequence():
+    # is allowed
+    model = get_model_str("model_cmplx_in_abstr_seq1")
+    ret1 = objects.model_parser.parse(model)
+    assert ret1.success
+
+    # should be allowed
+    model = get_model_str("model_cmplx_in_abstr_seq2")
+    ret2 = objects.model_parser.parse(model)
+    assert ret2.success
+    assert ret1.data == ret2.data
diff --git a/Testing/parsing/test_complex.py b/Testing/parsing/test_complex.py
new file mode 100644
index 0000000..e16dfb3
--- /dev/null
+++ b/Testing/parsing/test_complex.py
@@ -0,0 +1,69 @@
+import Testing.objects_testing as objects
+
+
+def test_parser():
+    ret = objects.rate_complex_parser.parse("B(T{s})::cell")
+    assert ret.success
+    assert ret.data.children[0] == objects.c1
+
+    ret = objects.rate_complex_parser.parse("B(T{s}).D().K(T{s},S{s},U{a})::cell")
+    assert ret.success
+    assert ret.data.children[0] == objects.c2
+
+    ret = objects.rate_complex_parser.parse(
+        "B(T{s}).K(T{s}, S{s}, U{a}).D(S{_},T{p})::cyt"
+    )
+    assert ret.success
+    assert ret.data.children[0] == objects.c3
+
+    ret = objects.rate_complex_parser.parse("B(T{s}).T{s}::cyt")
+    assert ret.success
+    assert ret.data.children[0] == objects.c4
+
+    ret = objects.rate_complex_parser.parse("D().D().D().B(T{_}).B(T{_}).B(T{_})::cell")
+    assert ret.success
+    assert ret.data.children[0] == objects.c5
+
+    ret = objects.rate_complex_parser.parse("K().K().T{p}::cyt")
+    assert ret.success
+    assert ret.data.children[0] == objects.c6
+
+    ret = objects.rate_complex_parser.parse("T{u}.T{_}::cell")
+    assert ret.success
+    assert ret.data.children[0] == objects.c7
+
+    ret = objects.rate_complex_parser.parse("(T{s})::cell")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("()::cell")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("x::cell")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("BT{s})::cell")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("B(T{s}::cell")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("B(T{s}::cell)")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("B(T{})::cell")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("B(T{s}))::cell")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("B(T{s})::")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("B(T{s}, T{_})::cell")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("B(T{s}, T{s})::cell")
+    assert not ret.success
+
+    ret = objects.rate_complex_parser.parse("B(T{s}, T{a})::cell")
+    assert not ret.success
diff --git a/Testing/parsing/test_insert_atomic_to_complex.py b/Testing/parsing/test_insert_atomic_to_complex.py
new file mode 100644
index 0000000..c32f5d8
--- /dev/null
+++ b/Testing/parsing/test_insert_atomic_to_complex.py
@@ -0,0 +1,42 @@
+import pytest
+
+
+from eBCSgen.Errors.ComplexParsingError import ComplexParsingError
+from eBCSgen.Parsing.ParseBCSL import Parser, TransformAbstractSyntax
+
+
+@pytest.mark.parametrize('string_atomic, string_value, string_expected', [
+    ["A{i}", "B(S{p}).A{_}.A{d}", "B(S{p}).A{i}.A{d}"],
+    ["A{i}", "B(S{p}).A{d}.A{_}", "B(S{p}).A{d}.A{i}"]
+])
+def test_insert_atomic_to_complex_correct_input(string_atomic, string_value, string_expected):
+    parser_atomic = Parser("atomic")
+    parser_value = Parser("value")
+
+    result_atomic = parser_atomic.syntax_check(string_atomic).data.children[0]
+    result_value = parser_value.syntax_check(string_value).data.children[0]
+
+    transformer = TransformAbstractSyntax(set())
+    transformed_value = transformer.insert_atomic_to_complex(result_atomic, result_value)
+
+    result_expected = parser_value.syntax_check(string_expected).data.children[0]
+
+    assert transformed_value == result_expected
+
+
+@pytest.mark.parametrize('string_atomic, string_value', [
+    ["A{i}", "B(S{p}).A{a}.A{d}"],
+    ["A{_}", "B(S{p}).A{d}.A{a}"],
+    ["A{i}", "B(S{p}).D{d}"]
+])
+def test_insert_atomic_to_complex_incorrect_input(string_atomic, string_value):
+    parser_atomic = Parser("atomic")
+    parser_value = Parser("value")
+
+    result_atomic = parser_atomic.syntax_check(string_atomic).data.children[0]
+    result_value = parser_value.syntax_check(string_value).data.children[0]
+
+    transformer = TransformAbstractSyntax(set())
+
+    with pytest.raises(ComplexParsingError):
+        _ = transformer.insert_atomic_to_complex(result_atomic, result_value)
diff --git a/Testing/parsing/test_insert_atomic_to_struct.py b/Testing/parsing/test_insert_atomic_to_struct.py
new file mode 100644
index 0000000..d5d15dd
--- /dev/null
+++ b/Testing/parsing/test_insert_atomic_to_struct.py
@@ -0,0 +1,41 @@
+import pytest
+
+
+from eBCSgen.Errors.ComplexParsingError import ComplexParsingError
+from eBCSgen.Parsing.ParseBCSL import Parser, TransformAbstractSyntax
+
+
+@pytest.mark.parametrize('string_atomic, string_structure, string_expected', [
+    ["A{i}", "B(S{p})", "B(S{p},A{i})"],
+    ["A{i}", "B()", "B(A{i})"]
+])
+def test_insert_atomic_to_struct_correct_input(string_atomic, string_structure, string_expected):
+    parser_atomic = Parser("atomic")
+    parser_structure = Parser("structure")
+
+    result_atomic = parser_atomic.syntax_check(string_atomic).data.children[0]
+    result_structure = parser_structure.syntax_check(string_structure).data.children[0]
+
+    transformer = TransformAbstractSyntax(set())
+    transformed_value = transformer.insert_atomic_to_struct(result_atomic, result_structure)
+
+    result_expected = parser_structure.syntax_check(string_expected).data.children[0]
+
+    assert transformed_value == result_expected
+
+
+@pytest.mark.parametrize('string_atomic, string_structure', [
+    ["A{i}", "B(A{a})"],
+    ["A{i}", "B(A{i})"]
+])
+def test_insert_atomic_to_struct_incorrect_input(string_atomic, string_structure):
+    parser_atomic = Parser("atomic")
+    parser_structure = Parser("structure")
+
+    result_atomic = parser_atomic.syntax_check(string_atomic).data.children[0]
+    result_structure = parser_structure.syntax_check(string_structure).data.children[0]
+
+    transformer = TransformAbstractSyntax(set())
+
+    with pytest.raises(ComplexParsingError):
+        _ = transformer.insert_atomic_to_struct(result_atomic, result_structure)
diff --git a/Testing/parsing/test_insert_struct_to_complex.py b/Testing/parsing/test_insert_struct_to_complex.py
new file mode 100644
index 0000000..f028fb3
--- /dev/null
+++ b/Testing/parsing/test_insert_struct_to_complex.py
@@ -0,0 +1,44 @@
+import pytest
+
+
+from eBCSgen.Errors.ComplexParsingError import ComplexParsingError
+from eBCSgen.Parsing.ParseBCSL import Parser, TransformAbstractSyntax
+
+
+@pytest.mark.parametrize('string_structure, string_value, string_expected', [
+    ["B(S{i})", "B().A{_}.B(S{a}).A{d}", "B(S{i}).A{_}.B(S{a}).A{d}"],
+    ["B()", "B().A{_}.B(S{a}).A{d}", "B().A{_}.B(S{a}).A{d}"],
+    ["B(S{i})", "B(S{a}).A{_}.B().A{d}", "B(S{a}).A{_}.B(S{i}).A{d}"],
+    ["B(S{i})", "B(T{a}).A{_}.B().A{d}", "B(T{a},S{i}).A{_}.B().A{d}"]
+])
+def test_insert_struct_to_complex_correct_input(string_structure, string_value, string_expected):
+    parser_structure = Parser("structure")
+    parser_value = Parser("value")
+
+    result_structure = parser_structure.syntax_check(string_structure).data.children[0]
+    result_value = parser_value.syntax_check(string_value).data.children[0]
+
+    transformer = TransformAbstractSyntax(set())
+    transformed_value = transformer.insert_struct_to_complex(result_structure, result_value)
+
+    result_expected = parser_value.syntax_check(string_expected).data.children[0]
+
+    assert transformed_value == result_expected
+
+
+@pytest.mark.parametrize('string_structure, string_value', [
+    ["B()", "F(S{p}).A{a}.A{d}"],
+    ["B(S{i})", "B(S{p}).A{d}.A{a}"],
+    ["B(S{i})", "B(S{p}).A{d}.B(S{i}).A{a}"]
+])
+def test_insert_struct_to_complex_incorrect_input(string_structure, string_value):
+    parser_structure = Parser("structure")
+    parser_value = Parser("value")
+
+    result_structure = parser_structure.syntax_check(string_structure).data.children[0]
+    result_value = parser_value.syntax_check(string_value).data.children[0]
+
+    transformer = TransformAbstractSyntax(set())
+
+    with pytest.raises(ComplexParsingError):
+        _ = transformer.insert_struct_to_complex(result_structure, result_value)
diff --git a/Testing/parsing/test_model.py b/Testing/parsing/test_model.py
new file mode 100644
index 0000000..dfeae4a
--- /dev/null
+++ b/Testing/parsing/test_model.py
@@ -0,0 +1,46 @@
+import collections
+import pytest
+
+from eBCSgen.Core.Model import Model
+
+from Testing.models.get_model_str import get_model_str
+import Testing.objects_testing as objects
+
+
+# inits
+
+inits = collections.Counter({objects.c27: 2, objects.c28: 1})
+
+# defs
+
+defs = {'k1': 0.05, 'k2': 0.12}
+
+model = Model({objects.r1, objects.r2, objects.r3}, inits, defs, set())
+# model
+
+
+model_str_1 = get_model_str("model1")
+
+model_str_2 = get_model_str("model2")
+
+model_with_comments_str = get_model_str("model_with_comments")
+
+model_wrong_1 = get_model_str("model_wrong1")
+
+model_wrong_2 = get_model_str("model_wrong2")
+
+def test_parser():
+    assert objects.model_parser.parse(model_str_1).data == model
+
+
+def test_parser_errors():
+    assert not objects.model_parser.parse(model_wrong_1).success
+
+    assert not objects.model_parser.parse(model_wrong_2).success
+        
+
+def test_comments():
+    model_with_comments = objects.model_parser.parse(model_with_comments_str).data
+    model_without_comments = objects.model_parser.parse(model_str_2).data
+
+    assert model_with_comments == model_without_comments
\ No newline at end of file
diff --git a/Testing/parsing/test_observables.py b/Testing/parsing/test_observables.py
new file mode 100644
index 0000000..75340a9
--- /dev/null
+++ b/Testing/parsing/test_observables.py
@@ -0,0 +1,67 @@
+import pytest
+
+import Testing.objects_testing as objects
+
+
+def test_parser():
+    observable_expr1 = "abc: A()::cell"
+    assert objects.observable_parser.parse(observable_expr1)
+
+    observable_expr2 = "efg: E(F{_})::cell"
+    assert objects.observable_parser.parse(observable_expr2)
+
+    observable_expr3 = "hij: H()::cell"
+    assert objects.observable_parser.parse(observable_expr3)
+
+    observable_expr4 = "klm: K()::cyt * L()::cell + M()::cell"
+    assert objects.observable_parser.parse(observable_expr4)
+
+    observable_expr5 = "nop: N()::cell"
+    assert objects.observable_parser.parse(observable_expr5).success
+
+    observable_expr6 = "qrs: Q().R().S()::cell"
+    assert objects.observable_parser.parse(observable_expr6).success
+
+    observable_expr7 = "tuv: T(U{v})::cell "
+    assert objects.observable_parser.parse(observable_expr7).success
+
+    observable_expr8 = "wx: 2 * W{x}::cell"
+    assert objects.observable_parser.parse(observable_expr8).success
+
+    observable_expr9 = "z: Y{z}::cyt + Z{y}::ext"
+    assert objects.observable_parser.parse(observable_expr9).success
+
+    observable_expr10 = "z: 2 * Y{z}::cyt + Z{y}::ext ** 2"
+    assert objects.observable_parser.parse(observable_expr10).success
+
+    observable_expr10 = "z: (Y{z}::cell + Z{y}::cyt) / 2.1 ** 10"
+    assert objects.observable_parser.parse(observable_expr10).success
+
+    observable_expr11 = "scaled_A: 1000 * A{i}::cell"
+    assert objects.observable_parser.parse(observable_expr11).success
+
+    observable_expr12 = "obs_A_all: A{i}::cell + A{a}::cell"
+    assert objects.observable_parser.parse(observable_expr12).success
+
+    observables_expr = (
+        "#! observables\n"
+        + observable_expr1
+        + "\n"
+        + observable_expr2
+        + "\n"
+        + observable_expr3
+        + "\n"
+        + observable_expr4
+        + "\n"
+        + observable_expr5
+        + "\n"
+        + observable_expr6
+    )
+    assert objects.observables_parser.parse(observables_expr).success
+
+    assert not objects.observable_parser.parse("A()::cell > 2").success
+    assert not objects.observable_parser.parse("a: A(::cell").success
+    assert not objects.observable_parser.parse("a: b: A():cell > 2").success
+    assert not objects.observable_parser.parse("a: 2 > A():cell").success
+    assert not objects.observable_parser.parse("a: A()::cell$").success
+    assert not objects.observable_parser.parse("a: A{}::cell").success
diff --git a/Testing/parsing/test_rate.py b/Testing/parsing/test_rate.py
new file mode 100644
index 0000000..39efb9d
--- /dev/null
+++ b/Testing/parsing/test_rate.py
@@ -0,0 +1,98 @@
+import Testing.objects_testing as objects
+
+
+def test_parser():
+    ret = objects.rate_parser.parse("1")
+    assert ret.success
+
+    ret = objects.rate_parser.parse("0.5")
+    assert ret.success
+
+    ret = objects.rate_parser.parse("2*3")
+    assert ret.success
+
+    ret = objects.rate_parser.parse("2**3")
+    assert ret.success
+
+    ret = objects.rate_parser.parse("2*[K()::cyt]")
+    assert ret.success
+
+    ret = objects.rate_parser.parse("[K()::cyt]**2")
+    assert ret.success
+
+    ret = objects.rate_parser.parse(
+        "(2 * [K()::cyt] + 3 * [B(T{s})::cell]) / (4 * [T{u}.T{_}::cell] - 2 * [B(T{s})::cyt])"
+    )
+    assert ret.success
+
+    expr1 = "([B(T{s})::cell] + [B(T{s})::cyt]) / ([B(T{_}).C()::cell] + [B()::cell])"
+    expr2 = (
+        "([B(T{s}).T{s}::cyt] - [T{_}::cell]) / ([T{u}.T{_}::cell] - [B(T{s})::cell])"
+    )
+    ret = objects.rate_parser.parse(expr1)
+    assert ret.success
+    ret = objects.rate_parser.parse(expr2)
+    assert ret.success
+
+    # currently does not pass
+    ret = objects.rate_parser.parse(f"({expr1}/{expr2})")
+    assert ret.success
+
+    ret = objects.rate_parser.parse("0,5")
+    assert not ret.success
+
+    # missing []
+    ret = objects.rate_parser.parse("2*K()::cyt")
+    assert not ret.success
+
+    # invalid syntax
+    ret = objects.rate_parser.parse("2 +")
+    assert not ret.success
+
+    # unsupported operators
+    ret = objects.rate_parser.parse("2 & 3")
+    assert not ret.success
+
+    # mismatched parentheses
+    ret = objects.rate_parser.parse("(2 * 3")
+    assert not ret.success
+
+    # incorrect format
+    ret = objects.rate_parser.parse("2,,3")
+    assert not ret.success
+
+    # nested complexes
+    ret = objects.rate_parser.parse("[K([L()::cell])::cyt]")
+    assert not ret.success
+
+    # nested parenthesis
+    assert objects.rate_parser.parse("(2 * (3 + 4)) / (5 - (6 * 7))").success
+    assert objects.rate_parser.parse("((2 + 3) * (4 - 5)) / ((6 * 7) + 8)").success
+
+    # complex mathematical operations
+    assert objects.rate_parser.parse("2 * [A()::cyt] + 3 / [B(T{s})::cell]").success
+    assert objects.rate_parser.parse("[X()::cyt] ** 2 - (4 * [Y()::cell])").success
+
+    # multiple level nesting
+    assert objects.rate_parser.parse("((2 * [A()::cyt]) + (3 / [B(T{s})::cell])) / ((4 * [C()::cell]) - (2 * [D()::cyt]))").success
+    assert objects.rate_parser.parse("([X()::cyt] ** 2) / ([Y()::cell] + ([Z()::cyt] * 3))").success
+    assert objects.rate_parser.parse("1/([X()::cyt] ** 2) + ([Y()::cell] + ([Z()::cyt] * 3))/(5 + [KaiC()::cyt])").success
+    
+    # complex rate_agents
+    assert objects.rate_parser.parse("[A().B()::cyt] + [C().D()::nuc]").success
+    assert objects.rate_parser.parse("[X(Y{a}).Z(W{b})::cell]").success
+    assert objects.rate_parser.parse("2 * [A()::cyt] + [B(T{s})::cell] ** 3").success
+    assert objects.rate_parser.parse("[X()::cyt] / ([Y{a}::cell] + [Z{b}::cyt])").success
+
+    # invalid syntax
+    assert not objects.rate_parser.parse("2 * [A()::cyt] + 3 / [B(T{s})::cell] +").success # extra +
+    assert not objects.rate_parser.parse("[A().B()::cyt] + C().D()::nuc]").success # missing [
+    assert not objects.rate_parser.parse("[X(Y{a}).Z(W{b})::cell").success # missing ]
+    assert not objects.rate_parser.parse("2 * [A()::cyt] + [B(T{s})::cell * 3").success # missing ]
+    assert not objects.rate_parser.parse("[X()::cyt] ** 2 - (4 * [Y()::cell]))").success # extra )
+    assert not objects.rate_parser.parse("2 & [A()::cyt] + 3 | [B(T{s})::cell]").success # invalid operators
+    assert not objects.rate_parser.parse("2 ** [X()::cyt]").success # invalid operation - exponentiation
+    assert not objects.rate_parser.parse("2 * K()::cyt").success # missing required brackets for rate agent
+    assert not objects.rate_parser.parse("[X()::cyt] ** 2 - 4 * Y()::cell").success # missing required brackets for rate agent
+    assert not objects.rate_parser.parse("[X()::cyt] / {Y} + [Z()::cell]").success
+    assert not objects.rate_parser.parse("[K([L()::cell])::cyt] + [A(B(C{d})).D(E{e}::cyt)]").success # nested rate agents and incorrect syntax
diff --git a/Testing/test_remote_parsing.py b/Testing/parsing/test_remote_parsing.py
similarity index 56%
rename from Testing/test_remote_parsing.py
rename to Testing/parsing/test_remote_parsing.py
index e8ac8bf..63af865 100644
--- a/Testing/test_remote_parsing.py
+++ b/Testing/parsing/test_remote_parsing.py
@@ -1,51 +1,17 @@
 import unittest
 import requests
 
-from eBCSgen.Parsing.ParseBCSL import Parser
+from Testing.models.get_model_str import get_model_str
+import Testing.objects_testing as objects
 
 
 class TestModel(unittest.TestCase):
     def setUp(self):
-        self.model_parser = Parser("model")
         self.url = 'http://biodivine-vm.fi.muni.cz/BCSLparser/'
 
-        self.model_wrong_1 = \
-            """#! rules
-            X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep]
-            X(T{a})::rep => X(T{o}):;rep @ k2*[Z()::rep]
-            => Y(P{f})::rep @ 1/(1+([X()::rep])**4)
+        self.model_wrong_1 = get_model_str("model_wrong1")
 
-            #! inits
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-            Y(P{g}, N{l})::rep
-
-            #! definitions
-            k1 = 0.05
-            k2 = 0.12
-            """
-
-        self.model_with_complexes = """
-            #! rules
-            // commenting
-            X(T{a}):XX::rep => X(T{o}):XX::rep @ k2*[X().X()::rep]
-            K{i}:X():XYZ::rep => K{p}:X():XYZ::rep @ k1*[X().Y().Z()::rep] // also here
-            => P{f}:XP::rep @ 1/(1+([X().P{_}::rep])**4) // ** means power (^)
-
-            #! inits
-            // here
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-            Y(P{g}, N{l})::rep // comment just 1 item
-
-            #! definitions
-            // and
-            k1 = 0.05 // also
-            k2 = 0.12
-
-            #! complexes
-            XYZ = X().Y().Z() // a big complex
-            XX = X().X()
-            XP = X().P{_}
-            """
+        self.model_with_complexes = get_model_str("model_with_complexes.txt")
 
     def test_remote_request(self):
         try:
@@ -57,7 +23,7 @@ def test_remote_request(self):
     def test_remote_parsing_model(self):
         try:
             # correct one
-            result_local = self.model_parser.syntax_check(self.model_with_complexes)
+            result_local = objects.model_parser.syntax_check(self.model_with_complexes)
 
             headers = {'value-type': 'application/json'}
 
@@ -69,7 +35,7 @@ def test_remote_parsing_model(self):
             self.assertEqual(result_remote['success'], result_local.success)
 
             # wrong one
-            result_local = self.model_parser.syntax_check(self.model_wrong_1)
+            result_local = objects.model_parser.syntax_check(self.model_wrong_1)
 
             response = requests.post(self.url + 'parse',
                                      json={'start': 'model', 'expression': self.model_wrong_1},
diff --git a/Testing/parsing/test_rule.py b/Testing/parsing/test_rule.py
new file mode 100644
index 0000000..3ce55df
--- /dev/null
+++ b/Testing/parsing/test_rule.py
@@ -0,0 +1,97 @@
+import pytest
+
+import Testing.objects_testing as objects
+from eBCSgen.Core.Rate import Rate
+
+
+def test_parser():
+    rule_expr = "K(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt + D(B{_})::cell @ 3*[K()::cyt]/2*v_1"
+    assert objects.rule_parser.parse(rule_expr).data[1] == objects.r5
+
+    rule_expr = "K(B{-}).B()::cyt + D(B{_})::cell => K(B{+})::cyt + B()::cyt @ 3*[K(T{3+})::cyt]/2*v_1"
+    assert objects.rule_parser.parse(rule_expr).data[1] == objects.r6
+
+    rule_expr = "X()::rep => @ k1*[X()::rep]"
+    assert objects.rule_parser.parse(rule_expr).data[1] == objects.r1
+
+    rule_expr = "=> Y()::rep @ 1/(1+([X()::rep])**4)"
+    assert objects.rule_parser.parse(rule_expr).data[1] == objects.r3
+
+    rule_expr = "K(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt"
+    assert objects.rule_parser.parse(rule_expr).data[1] == objects.rule_no_rate
+
+
+def test_bidirectional():
+    rule_expr = "#! rules\nK(S{u}).B()::cyt <=> K(S{p})::cyt + B()::cyt"
+    parsed = objects.rules_parser.parse(rule_expr)
+    assert parsed.success
+    assert objects.rule_no_rate in parsed.data["rules"]
+    assert objects.reversed_no_rate in parsed.data["rules"]
+
+    rule_expr = "#! rules\nK(S{u}).B()::cyt <=> K(S{p})::cyt + B()::cyt @ 3*[K()::cyt]/2*v_1 | 3*[K()::cyt]/2*v_1"
+    parsed = objects.rules_parser.parse(rule_expr)
+    assert parsed.success
+    assert objects.r4 in parsed.data["rules"]
+    assert objects.reversed_r4a in parsed.data["rules"]
+
+    rule_expr = "#! rules\nK(S{u}).B()::cyt <=> K(S{p})::cyt + B()::cyt @ 3*[K()::cyt]/2*v_1 | 2*[K()::cyt]/3*v_1"
+    parsed = objects.rules_parser.parse(rule_expr)
+    assert parsed.success
+    assert objects.r4 in parsed.data["rules"]
+    assert objects.reversed_r4b in parsed.data["rules"]
+
+    rule_expr = "#! rules\n <=> K(S{p})::cyt + B()::cyt @ 3*[K()::cyt]/2*v_1 | 3*[K()::cyt]/2*v_1"
+    parsed = objects.rules_parser.parse(rule_expr)
+    assert parsed.success
+    assert objects.one_side_bidirectional_a in parsed.data["rules"]
+    assert objects.one_side_bidirectional_b in parsed.data["rules"]
+
+    rule_expr = "#! rules\n K(S{p})::cyt + B()::cyt <=> @ 3*[K()::cyt]/2*v_1 | 3*[K()::cyt]/2*v_1"
+    parsed = objects.rules_parser.parse(rule_expr)
+    assert parsed.success
+    assert objects.one_side_bidirectional_a in parsed.data["rules"]
+    assert objects.one_side_bidirectional_b in parsed.data["rules"]
+
+    rule_expr = "#! rules\n K(S{p})::cyt + B()::cyt <=> @ 3*[K()::cyt]/2*v_1 | 2*[K()::cyt]/3*v_1"
+    parsed = objects.rules_parser.parse(rule_expr)
+    assert parsed.success
+    assert objects.one_side_bidirectional_a in parsed.data["rules"]
+    assert objects.one_side_bidirectional_b_reversed_rate in parsed.data["rules"]
+
+    rule_expr = "#! rules\n K(S{p})::cyt + B()::cyt <=>"
+    parsed = objects.rules_parser.parse(rule_expr)
+    assert parsed.success
+    assert objects.one_side_bidirectional_a_no_rate in parsed.data["rules"]
+    assert objects.one_side_bidirectional_b_no_rate in parsed.data["rules"]
+
+    rule_expr = (
+        "#! rules\nK(S{u}).B()::cyt <=> K(S{p})::cyt + B()::cyt @ 3*[K()::cyt]/2*v_1"
+    )
+    assert not objects.rules_parser.parse(rule_expr).success
+
+    rule_expr = "#! rules\nK(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt @ 3*[K()::cyt]/2*v_1 | 2*[K()::cyt]/3*v_1"
+    assert not objects.rules_parser.parse(rule_expr).success
+
+
+def test_replication():
+    rule_expr = "X()::rep =*> X()::rep + X()::rep"
+    result = objects.rule_parser.parse(rule_expr)
+    assert result.success
+    assert result.data[1] == objects.rule_repl1
+
+    rule_expr = "X()::rep =*> X()::rep + X()::rep @ 3*[X()::rep]/2*v_1"
+    result = objects.rule_parser.parse(rule_expr)
+    assert result.success
+    rate_repl1 = Rate("3.0*[X()::rep]/2*v_1")
+    assert result.data[1] == objects.rule_repl1_rate
+
+    rule_expr = "X()::rep =*> X()::rep + X()::rep + X()::rep"
+    result = objects.rule_parser.parse(rule_expr)
+    assert result.success
+    assert result.data[1] == objects.rule_repl2
+
+    rule_expr = "X()::rep + Y()::rep =*> X()::rep + X()::rep"
+    assert not objects.rule_parser.parse(rule_expr).success
+
+    rule_expr = "X()::rep =*> X()::rep + X()::rep + Y()::rep"
+    assert not objects.rule_parser.parse(rule_expr).success
diff --git a/Testing/parsing/test_side.py b/Testing/parsing/test_side.py
new file mode 100644
index 0000000..908aa9b
--- /dev/null
+++ b/Testing/parsing/test_side.py
@@ -0,0 +1,56 @@
+import Testing.objects_testing as objects
+
+
+def test_parser():
+    ret = objects.side_parser.parse("B(T{s})::cell")
+    assert ret.success
+    assert ret.data.to_side() == objects.side1
+
+    ret = objects.side_parser.parse(
+        "B(T{s}).T{s}::cyt + D().D().D().B(T{_}).B(T{_}).B(T{_})::cell"
+    )
+    assert ret.success
+    assert ret.data.to_side() == objects.side2
+
+    ret = objects.side_parser.parse(
+        "B(T{s})::cell + B(T{s}).D().K(T{s},S{s},U{a})::cell + B(T{s}).D().K(T{s},S{s},U{a})::cell"
+    )
+    assert ret.success
+    assert ret.data.to_side() == objects.side3
+
+    ret = objects.side_parser.parse(
+        "B(T{s})::cell + 2 B(T{s}).D().K(T{s},S{s},U{a})::cell"
+    )
+    assert ret.success
+    assert ret.data.to_side() == objects.side3
+
+    ret = objects.side_parser.parse(
+        "T{u}.T{_}::cell + K().K().T{p}::cyt + D().D().D().B(T{_}).B(T{_}).B(T{_})::cell"
+    )
+    assert ret.success
+    assert ret.data.to_side() == objects.side4
+
+    ret = objects.side_parser.parse("")
+    assert ret.success
+    assert ret.data.to_side() == objects.side5
+
+    ret = objects.side_parser.parse("(T{s})::cell")
+    assert not ret.success
+
+    ret = objects.side_parser.parse("BT{s})::cell")
+    assert not ret.success
+
+    ret = objects.side_parser.parse("B(T{s}::cell")
+    assert not ret.success
+
+    ret = objects.side_parser.parse("B(T{s})::")
+    assert not ret.success
+
+    ret = objects.side_parser.parse("B(T{s}")
+    assert not ret.success
+
+    # not unique atomics in structure
+    ret = objects.side_parser.parse(
+        "B(T{s})::cell + B(T{s}).D().K(T{s},S{s},S{_})::cell + B(T{s}).D().K(T{s},S{s},S{_})::cell"
+    )
+    assert not ret.success
diff --git a/Testing/parsing/test_structure.py b/Testing/parsing/test_structure.py
new file mode 100644
index 0000000..a7fbfd4
--- /dev/null
+++ b/Testing/parsing/test_structure.py
@@ -0,0 +1,23 @@
+import Testing.objects_testing as objects
+
+
+def test_parser():
+    assert objects.structure_parser.parse("B(T{s})").data == objects.s1
+    assert objects.structure_parser.parse("D()").data == objects.s2
+    assert objects.structure_parser.parse("K(T{s}, S{s}, U{a})").data == objects.s3
+    assert objects.structure_parser.parse("B(T{_})").data == objects.s4
+    assert objects.structure_parser.parse("D(S{_},T{p})").data == objects.s5
+    assert objects.structure_parser.parse("K()").data == objects.s6
+    assert objects.structure_parser.parse("K(S{a},T{_})").data == objects.s7
+
+    assert not objects.structure_parser.parse("B(T{})").success
+    assert not objects.structure_parser.parse("(T{s})").success
+    assert not objects.structure_parser.parse("BT{s})").success
+    assert not objects.structure_parser.parse("B(T{s}").success
+    assert not objects.structure_parser.parse("(B(T{s}))").success
+    assert not objects.structure_parser.parse("[B(T{s})]").success
+    assert not objects.structure_parser.parse("").success
+    assert not objects.structure_parser.parse("B({s})").success
+    assert not objects.structure_parser.parse("B(S{s}, S{a})").success
+    assert not objects.structure_parser.parse("B(S{a}, S{a})").success
+    assert not objects.structure_parser.parse("B(S{_}, S{a})").success
diff --git a/Testing/test_PCTL.py b/Testing/test_PCTL.py
index 174162b..58e6fd8 100644
--- a/Testing/test_PCTL.py
+++ b/Testing/test_PCTL.py
@@ -31,18 +31,6 @@ def test_str(self):
         formula = self.parser.parse(self.formula_3)
         self.assertEqual(self.formula_3, str(formula))
 
-    def test_parse(self):
-        formula = "P =? [F T(P{m})::x >= 2 & T(P{i})::x = 0]"
-        self.assertTrue(self.parser.parse(formula).success)
-        formula = "P >= 0.5 [G T(P{m})::x = 2 & T(P{i})::x <= 0]"
-        self.assertTrue(self.parser.parse(formula).success)
-        formula = "P <= 0.5 [T(P{m})::x < 2 U T(P{m})::x > 0]"
-        self.assertTrue(self.parser.parse(formula).success)
-        formula = "P > 0.5 [F T(P{m})::x < 2 & ( T(P{m})::x = 0 | T()::x >= 7)]"
-        self.assertTrue(self.parser.parse(formula).success)
-        formula = "P < 0.5 [G (T(P{m})::x = 2 & T(P{m})::x = 0) | (T(P{m})::x <= 2 & T(P{m})::x >= 0) ]"
-        self.assertTrue(self.parser.parse(formula).success)
-
     def test_get_complexes(self):
         formula = self.parser.parse(self.formula_1)
         self.assertEqual(formula.get_complexes(), [self.complex_1])
diff --git a/Testing/test_TS.py b/Testing/test_TS.py
index 526d184..99c1dcd 100644
--- a/Testing/test_TS.py
+++ b/Testing/test_TS.py
@@ -7,22 +7,14 @@
 from eBCSgen.TS.State import State, Vector, Memory
 from eBCSgen.TS.TransitionSystem import TransitionSystem
 
+import Testing.objects_testing as objects
+
 
 class TestTransitionSystem(unittest.TestCase):
     def setUp(self):
-        self.str1 = StructureAgent("X", set())
-        self.str2 = StructureAgent("Y", set())
-        self.str3 = StructureAgent("Z", set())
-        self.str4 = StructureAgent("W", set())
-
-        self.c1 = Complex([self.str1], "rep")
-        self.c2 = Complex([self.str2], "rep")
-        self.c3 = Complex([self.str3], "rep")
-        self.c4 = Complex([self.str4], "rep")
-
-        ordering = (self.c1, self.c2, self.c3)
-        ordering_wrong = (self.c1, self.c2, self.c3, self.c4),
-        ordering_reordered = (self.c3, self.c1, self.c2)
+        ordering = (objects.c27, objects.c28, objects.c29)
+        ordering_wrong = (objects.c27, objects.c28, objects.c29, objects.c30),
+        ordering_reordered = (objects.c29, objects.c27, objects.c28)
 
         self.s1 = State(Vector(np.array((1, 2, 3))), Memory(0))
         self.s2 = State(Vector(np.array((1, 2, 4))), Memory(0))
@@ -69,7 +61,7 @@ def setUp(self):
 
         # test bigger TS
 
-        ordering = (self.c1, self.c2, self.c3, self.c4)
+        ordering = (objects.c27, objects.c28, objects.c29, objects.c30)
 
         self.ts_bigger = TransitionSystem(ordering, 5)
         self.ts_bigger.states_encoding = {1: self.s1, 2: self.s2, 0: self.s3, 3: self.s4}
@@ -106,7 +98,7 @@ def test_sort_edges(self):
         self.assertEqual(sorted(edges), edges_sorted)
 
     def test_change_hell(self):
-        ordering = (self.c1, self.c2, self.c3, self.c4)
+        ordering = (objects.c27, objects.c28, objects.c29, objects.c30)
         ts = TransitionSystem(ordering, 4)
         ts.states_encoding = {1: self.s1, 2: self.s2, 0: self.s3, 3: self.hell}
         ts.change_hell()
diff --git a/Testing/test_atomic.py b/Testing/test_atomic.py
index f9cb7a1..16529bf 100644
--- a/Testing/test_atomic.py
+++ b/Testing/test_atomic.py
@@ -1,44 +1,41 @@
 import unittest
-from eBCSgen.Core.Atomic import AtomicAgent
+import Testing.objects_testing as objects
 
 
 class TestAtomic(unittest.TestCase):
-    def setUp(self):
-        self.a1 = AtomicAgent("T", "s")
-        self.a2 = AtomicAgent("S", "a")
-        self.a3 = AtomicAgent("T", "s")
-        self.a4 = AtomicAgent("T", "_")
-        self.a5 = AtomicAgent("T", "_")
-        self.a6 = AtomicAgent("T", "p")
-        self.a7 = AtomicAgent("T", "u")
-
     def test_eq(self):
-        self.assertEqual(self.a1, self.a3)
-        self.assertNotEqual(self.a1, self.a2)
-        self.assertNotEqual(self.a1, self.a4)
+        self.assertEqual(objects.a1, objects.a8)
+        self.assertNotEqual(objects.a1, objects.a2)
+        self.assertNotEqual(objects.a1, objects.a4)
 
     def test_print(self):
-        self.assertEqual(str(self.a1), "T{s}")
-        self.assertEqual(str(self.a4), "T{_}")
+        self.assertEqual(str(objects.a1), "T{s}")
+        self.assertEqual(str(objects.a4), "T{_}")
 
     def test_compatibility(self):
-        self.assertTrue(self.a4.compatible(self.a1))
-        self.assertFalse(self.a2.compatible(self.a1))
-        self.assertFalse(self.a1.compatible(self.a4))
+        self.assertTrue(objects.a4.compatible(objects.a1))
+        self.assertFalse(objects.a2.compatible(objects.a1))
+        self.assertFalse(objects.a1.compatible(objects.a4))
 
     def test_add_context(self):
         atomic_signature = {"T": {"u", "p"}}
         structure_signature = dict()
-        self.assertEqual(self.a4.add_context(self.a5, atomic_signature, structure_signature),
-                         {(self.a6, self.a6), (self.a7, self.a7)})
-        self.assertEqual(self.a6.add_context(self.a6, atomic_signature, structure_signature),
-                         {(self.a6, self.a6)})
-        self.assertEqual(self.a4.add_context(-1, atomic_signature, structure_signature),
-                         {(None, self.a6), (None, self.a7)})
+        self.assertEqual(
+            objects.a4.add_context(objects.a4, atomic_signature, structure_signature),
+            {(objects.a6, objects.a6), (objects.a7, objects.a7)},
+        )
+        self.assertEqual(
+            objects.a6.add_context(objects.a6, atomic_signature, structure_signature),
+            {(objects.a6, objects.a6)},
+        )
+        self.assertEqual(
+            objects.a4.add_context(-1, atomic_signature, structure_signature),
+            {(None, objects.a6), (None, objects.a7)},
+        )
 
     def test_reduce_context(self):
-        self.assertEqual(self.a3.reduce_context(), self.a4)
+        self.assertEqual(objects.a8.reduce_context(), objects.a4)
 
     def test_replace(self):
-        self.assertEqual(self.a4.replace(self.a6), self.a6)
-        self.assertEqual(self.a6.replace(self.a7), self.a6)
+        self.assertEqual(objects.a4.replace(objects.a6), objects.a6)
+        self.assertEqual(objects.a6.replace(objects.a7), objects.a6)
diff --git a/Testing/test_complex.py b/Testing/test_complex.py
index 3806486..d74f40b 100644
--- a/Testing/test_complex.py
+++ b/Testing/test_complex.py
@@ -1,80 +1,52 @@
 import collections
 import unittest
 
-from eBCSgen.Core.Structure import StructureAgent
-from eBCSgen.Core.Atomic import AtomicAgent
-from eBCSgen.Core.Complex import Complex, align_agents
-from eBCSgen.Parsing.ParseBCSL import Parser
+from eBCSgen.Core.Complex import align_agents
+import Testing.objects_testing as objects
 
 
 class TestComplex(unittest.TestCase):
-    def setUp(self):
-        self.a1 = AtomicAgent("T", "s")
-        self.a2 = AtomicAgent("S", "i")
-        self.a3 = AtomicAgent("U", "a")
-        self.a4 = AtomicAgent("T", "_")
-        self.a5 = AtomicAgent("U", "_")
-        self.a6 = AtomicAgent("S", "_")
-        self.a7 = AtomicAgent("U", "b")
-
-        self.s1 = StructureAgent("X", {self.a1})
-        self.s2 = StructureAgent("A", {self.a2, self.a3})
-        self.s3 = StructureAgent("X", {self.a4})
-        self.s4 = StructureAgent("A", {self.a2, self.a5})
-        self.s5 = StructureAgent("A", {self.a6, self.a3})
-        self.s6 = StructureAgent("A", set())
-        self.s7 = StructureAgent("A", {self.a2, self.a7})
-
-        self.c1 = Complex([self.s1, self.s2, self.s2], "cyt")
-        self.c2 = Complex([self.s3, self.s4, self.s5], "cyt")
-        self.c3 = Complex([self.s2, self.s2, self.s1], "cyt")
-        self.c4 = Complex([self.s2, self.s2, self.s1], "cell")
-
-        self.c5 = Complex([self.s2, self.s4, self.a1], "cell")
-        self.c6 = Complex([self.s6, self.s6, self.a4], "cell")
-
-        self.large_c1 = Complex([self.s4] * 6 + [self.s3] * 5, "cell")
-        self.large_c2 = Complex([self.s5] * 7 + [self.s3] * 6, "cell")
-
     def test_eq(self):
-        self.assertEqual(self.c1, self.c3)
-        self.assertNotEqual(self.c1, self.c4)
+        self.assertEqual(objects.c8, objects.c10)
+        self.assertNotEqual(objects.c8, objects.c11)
 
     def test_print(self):
-        self.assertEqual(str(self.c1), "X(T{s}).A(S{i},U{a}).A(S{i},U{a})::cyt")
-        self.assertEqual(str(self.c2), "X(T{_}).A(S{i},U{_}).A(S{_},U{a})::cyt")
+        self.assertEqual(str(objects.c8), "X(T{s}).A(S{i},U{a}).A(S{i},U{a})::cyt")
+        self.assertEqual(str(objects.c9), "X(T{_}).A(S{i},U{_}).A(S{_},U{a})::cyt")
 
     def test_compatibility(self):
-        self.assertTrue(self.c2.compatible(self.c1))
-        self.assertFalse(self.c1.compatible(self.c2))
-        self.assertFalse(self.c2.compatible(self.c4))
-        self.assertFalse(self.large_c1.compatible(self.large_c2))
+        self.assertTrue(objects.c9.compatible(objects.c8))
+        self.assertFalse(objects.c8.compatible(objects.c9))
+        self.assertFalse(objects.c9.compatible(objects.c11))
+        self.assertFalse(objects.large_c1.compatible(objects.large_c2))
 
     def test_to_PRISM_code(self):
-        self.assertEqual(self.c1.to_PRISM_code(5), "VAR_5")
+        self.assertEqual(objects.c8.to_PRISM_code(5), "VAR_5")
 
     def test_reduce_context(self):
-        self.assertEqual(self.c5.reduce_context(), self.c6)
+        self.assertEqual(objects.c12.reduce_context(), objects.c13)
 
     def test_create_all_compatible(self):
         atomic_signature = {"S": {"a", "i"}, "T": {"u", "p"}}
         structure_signature = {"KaiC": {"S"}}
 
-        complex_parser = Parser("rate_complex")
-
-        complex1 = complex_parser.parse("KaiC(S{a}).T{u}::cyt").data.children[0]
-        complex2 = complex_parser.parse("KaiC(S{a}).T{p}::cyt").data.children[0]
-        complex3 = complex_parser.parse("KaiC(S{i}).T{u}::cyt").data.children[0]
-        complex4 = complex_parser.parse("KaiC(S{i}).T{p}::cyt").data.children[0]
+        complex1 = objects.rate_complex_parser.parse("KaiC(S{a}).T{u}::cyt").data.children[0]
+        complex2 = objects.rate_complex_parser.parse("KaiC(S{a}).T{p}::cyt").data.children[0]
+        complex3 = objects.rate_complex_parser.parse("KaiC(S{i}).T{u}::cyt").data.children[0]
+        complex4 = objects.rate_complex_parser.parse("KaiC(S{i}).T{p}::cyt").data.children[0]
         results_1 = {complex1, complex2, complex3, complex4}
         results_2 = {complex1, complex2}
 
-        complex = complex_parser.parse("KaiC().T{_}::cyt").data.children[0]
-        output_comples = complex.create_all_compatible(atomic_signature, structure_signature)
+        complex = objects.rate_complex_parser.parse("KaiC().T{_}::cyt").data.children[0]
+        output_comples = complex.create_all_compatible(
+            atomic_signature, structure_signature
+        )
         self.assertEqual(output_comples, results_1)
 
-        complex = complex_parser.parse("KaiC(S{a}).T{_}::cyt").data.children[0]
-        output_comples = complex.create_all_compatible(atomic_signature, structure_signature)
+        complex = objects.rate_complex_parser.parse("KaiC(S{a}).T{_}::cyt").data.children[0]
+        output_comples = complex.create_all_compatible(
+            atomic_signature, structure_signature
+        )
         self.assertEqual(output_comples, results_2)
 
         atomic_signature = {"S": {"a", "i"}, "T": {"u", "p"}}
@@ -83,10 +55,15 @@ def test_create_all_compatible(self):
         results = set()
         with open("Testing/complexes_1.txt") as file:
             for complex in file.readlines():
-                results.add(complex_parser.parse(complex).data.children[0])
-
-        complex = complex_parser.parse("KaiC().KaiC().KaiC().KaiC().KaiC().KaiC()::cyt").data.children[0]
-        output_comples = complex.create_all_compatible(atomic_signature, structure_signature)
+                complex = complex.strip()
+                results.add(objects.rate_complex_parser.parse(complex).data.children[0])
+
+        complex = objects.rate_complex_parser.parse(
+            "KaiC().KaiC().KaiC().KaiC().KaiC().KaiC()::cyt"
+        ).data.children[0]
+        output_comples = complex.create_all_compatible(
+            atomic_signature, structure_signature
+        )
         self.assertEqual(output_comples, results)
 
         atomic_signature = {"A": {"+", "-"}, "B": {"HA", "HE"}}
@@ -95,22 +72,28 @@ def test_create_all_compatible(self):
         results = set()
         with open("Testing/complexes_2.txt") as file:
             for complex in file.readlines():
-                results.add(complex_parser.parse(complex).data.children[0])
-
-        complex = complex_parser.parse("KaiB().KaiB().KaiB()::cyt").data.children[0]
-        output_comples = complex.create_all_compatible(atomic_signature, structure_signature)
+                complex = complex.strip()
+                results.add(objects.rate_complex_parser.parse(complex).data.children[0])
+
+        complex = objects.rate_complex_parser.parse(
+            "KaiB().KaiB().KaiB()::cyt"
+        ).data.children[0]
+        output_comples = complex.create_all_compatible(
+            atomic_signature, structure_signature
+        )
         self.assertEqual(output_comples, results)
 
     def test_align_agents(self):
         agent = "X(T{s}).A(S{i},U{a}).A(S{i},U{b})::cyt"
         pattern = "A().A(S{i}).X()::cyt"
 
-        complex_parser = Parser("rate_complex")
-
-        complex = complex_parser.parse(agent).data.children[0]
-        pattern = complex_parser.parse(pattern).data.children[0]
+        complex = objects.rate_complex_parser.parse(agent).data.children[0]
+        pattern = objects.rate_complex_parser.parse(pattern).data.children[0]
 
-        expected_result = {(self.s2, self.s7, self.s1), (self.s7, self.s2, self.s1)}
+        expected_result = {
+            (objects.s9, objects.s14, objects.s8),
+            (objects.s14, objects.s9, objects.s8),
+        }
         result = align_agents(pattern.agents, collections.Counter(complex.agents))
         result = {tuple(item) for item in result}
         self.assertEqual(expected_result, result)
diff --git a/Testing/test_export_SBML.py b/Testing/test_export_SBML.py
index e8e19f1..c51a627 100644
--- a/Testing/test_export_SBML.py
+++ b/Testing/test_export_SBML.py
@@ -2,23 +2,22 @@
 import libsbml
 import os
 
-from eBCSgen.Parsing.ParseBCSL import Parser
+import Testing.objects_testing as objects
 
 
 class TestSBMLexport(unittest.TestCase):
     def setUp(self):
-        model_parser = Parser("model")
         self.models_to_test = {}
 
         if not os.path.exists("Testing/Output/"):
             os.mkdir("Testing/Output/")
 
         with open("Testing/BCSL_models_for_SBML_export/general.txt", "r") as model_exp:
-            self.models_to_test["general"] = model_parser.parse(model_exp.read()).data
+            self.models_to_test["general"] = objects.model_parser.parse(model_exp.read()).data
         with open("Testing/BCSL_models_for_SBML_export/izomorphic.txt", "r") as model_izo:
-            self.models_to_test["izomorphic"] = model_parser.parse(model_izo.read()).data
+            self.models_to_test["izomorphic"] = objects.model_parser.parse(model_izo.read()).data
         with open("Testing/BCSL_models_for_SBML_export/transition.txt", "r") as model_transition:
-            self.models_to_test["transition"] = model_parser.parse(model_transition.read()).data
+            self.models_to_test["transition"] = objects.model_parser.parse(model_transition.read()).data
 
     def test_by_validator(self):
         validator = libsbml.SBMLValidator()
diff --git a/Testing/test_formal_methods.py b/Testing/test_formal_methods.py
index 8a83a81..f67ac70 100644
--- a/Testing/test_formal_methods.py
+++ b/Testing/test_formal_methods.py
@@ -8,60 +8,18 @@
 from eBCSgen.Parsing.ParsePCTLformula import PCTLparser
 from eBCSgen.Parsing.ParseCTLformula import CTLparser
 
+from Testing.models.get_model_str import get_model_str
+
 
 class TestFormalMethods(unittest.TestCase):
     def setUp(self):
         self.parser = PCTLparser()
         self.model_parser = Parser("model")
-        self.model = """
-            #! rules
-            Y{i}::rep => Y{a}::rep @ p*[Y{i}::rep]
-            Y{i}::rep => Y{-}::rep @ (1-p)*[Y{i}::rep]
-            X()::rep + Y{a}::rep => X().Y{a}::rep @ q*[X()::rep]*[Y{a}::rep]
-            X(K{i}).Y{_}::rep => X(K{p}).Y{_}::rep @ p*[X(K{i}).Y{_}::rep] // also here
-
-            #! inits
-            2 X(K{i})::rep
-            1 Y{i}::rep
-
-            #! definitions
-            p = 0.3
-        """
-
-        self.tumor = """
-            #! rules
-            T(P{i})::x => T(P{m})::x @ a1*[T(P{i})::x]
-            T(P{m})::x => T(P{i})::x + T(P{i})::x @ a2*[T(P{m})::x]
-            T(P{i})::x => @ d1*[T(P{i})::x]
-            T(P{m})::x => @ d2*[T(P{m})::x]
-
-            #! inits
-            2 T(P{m})::x
-            1 T(P{i})::x
+        self.model = get_model_str("model3")
 
-            #! definitions
-            a1 = 1.2
-            a2 = 2
-            d1 = 0.8
-            d2 = 0.5
-        """
+        self.tumor = get_model_str("tumor_model")
 
-        self.tumor_parametric = """
-            #! rules
-            T(P{i})::x => T(P{m})::x @ a1*[T(P{i})::x]
-            T(P{m})::x => T(P{i})::x + T(P{i})::x @ a2*[T(P{m})::x]
-            T(P{i})::x => @ d1*[T(P{i})::x]
-            T(P{m})::x => @ d2*[T(P{m})::x]
-
-            #! inits
-            2 T(P{m})::x
-            1 T(P{i})::x
-
-            #! definitions
-            a1 = 1.2
-            a2 = 2
-            d1 = 0.8
-        """
+        self.tumor_parametric = get_model_str("tumor_parametric_model")
 
     def test_tumor_model_checking(self):
         model_parsed = self.model_parser.parse(self.tumor).data
@@ -109,19 +67,7 @@ def test_parametric_model(self):
         self.assertTrue("Result (initial states)" in str(result))
 
     def test_synthesis_out_of_scope(self):
-        model_str = """
-        #! rules
-        X()::rep => @ k1*[X()::rep]
-        Z()::rep => X()::rep @ k2
-        => Y()::rep @ 1/(1+([X()::rep])**4)
-
-        #! inits
-        2 X()::rep
-        Y()::rep
-
-        #! definitions
-        k2 = 5
-        """
+        model_str = get_model_str("model4")
         model = self.model_parser.parse(model_str).data
         formula = self.parser.parse('P <= 0.5[F X()::out = 1]')
         ts = model.generate_direct_transition_system()
@@ -131,19 +77,7 @@ def test_synthesis_out_of_scope(self):
         self.assertRaises(ComplexOutOfScope, PCTL.parameter_synthesis, ts, formula, '0<=k1<=1')
 
     def test_synthesis_simple(self):
-        model_str = """
-        #! rules
-        X()::rep => @ k1*[X()::rep]
-        Z()::rep => X()::rep @ k2
-        => Y()::rep @ 1/(1+([X()::rep])**4)
-
-        #! inits
-        2 X()::rep
-        Y()::rep
-
-        #! definitions
-        k2 = 5
-        """
+        model_str = get_model_str("model4")
         model = self.model_parser.parse(model_str).data
         ts = model.generate_direct_transition_system()
         ts.change_to_vector_backend()
@@ -159,20 +93,7 @@ def test_synthesis_simple(self):
         self.assertTrue("Result (initial states)" in str(output))
 
     def test_synthesis_advanced(self):
-        model_str = """
-        #! rules
-        // commenting
-        X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep] // also here
-        X(T{a})::rep => X(T{o})::rep @ k2*[Z()::rep]
-        => Y(P{f})::rep @ 1/(1+([X()::rep])**4) // ** means power (^)
-
-        #! inits
-        2 X(K{c}, T{e})::rep
-        Y(P{g}, N{l})::rep
-
-        #! definitions
-        k2 = 0.05 // also comment
-        """
+        model_str = get_model_str("model5")
         model = self.model_parser.parse(model_str).data
         ts = model.generate_direct_transition_system()
         ts.change_to_vector_backend()
@@ -188,20 +109,7 @@ def test_synthesis_advanced(self):
         self.assertTrue("Result (initial states)" in str(output))
 
     def test_model_checking_simple(self):
-        model_str = """
-        #! rules
-        X()::rep => @ k1*[X()::rep]
-        Z()::rep => X()::rep @ k2
-        => Y()::rep @ 1/(1+([X()::rep])**4)
-
-        #! inits
-        2 X()::rep
-        Y()::rep
-
-        #! definitions
-        k2 = 5
-        k1 = 2
-        """
+        model_str = get_model_str("model6")
         model = self.model_parser.parse(model_str).data
         ts = model.generate_direct_transition_system()
         ts.change_to_vector_backend()
@@ -217,21 +125,7 @@ def test_model_checking_simple(self):
         self.assertTrue("Result (for initial states)" in str(output))
 
     def test_model_checking_advanced(self):
-        model_str = """
-        #! rules
-        // commenting
-        X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep] // also here
-        X(T{a})::rep => X(T{o})::rep @ k2*[Z()::rep]
-        => Y(P{f})::rep @ 1/(1+([X()::rep])**4) // ** means power (^)
-
-        #! inits
-        2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-        Y(P{g}, N{l})::rep
-
-        #! definitions
-        k2 = 0.05 // also comment
-        k1 = 2
-        """
+        model_str = get_model_str("model7")
         model = self.model_parser.parse(model_str).data
         ts = model.generate_direct_transition_system()
         ts.change_to_vector_backend()
@@ -247,20 +141,7 @@ def test_model_checking_advanced(self):
         self.assertTrue("Result (for initial states)" in str(output))
 
     def test_CTL_model_checking(self):
-        model_str = """
-        #! rules
-        X()::rep => @ k1*[X()::rep]
-        Z()::rep => X()::rep @ k2
-        => Y()::rep @ 1/(1+([X()::rep])**4)
-
-        #! inits
-        2 X()::rep
-        Y()::rep
-
-        #! definitions
-        k2 = 5
-        k1 = 2
-        """
+        model_str = get_model_str("model6")
         model = self.model_parser.parse(model_str).data
         ts = model.generate_direct_transition_system()
         ts.change_to_vector_backend()
diff --git a/Testing/test_formal_methods_files.py b/Testing/test_formal_methods_files.py
index fa99ed9..4b5c200 100644
--- a/Testing/test_formal_methods_files.py
+++ b/Testing/test_formal_methods_files.py
@@ -5,10 +5,8 @@
 
 from eBCSgen.TS.Edge import Edge
 from eBCSgen.TS.TransitionSystem import TransitionSystem
-from eBCSgen.Core.Structure import StructureAgent
-from eBCSgen.Core.Complex import Complex
-from eBCSgen.Parsing.ParseBCSL import Parser
-from eBCSgen.TS.State import Vector, State, Memory
+
+import Testing.objects_testing as objects
 
 
 def get_storm_result(cmd: str):
@@ -24,40 +22,19 @@ def get_storm_result(cmd: str):
 
 class TestFormalMethods(unittest.TestCase):
     def setUp(self):
-        self.model_parser = Parser("model")
 
         """
         Model 1 - Transition system of die model
         Analysis of a PRISM example model from the Knuth-Yao
         source: storm website
         """
-        self.str1 = StructureAgent("S", set())
-        self.str2 = StructureAgent("D", set())
-
-        self.c1 = Complex([self.str1], "rep")
-        self.c2 = Complex([self.str2], "rep")
-
-        ordering = (self.c1, self.c2)
-
-        self.s1 = State(Vector(np.array((0, 0))), Memory(0))
-        self.s2 = State(Vector(np.array((1, 0))), Memory(0))
-        self.s3 = State(Vector(np.array((2, 0))), Memory(0))
-        self.s4 = State(Vector(np.array((3, 0))), Memory(0))
-        self.s5 = State(Vector(np.array((4, 0))), Memory(0))
-        self.s6 = State(Vector(np.array((5, 0))), Memory(0))
-        self.s7 = State(Vector(np.array((6, 0))), Memory(0))
-        self.s8 = State(Vector(np.array((7, 1))), Memory(0))
-        self.s9 = State(Vector(np.array((7, 2))), Memory(0))
-        self.s10 = State(Vector(np.array((7, 3))), Memory(0))
-        self.s11 = State(Vector(np.array((7, 4))), Memory(0))
-        self.s12 = State(Vector(np.array((7, 5))), Memory(0))
-        self.s13 = State(Vector(np.array((7, 6))), Memory(0))
+        ordering = (objects.c21, objects.c22)
 
         self.die_ts = TransitionSystem(ordering, 6)
         self.die_ts.init = 0
-        self.die_ts.states_encoding = {0: self.s1, 1: self.s2, 2: self.s3, 3: self.s4, 4: self.s5,
-                                       5: self.s6, 6: self.s7, 7: self.s8, 8: self.s9, 9: self.s10,
-                                       10: self.s11, 11: self.s12, 12: self.s13}
+        self.die_ts.states_encoding = {0: objects.state3, 1: objects.state4, 2: objects.state5, 3: objects.state6, 4: objects.state7,
+                                       5: objects.state8, 6: objects.state9, 7: objects.state10, 8: objects.state11, 9: objects.state12,
+                                       10: objects.state13, 11: objects.state14, 12: objects.state15}
         self.die_ts.edges = {Edge(0, 1, 0.5), Edge(0, 2, 0.5), Edge(1, 3, 0.5), Edge(1, 4, 0.5), Edge(2, 5, 0.5),
                              Edge(2, 6, 0.5), Edge(3, 1, 0.5), Edge(3, 7, 0.5), Edge(4, 8, 0.5), Edge(4, 9, 0.5),
                              Edge(5, 10, 0.5), Edge(5, 11, 0.5), Edge(6, 2, 0.5), Edge(6, 12, 0.5), Edge(7, 7, 1),
@@ -66,9 +43,9 @@ def setUp(self):
         # die parametric TS
         self.die_ts_parametric = TransitionSystem(ordering, 6)
         self.die_ts_parametric.init = 0
-        self.die_ts_parametric.states_encoding = {0: self.s1, 1: self.s2, 2: self.s3, 3: self.s4, 4: self.s5,
-                                                  5: self.s6, 6: self.s7, 7: self.s8, 8: self.s9, 9: self.s10,
-                                                  10: self.s11, 11: self.s12, 12: self.s13}
+        self.die_ts_parametric.states_encoding = {0: objects.state3, 1: objects.state4, 2: objects.state5, 3: objects.state6, 4: objects.state7,
+                                                  5: objects.state8, 6: objects.state9, 7: objects.state10, 8: objects.state11, 9: objects.state12,
+                                                  10: objects.state13, 11: objects.state14, 12: objects.state15}
         self.die_ts_parametric.edges = {Edge(0, 1, "p"), Edge(0, 2, "(1-p)"), Edge(1, 3, "p"), Edge(1, 4, "(1-p)"),
                                         Edge(2, 5, "p"),
                                         Edge(2, 6, "(1-p)"), Edge(3, 1, "p"), Edge(3, 7, "(1-p)"), Edge(4, 8, "p"),
diff --git a/Testing/test_model.py b/Testing/test_model.py
index 1148ece..b160d84 100644
--- a/Testing/test_model.py
+++ b/Testing/test_model.py
@@ -8,111 +8,41 @@
 from eBCSgen.Core.Formula import Formula
 from eBCSgen.Core.Model import Model
 from eBCSgen.Core.Rate import Rate
-from eBCSgen.Core.Structure import StructureAgent
-from eBCSgen.Core.Complex import Complex
-from eBCSgen.Core.Rule import Rule
-from eBCSgen.Parsing.ParseBCSL import Parser
 from eBCSgen.TS.State import Vector, State, Memory
 from eBCSgen.TS.VectorModel import VectorModel
 from eBCSgen.TS.VectorReaction import VectorReaction
 from eBCSgen.Core.Formula import AtomicProposition
 from eBCSgen.TS.TransitionSystem import TransitionSystem
 
+from Testing.models.get_model_str import get_model_str
+import Testing.objects_testing as objects
+
 
 class TestModel(unittest.TestCase):
     def setUp(self):
-        # agents
-
-        self.s1 = StructureAgent("X", set())
-        self.s2 = StructureAgent("Y", set())
-        self.s3 = StructureAgent("Z", set())
-
-        self.c1 = Complex([self.s1], "rep")
-        self.c2 = Complex([self.s2], "rep")
-        self.c3 = Complex([self.s3], "rep")
-
-        #  rules
-
-        sequence_1 = (self.s1,)
-        mid_1 = 1
-        compartments_1 = ["rep"]
-        complexes_1 = [(0, 0)]
-        pairs_1 = [(0, None)]
-        rate_1 = Rate("k1*[X()::rep]")
-
-        self.r1 = Rule(sequence_1, mid_1, compartments_1, complexes_1, pairs_1, rate_1)
-
-        sequence_2 = (self.s3, self.s1)
-        mid_2 = 1
-        compartments_2 = ["rep"] * 2
-        complexes_2 = [(0, 0), (1, 1)]
-        pairs_2 = [(0, 1)]
-
-        self.r2 = Rule(sequence_2, mid_2, compartments_2, complexes_2, pairs_2, None)
-
-        sequence_3 = (self.s2,)
-        mid_3 = 0
-        compartments_3 = ["rep"]
-        complexes_3 = [(0, 0)]
-        pairs_3 = [(None, 0)]
-        rate_3 = Rate("1.0/(1.0+([X()::rep])**4.0)")
-
-        self.r3 = Rule(sequence_3, mid_3, compartments_3, complexes_3, pairs_3, rate_3)
-
         # inits
-
-        self.inits = collections.Counter({self.c1: 2, self.c2: 1})
+        self.inits = collections.Counter({objects.c27: 2, objects.c28: 1})
 
         # defs
-
         self.defs = {'k1': 0.05, 'k2': 0.12}
 
-        self.model = Model({self.r1, self.r2, self.r3}, self.inits, self.defs, set())
+        self.model = Model({objects.r1, objects.r2, objects.r3}, self.inits, self.defs, set())
         # model
 
-        self.model_str_1 = """
-        #! rules
-        X()::rep => @ k1*[X()::rep]
-        Z()::rep => X()::rep
-        => Y()::rep @ 1/(1+([X()::rep])**4)
-
-        #! inits
-        2 X()::rep
-        Y()::rep
-        
-        #! definitions
-        k1 = 0.05
-        k2 = 0.12   
-        """
+        self.model_str_1 = get_model_str("model1")
 
-        self.model_parser = Parser("model")
-
-        self.model_str_2 = """
-        #! rules
-        X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep]
-        X(T{a})::rep => X(T{o})::rep @ k2*[Z()::rep]
-        => Y(P{f})::rep @ 1/(1+([X()::rep])**4)
-
-        #! inits
-        2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-        Y(P{g}, N{l})::rep
-
-        #! definitions
-        k1 = 0.05
-        k2 = 0.12
-        """
+        self.model_str_2 = get_model_str("model2")
 
         # vectors
 
-        ordering = (self.c1, self.c2, self.c3)
+        ordering = (objects.c27, objects.c28, objects.c29)
 
-        self.rate_parser = Parser("rate")
         rate_expr = "1/(1+([X()::rep])**4)"
-        rate_1 = Rate(self.rate_parser.parse(rate_expr).data)
+        rate_1 = Rate(objects.rate_parser.parse(rate_expr).data)
         rate_1.vectorize(ordering, dict())
 
         rate_expr = "k1*[X()::rep]"
-        rate_2 = Rate(self.rate_parser.parse(rate_expr).data)
+        rate_2 = Rate(objects.rate_parser.parse(rate_expr).data)
         rate_2.vectorize(ordering, {"k1": 0.05})
 
         rate_3 = Rate(Tree('rate', [Tree('fun', [1.0])]))
@@ -130,415 +60,97 @@ def setUp(self):
                                            rate_3)}
 
         self.vm_1 = VectorModel(vector_reactions, init, ordering, None)
-
-        # wrong models
-
-        self.model_wrong_1 = \
-            """#! rules
-            X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep]
-            X(T{a})::rep => X(T{o}):;rep @ k2*[Z()::rep]
-            => Y(P{f})::rep @ 1/(1+([X()::rep])**4)
-
-            #! inits
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-            Y(P{g}, N{l})::rep
-
-            #! definitions
-            k1 = 0.05
-            k2 = 0.12
-            """
-
-        self.model_wrong_2 = \
-            """#! rules
-            X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep]
-            X(T{a})::rep = X(T{o})::rep @ k2*[Z()::rep]
-            => Y(P{f})::rep @ 1/(1+([X()::rep])**4)
-
-            #! inits
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-            Y(P{g}, N{l})::rep
-
-            #! definitions
-            k1 = 0.05
-            k2 = 0.12
-            """
-
-        self.model_with_comments = """
-            #! rules
-            // commenting
-            X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep] // also here
-            X(T{a})::rep => X(T{o})::rep @ k2*[Z()::rep]
-            => Y(P{f})::rep @ 1/(1+([X()::rep])**4) // ** means power (^)
-
-            #! inits
-            // here
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-            Y(P{g}, N{l})::rep // comment just 1 item
-
-            #! definitions
-            // and
-            k1 = 0.05 // also
-            k2 = 0.12
-            """
         
-        self.model_with_complexes = """
-            #! rules
-            // commenting
-            X(T{a}):XX::rep => X(T{o}):XX::rep @ k2*[X().X()::rep]
-            K{i}:X():XYZ::rep => K{p}:X():XYZ::rep @ k1*[X().Y().Z()::rep] // also here
-            => P{f}:XP::rep @ 1/(1+([X().P{_}::rep])**4) // ** means power (^)
-
-            #! inits
-            // here
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-            Y(P{g}, N{l})::rep // comment just 1 item
-
-            #! definitions
-            // and
-            k1 = 0.05 // also
-            k2 = 0.12
-
-            #! complexes
-            XYZ = X().Y().Z() // a big complex
-            XX = X().X()
-            XP = X().P{_}
-            """
-
-        self.model_without_complexes = """
-            #! rules
-            // commenting
-            X(T{a}).X()::rep => X(T{o}).X()::rep @ k2*[X().X()::rep]
-            X(K{i}).Y().Z()::rep => X(K{p}).Y().Z()::rep @ k1*[X().Y().Z()::rep] // also here
-            => X().P{f}::rep @ 1/(1+([X().P{_}::rep])**4) // ** means power (^)
-
-            #! inits
-            // here
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-            Y(P{g}, N{l})::rep // comment just 1 item
-
-            #! definitions
-            // and
-            k1 = 0.05 // also
-            k2 = 0.12
-            """
-
-        self.model_with_variable = """
-            #! rules
-            // commenting
-            T{a}:X():?::rep => T{o}:X():?::rep @ k2*[X().X()::rep] ; ? = { XX, XY }
-            K{i}:X():XY::rep => K{p}:X():XY::rep @ k1*[X().Y().Z().X()::rep] // also here
-
-            #! inits
-            // here
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-
-            #! definitions
-            // and
-            k1 = 0.05 // also
-            k2 = 0.12
-
-            #! complexes
-            XX = X().X()
-            XY = X().Y()
-            """
-
-        self.model_without_variable = """
-            #! rules
-            // commenting
-            X(K{i}).Y()::rep => X(K{p}).Y()::rep @ k1*[X().Y().Z().X()::rep]
-            X(T{a}).X()::rep => X(T{o}).X()::rep @ k2*[X().X()::rep]
-            X(T{a}).Y()::rep => X(T{o}).Y()::rep @ k2*[X().X()::rep]
-
-            #! inits
-            // here
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-
-            #! definitions
-            // and
-            k1 = 0.05 // also
-            k2 = 0.12
-            """
-
-        self.model_with_redundant = """
-            #! rules
-            K(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt + D(A{_})::cell @ 3*[K().B()::cyt]/2*v_1
-            K().B()::cyt => K()::cyt + B()::cyt + D(A{_})::cell @ 3*[K().B()::cyt]/2*v_1
-            K().K()::cyt => K()::cyt + K()::cyt
-            K(S{i}).K()::cyt => K(S{a})::cyt + K()::cyt
-            K(S{i}, T{p}).K()::cyt => K(S{a}, T{p})::cyt + K()::cyt
-
-            #! inits
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-
-            #! definitions
-            v_1 = 0.05
-            k2 = 0.12
-            """
-
-        self.model_without_redundant = """
-            #! rules
-            K().B()::cyt => K()::cyt + B()::cyt + D(A{_})::cell @ 3*[K().B()::cyt]/2*v_1
-            K().K()::cyt => K()::cyt + K()::cyt
-
-            #! inits
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-
-            #! definitions
-            v_1 = 0.05
-            k2 = 0.12
-            """
-
-        self.model_with_context = """
-            #! rules
-            K(S{i}).B(T{a})::cyt => K(S{i})::cyt + B(T{a})::cyt @ 3*[K(S{i}).B(T{a})::cyt]/2*v_1
-            A{p}.K(S{i},T{i})::cyt => A{i}::cyt + K(S{a},T{a})::cyt
-            K(S{i},T{i})::cyt => K(S{a},T{i})::cyt
-
-            #! inits
-            2 K(S{i}).B(T{a})::cyt
-            1 A{p}.K(S{i},T{i})::cyt
-
-            #! definitions
-            v_1 = 0.05
-            k2 = 0.12
-            """
-
-        self.model_without_context = """
-            #! rules
-            K().B()::cyt => K()::cyt + B()::cyt @ 3*[K().B()::cyt]/2*v_1
-            A{_}.K()::cyt => A{_}::cyt + K()::cyt
-
-            #! inits
-            2 K().B()::cyt
-            1 A{_}.K()::cyt
-
-            #! definitions
-            v_1 = 0.05
-            k2 = 0.12
-            """
-
-        self.model_reachable = """
-            #! rules
-            K(S{i}).B()::cyt => K(S{a})::cyt + B()::cyt @ 3*[K(S{i}).B()::cyt]/2*v_1
-            K(S{a})::cyt + A{i}::cyt => K(S{a}).A{i}::cyt
-            K().A{i}::cyt => K().A{a}::cyt
-
-            #! inits
-            2 K(S{i}).B()::cyt
-            1 A{i}::cyt
-
-            #! definitions
-            v_1 = 0.05
-            k2 = 0.12
-            """
-
-        self.model_nonreachable = """
-            #! rules
-            K(S{i}).B()::cyt => K(S{a})::cyt + B()::cyt @ 3*[K(S{i}).B()::cyt]/2*v_1
-            K(S{a})::cyt + A{i}::cyt => K(S{a}).A{i}::cyt
-
-            #! inits
-            2 K(S{i}).B()::cyt
-            1 A{i}::cyt
-
-            #! definitions
-            v_1 = 0.05
-            k2 = 0.12
-            """
+        self.model_with_complexes = get_model_str("model_with_complexes")
+
+        self.model_without_complexes = get_model_str("model_without_complexes")
+
+        self.model_with_variable = get_model_str("model_with_variable")
+
+        self.model_without_variable = get_model_str("model_without_variable")
+
+        self.model_with_redundant = get_model_str("model_with_redundant")
+
+        self.model_without_redundant = get_model_str("model_without_redundant")
+
+        self.model_with_context = get_model_str("model_with_context")
+
+        self.model_without_context = get_model_str("model_without_context")
+
+        self.model_reachable = get_model_str("model_reachable")
+
+        self.model_nonreachable = get_model_str("model_nonreachable")
         
-        self.model_parametrised = """
-            #! rules
-            // commenting
-            X(K{i})::rep => X(K{p})::rep @ k1*[X()::rep] // also here
-            X(T{a})::rep => X(T{o})::rep @ k2*[Z()::rep]
-            => Y(P{f})::rep @ 1/(v_3+([X()::rep])**4) // ** means power (^)
-
-            #! inits
-            2 X(K{c}, T{e}).X(K{c}, T{j})::rep
-            Y(P{g}, N{l})::rep // comment just 1 item
-
-            #! definitions
-            k1 = 0.05
-            """
-
-        self.miyoshi = """
-            #! rules
-            S{u}:KaiC():KaiC6::cyt => S{p}:KaiC():KaiC6::cyt @ (kcat1*[KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
-            S{p}:KaiC():KaiC6::cyt => S{u}:KaiC():KaiC6::cyt @ (kcat2*[KaiB4{a}.KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
-            T{u}:KaiC():KaiC6::cyt => T{p}:KaiC():KaiC6::cyt @ (kcat3*[KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
-            T{p}:KaiC():KaiC6::cyt => T{u}:KaiC():KaiC6::cyt @ (kcat4*[KaiB4{a}.KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
-            KaiB4{i}::cyt => KaiB4{a}::cyt @ (kcatb2*[KaiB4{i}::cyt])/(Kmb2 + [KaiB4{i}::cyt])
-            KaiB4{a}::cyt => KaiB4{i}::cyt @ (kcatb1*[KaiB4{a}::cyt])/(Kmb1 + [KaiB4{a}::cyt])
-            KaiB4{a}.KaiA2()::cyt => KaiB4{a}::cyt + KaiA2()::cyt @ k12*[KaiB4{a}.KaiA2()::cyt]
-            KaiC6::cyt => 6 KaiC()::cyt @ kdimer*[KaiC6::cyt]
-            6 KaiC()::cyt => KaiC6::cyt @ kdimer*[KaiC()::cyt]*([KaiC()::cyt] - 1)*([KaiC()::cyt] - 2)*([KaiC()::cyt] - 3)*([KaiC()::cyt] - 4)*([KaiC()::cyt] - 5)
-
-            #! inits
-            6 KaiC(S{p},T{p})::cyt
-            1 KaiB4{a}.KaiA2()::cyt
-
-            #! definitions
-            kcat1 = 0.539
-            kcat3 = 0.89
-            Km = 0.602
-            kcatb2 = 0.346
-            kcatb1 = 0.602
-            Kmb2 = 66.75
-            Kmb1 = 2.423
-            k12 = 0.0008756
-            kdimer = 1.77
-
-            #! complexes
-            KaiC6 = KaiC().KaiC().KaiC().KaiC().KaiC().KaiC()
-            """
-
-        self.miyoshi_non_param = """
-            #! rules
-            S{u}:KaiC():KaiC6::cyt => S{p}:KaiC():KaiC6::cyt @ (kcat1*[KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
-            S{p}:KaiC():KaiC6::cyt => S{u}:KaiC():KaiC6::cyt @ (kcat2*[KaiB4{a}.KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
-            T{u}:KaiC():KaiC6::cyt => T{p}:KaiC():KaiC6::cyt @ (kcat3*[KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
-            T{p}:KaiC():KaiC6::cyt => T{u}:KaiC():KaiC6::cyt @ (kcat4*[KaiB4{a}.KaiA2()::cyt]*[KaiC6::cyt])/(Km + [KaiC6::cyt])
-            KaiB4{i}::cyt => KaiB4{a}::cyt @ (kcatb2*[KaiB4{i}::cyt])/(Kmb2 + [KaiB4{i}::cyt])
-            KaiB4{a}::cyt => KaiB4{i}::cyt @ (kcatb1*[KaiB4{a}::cyt])/(Kmb1 + [KaiB4{a}::cyt])
-            KaiB4{a}.KaiA2()::cyt => KaiB4{a}::cyt + KaiA2()::cyt @ k12*[KaiB4{a}.KaiA2()::cyt]
-            KaiC6::cyt => 6 KaiC()::cyt @ kdimer*[KaiC6::cyt]
-            6 KaiC()::cyt => KaiC6::cyt @ kdimer*[KaiC()::cyt]
-            KaiC().KaiC()::cyt => KaiC()::cyt + KaiC()::cyt @ kcat1
-
-            #! inits
-            6 KaiC(S{p},T{p})::cyt
-            1 KaiB4{a}.KaiA2()::cyt
-            1 KaiC(S{u},T{p})::cyt
-            1 KaiC(S{u},T{p}).KaiC(S{p},T{p})::cyt
-
-            #! definitions
-            kcat2 = 0.539
-            kcat4 = 0.89
-            kcat1 = 0.539
-            kcat3 = 0.89
-            Km = 0.602
-            kcatb2 = 0.346
-            kcatb1 = 0.602
-            Kmb2 = 66.75
-            Kmb1 = 2.423
-            k12 = 0.0008756
-            kdimer = 1.77
-
-            #! complexes
-            KaiC6 = KaiC().KaiC().KaiC().KaiC().KaiC().KaiC()
-            """
-
-        self.model_for_matching = """
-            #! rules
-            K(S{i}).B()::cyt + C{_}::cell => K(S{i})::cyt + B()::cyt + C{_}::cell @ 3*[K(S{i}).B()::cyt]/2*v_1
-            A{p}.K(S{i})::cyt => A{i}::cyt + K(S{a})::cyt + C{a}::cyt @ 0.3*[A{p}.K(S{i})::cyt]
-            K(S{i},T{i})::cyt + C{_}::cell => K(S{a},T{i})::cyt @ k2*[K(S{i},T{i})::cyt]
-            C{_}::cell + C{_}::cell => C{_}.C{_}::cell @ v_1*[C{_}::cell]**2
-            C{_}::cell + K()::cell => C{_}.K()::cell @ v_1*[C{_}::cell]**2
-    
-            #! inits
-            2 K(S{i},T{i}).B(T{a})::cyt
-            1 A{p}.K(S{i},T{i})::cyt
-            2 C{i}::cell
-            1 C{a}::cell
-    
-            #! definitions
-            v_1 = 0.05
-            k2 = 0.12
-            """
-
-        self.model_for_bound = """
-            #! rules
-             => A{i}::cyt @ k
-             => C{i}::cyt @ k
-             
-            #! inits
-            2 B{i}::cyt
-
-            #! definitions
-            k = 1
-            """
+        self.model_parametrised = get_model_str("model_parametrised")
 
-    def test_str(self):
-        model = self.model_parser.parse(self.model_str_1).data
-        back_to_str = repr(model)
-        parsed_again = self.model_parser.parse(back_to_str).data
-        self.assertEqual(model, parsed_again)
+        self.miyoshi = get_model_str("model_miyoshi")
 
-    def test_comments(self):
-        model_with_comments = self.model_parser.parse(self.model_with_comments)
-        model_without_comments = self.model_parser.parse(self.model_str_2).data
+        self.miyoshi_non_param = get_model_str("model_miyoshi_non_param")
 
-        self.assertEqual(model_with_comments.data, model_without_comments)
+        self.model_for_matching = get_model_str("model_for_matching")
 
-    def test_parser(self):
-        self.assertEqual(self.model_parser.parse(self.model_str_1).data, self.model)
+        self.model_for_bound = get_model_str("model_for_bound")
+
+    def test_str(self):
+        model = objects.model_parser.parse(self.model_str_1).data
+        back_to_str = repr(model)
+        parsed_again = objects.model_parser.parse(back_to_str).data
+        self.assertEqual(model, parsed_again)
 
     def test_signatures(self):
-        model = self.model_parser.parse(self.model_str_2).data
+        model = objects.model_parser.parse(self.model_str_2).data
         self.assertEqual(model.atomic_signature, {'K': {'c', 'i', 'p'}, 'T': {'e', 'a', 'o', 'j'},
                                                   'P': {'g', 'f'}, 'N': {'l'}})
         self.assertEqual(model.structure_signature, {'X': {'K', 'T'}, 'Y': {'P', 'N'}})
 
     def test_to_vector_model(self):
-        model = self.model_parser.parse(self.model_str_1).data
+        model = objects.model_parser.parse(self.model_str_1).data
         self.assertTrue(model.to_vector_model() == self.vm_1)
 
-    def test_parser_errors(self):
-        self.assertEqual(self.model_parser.parse(self.model_wrong_1).data,
-                         {"unexpected": ";", "expected": {'?', 'name'}, "line": 3, "column": 37})
-
-        self.assertEqual(self.model_parser.parse(self.model_wrong_2).data,
-                         {"expected": {'decimal', '#! inits', ']', '#! definitions', '=>, <=>', '@', 'int',
-                                       '+', 'name', ';', '}', ',', '#! complexes', '#! regulation'},
-                          "line": 3, "column": 26, "unexpected": "="})
-
     def test_zooming_syntax(self):
-        model_abstract = self.model_parser.parse(self.model_with_complexes).data
-        model_base = self.model_parser.parse(self.model_without_complexes).data
+        model_abstract = objects.model_parser.parse(self.model_with_complexes).data
+        model_base = objects.model_parser.parse(self.model_without_complexes).data
         self.assertEqual(model_abstract, model_base)
 
     def test_variables(self):
-        model_abstract = self.model_parser.parse(self.model_with_variable).data
-        model_base = self.model_parser.parse(self.model_without_variable).data
+        model_abstract = objects.model_parser.parse(self.model_with_variable).data
+        model_base = objects.model_parser.parse(self.model_without_variable).data
         self.assertEqual(model_abstract, model_base)
 
     def test_redundant(self):
-        model = self.model_parser.parse(self.model_with_redundant).data
+        model = objects.model_parser.parse(self.model_with_redundant).data
         model.eliminate_redundant()
 
-        model_eliminated = self.model_parser.parse(repr(model)).data
-        model_check = self.model_parser.parse(self.model_without_redundant).data
+        model_eliminated = objects.model_parser.parse(repr(model)).data
+        model_check = objects.model_parser.parse(self.model_without_redundant).data
         self.assertEqual(model_eliminated, model_check)
 
     def test_reduce_context(self):
-        model = self.model_parser.parse(self.model_with_context).data
+        model = objects.model_parser.parse(self.model_with_context).data
         model.reduce_context()
 
-        model_check = self.model_parser.parse(self.model_without_context).data
+        model_check = objects.model_parser.parse(self.model_without_context).data
         self.assertEqual(model, model_check)
 
     def test_nonreachability(self):
-        complex_parser = Parser("rate_complex")
         agent = "K(S{a}).A{a}::cyt"
-        complex = complex_parser.parse(agent).data.children[0]
+        complex = objects.rate_complex_parser.parse(agent).data.children[0]
 
-        model_reach = self.model_parser.parse(self.model_reachable).data
-        model_nonreach = self.model_parser.parse(self.model_nonreachable).data
+        model_reach = objects.model_parser.parse(self.model_reachable).data
+        model_nonreach = objects.model_parser.parse(self.model_nonreachable).data
 
         self.assertTrue(model_reach.static_non_reachability(complex))
         self.assertFalse(model_nonreach.static_non_reachability(complex))
 
     def test_parametrised_model(self):
-        model = self.model_parser.parse(self.model_parametrised).data
+        model = objects.model_parser.parse(self.model_parametrised).data
         self.assertTrue(len(model.params) == 2)
 
     def test_create_complex_labels(self):
-        complex_parser = Parser("rate_complex")
-        complex_1 = complex_parser.parse("K(S{i},T{a}).B{o}::cyt").data.children[0]
-        complex_2 = complex_parser.parse("K(S{a},T{a}).B{o}::cyt").data.children[0]
-        complex_3 = complex_parser.parse("K(S{a},T{i}).B{o}::cyt").data.children[0]
-        complex_abstract = complex_parser.parse("K(S{a}).B{_}::cyt").data.children[0]
+        complex_1 = objects.rate_complex_parser.parse("K(S{i},T{a}).B{o}::cyt").data.children[0]
+        complex_2 = objects.rate_complex_parser.parse("K(S{a},T{a}).B{o}::cyt").data.children[0]
+        complex_3 = objects.rate_complex_parser.parse("K(S{a},T{i}).B{o}::cyt").data.children[0]
+        complex_abstract = objects.rate_complex_parser.parse("K(S{a}).B{_}::cyt").data.children[0]
 
         ordering = (complex_1, complex_2, complex_3)
         complexes = [complex_2, complex_abstract, complex_1]
@@ -554,12 +166,11 @@ def test_create_complex_labels(self):
         self.assertEqual(prism_formulas, result_formulas)
 
     def test_create_AP_labels(self):
-        complex_parser = Parser("rate_complex")
-        complex_1 = complex_parser.parse("K(S{i},T{a}).B{o}::cyt").data.children[0]
-        complex_2 = complex_parser.parse("K(S{a},T{a}).B{o}::cyt").data.children[0]
-        complex_3 = complex_parser.parse("K(S{a},T{i}).B{o}::cyt").data.children[0]
+        complex_1 = objects.rate_complex_parser.parse("K(S{i},T{a}).B{o}::cyt").data.children[0]
+        complex_2 = objects.rate_complex_parser.parse("K(S{a},T{a}).B{o}::cyt").data.children[0]
+        complex_3 = objects.rate_complex_parser.parse("K(S{a},T{i}).B{o}::cyt").data.children[0]
 
-        complex_abstract = complex_parser.parse("K(S{a}).B{_}::cyt").data.children[0]
+        complex_abstract = objects.rate_complex_parser.parse("K(S{a}).B{_}::cyt").data.children[0]
 
         ordering = (complex_1, complex_2, complex_3)
 
@@ -586,7 +197,7 @@ def test_create_AP_labels(self):
         self.assertEqual(AP_lables, result_AP_lables)
 
     def test_create_unique_agents(self):
-        model = self.model_parser.parse(self.miyoshi).data
+        model = objects.model_parser.parse(self.miyoshi).data
 
         reactions = set()
         unique_complexes = set()
@@ -603,14 +214,14 @@ def test_create_unique_agents(self):
         self.assertEqual(unique_complexes, set(ordering))
 
     def test_compute_bound(self):
-        model = self.model_parser.parse(self.model_for_bound).data
+        model = objects.model_parser.parse(self.model_for_bound).data
         for rule in model.rules:
             rule.lhs, rule.rhs = rule.create_complexes()
 
         self.assertEqual(model.compute_bound(), 2)
 
     def test_direct_ts_bound(self):
-        model = self.model_parser.parse(self.model_for_bound).data
+        model = objects.model_parser.parse(self.model_for_bound).data
         rules = set()
         for i, rule in enumerate(model.rules):
             rule.label = 'r{}'.format(i)
diff --git a/Testing/test_rate.py b/Testing/test_rate.py
index eb10ed9..e961eaf 100644
--- a/Testing/test_rate.py
+++ b/Testing/test_rate.py
@@ -9,6 +9,8 @@
 from eBCSgen.TS.State import Vector, State, Memory
 from eBCSgen.Core.Rate import Rate
 
+import Testing.objects_testing as objects
+
 
 CORRECT_MATHML = """
 <kineticLaw><math xmlns="http://www.w3.org/1998/Math/MathML"><apply><divide><times><cn>3.0</cn>\
@@ -18,47 +20,20 @@
 
 class TestRate(unittest.TestCase):
     def setUp(self):
-        # agents
-
-        self.a1 = AtomicAgent("T", "i")
-        self.a2 = AtomicAgent("S", "i")
-        self.a3 = AtomicAgent("T", "a")
-        self.a4 = AtomicAgent("S", "a")
-
-        self.s1 = StructureAgent("X", {self.a4})
-        self.s2 = StructureAgent("X", {self.a2})
-        self.s3 = StructureAgent("K", {self.a3})
-        self.s4 = StructureAgent("K", {self.a1})
-        self.s5 = StructureAgent("X", set())
-        self.s6 = StructureAgent("K", set())
-
-        self.c1 = Complex([self.s6], "cyt")  # K()::cyt
-        self.c2 = Complex([self.s3], "cyt")  # K(T{a})::cyt
-        self.c3 = Complex([self.s4], "cyt")  # K(T{i})::cyt
-        self.c4 = Complex([self.s4, self.s1], "cyt")  # K(T{i}).X(S{a})::cyt
-        self.c5 = Complex([self.s4, self.s2], "cyt")  # K(T{i}).X(S{i})::cyt
-        self.c6 = Complex([self.s3, self.s1], "cyt")  # K(T{a}).X(S{a})::cyt
-        self.c7 = Complex([self.s3, self.s2], "cyt")  # K(T{a}).X(S{i})::cyt
-
         # rates
 
         self.parser = Parser("rate")
         rate_expr = "3.0*[K()::cyt]/2.0*v_1"
 
-        self.rate_1 = Rate(self.parser.parse(rate_expr).data)
+        self.rate_1 = Rate(objects.rate_parser.parse(rate_expr).data)
 
         rate_expr = "3.0*[K(T{i}).X()::cyt] + [K()::cyt]"
 
-        self.rate_2 = Rate(self.parser.parse(rate_expr).data)
+        self.rate_2 = Rate(objects.rate_parser.parse(rate_expr).data)
 
         rate_expr = "(3.0*[K()::cyt])/(2.0*v_1)"
 
-        self.rate_3 = Rate(self.parser.parse(rate_expr).data)
-
-        # states
-
-        self.state_1 = State(Vector(np.array([2, 3])), Memory(0))
-        self.state_2 = State(Vector(np.array([2, 0, 3, 1, 6, 2])), Memory(0))
+        self.rate_3 = Rate(objects.rate_parser.parse(rate_expr).data)
 
     def test_to_str(self):
         self.assertEqual(str(self.rate_1), "3.0*[K()::cyt]/2.0*v_1")
@@ -68,46 +43,46 @@ def test_eq(self):
         self.assertNotEqual(self.rate_1, self.rate_2)
 
     def test_vectorize(self):
-        ordering = (self.c2, self.c3)
+        ordering = (objects.c15, objects.c16)
         self.assertEqual(self.rate_1.vectorize(ordering, dict()), [Vector(np.array([1, 1]))])
-        ordering = (self.c2, self.c3, self.c4, self.c5, self.c6, self.c7)
+        ordering = (objects.c15, objects.c16, objects.c17, objects.c18, objects.c19, objects.c20)
         self.assertEqual(self.rate_2.vectorize(ordering, dict()),
                          [Vector(np.array([0, 0, 1, 1, 0, 0])), Vector(np.array([1, 1, 0, 0, 0, 0]))])
 
     def test_evaluate(self):
-        ordering = (self.c2, self.c3)
+        ordering = (objects.c15, objects.c16)
         self.rate_1.vectorize(ordering, {"v_1": 5})
-        self.assertEqual(self.rate_1.evaluate(self.state_1), sympy.sympify("3*5.0/2*5"))
+        self.assertEqual(self.rate_1.evaluate(objects.state1), sympy.sympify("3*5.0/2*5"))
 
-        ordering = (self.c2, self.c3, self.c4, self.c5, self.c6, self.c7)
+        ordering = (objects.c15, objects.c16, objects.c17, objects.c18, objects.c19, objects.c20)
         self.rate_2.vectorize(ordering, dict())
-        self.assertEqual(self.rate_2.evaluate(self.state_2), sympy.sympify("3*4.0 + 2"))
+        self.assertEqual(self.rate_2.evaluate(objects.state2), sympy.sympify("3*4.0 + 2"))
 
     def test_to_symbolic(self):
-        ordering = (self.c2, self.c3)
+        ordering = (objects.c15, objects.c16)
         self.rate_1.vectorize(ordering, dict())
         self.rate_1.to_symbolic()
         self.assertEqual(str(self.rate_1), "3.0*(y[0] + y[1])/2.0*v_1")
 
-        ordering = (self.c2, self.c3, self.c4, self.c5, self.c6, self.c7)
+        ordering = (objects.c15, objects.c16, objects.c17, objects.c18, objects.c19, objects.c20)
         self.rate_2.vectorize(ordering, dict())
         self.rate_2.to_symbolic()
         self.assertEqual(str(self.rate_2), "3.0*(y[2] + y[3])+(y[0] + y[1])")
 
     def test_reduce_context(self):
         rate_expr = "3.0*[K(S{i})::cyt]/2.0*v_1"
-        rate = Rate(self.parser.parse(rate_expr).data)
+        rate = Rate(objects.rate_parser.parse(rate_expr).data)
 
         self.assertEqual(rate.reduce_context(), self.rate_1)
 
     def test_evaluate_direct(self):
-        values = {self.c1: 3}
+        values = {objects.c14: 3}
         params = {"v_1": 5}
         self.assertEqual(self.rate_3.evaluate_direct(values, params), 9/10)
 
     def test_mathML(self):
         rate_expr = "(3.0*[K(T{i}).X()::cyt])/([K()::cyt]**2.0+4.0*p)"
-        rate = Rate(self.parser.parse(rate_expr).data)
+        rate = Rate(objects.rate_parser.parse(rate_expr).data)
         expression = rate.to_mathML()
         agents, params = rate.get_params_and_agents()
         str_to_code = {"[" + str(agent) + "]": agent.to_SBML_species_code() for agent in agents}
diff --git a/Testing/test_regulations.py b/Testing/test_regulations.py
index b697a22..d34bf54 100644
--- a/Testing/test_regulations.py
+++ b/Testing/test_regulations.py
@@ -1,6 +1,7 @@
 import unittest
 
 from eBCSgen.Parsing.ParseBCSL import Parser
+from Testing.models.get_model_str import get_model_str
 
 
 class TestRegulations(unittest.TestCase):
@@ -8,29 +9,10 @@ def setUp(self):
         self.model_parser = Parser("model")
         self.complex_parser = Parser("rate_complex")
 
-        self.model_with_labels = """
-            #! rules
-            r1_S ~ A(S{i})::cell => A(S{a})::cell @ k1*[A(S{i})::cell]
-            r1_T ~ A(T{i})::cell => A(T{a})::cell @ k2*[A(T{i})::cell]
-            r2 ~ A()::cell => A()::out @ k3*[A()::cell]
-
-            #! inits
-            1 A(S{i},T{i})::cell
-
-            #! definitions
-            k1 = 0.3
-            k2 = 0.5
-            k3 = 0.1
-            """
+        self.model_with_labels = get_model_str("model_with_labels")
 
     def test_programmed(self):
-        regulation = """
-
-        #! regulation
-        type programmed
-        r1_S: {r1_T, r2}
-        r1_T: {r1_S}
-        """
+        regulation = get_model_str("regulation1")
         model = self.model_parser.parse(self.model_with_labels + regulation).data
 
         direct_ts = model.generate_direct_transition_system()
@@ -42,12 +24,7 @@ def test_programmed(self):
         self.assertEqual(direct_ts, indirect_ts)
 
     def test_ordered(self):
-        regulation = """
-
-        #! regulation
-        type ordered
-        (r1_S, r2), (r1_T, r2)
-        """
+        regulation = get_model_str("regulation2")
 
         model = self.model_parser.parse(self.model_with_labels + regulation).data
 
@@ -60,12 +37,7 @@ def test_ordered(self):
         self.assertEqual(direct_ts, indirect_ts)
 
     def test_conditional(self):
-        regulation = """
-
-        #! regulation
-        type conditional
-        r2: {A(S{a},T{i})::cell}
-        """
+        regulation = get_model_str("regulation3")
 
         model = self.model_parser.parse(self.model_with_labels + regulation).data
 
@@ -78,12 +50,7 @@ def test_conditional(self):
         self.assertEqual(direct_ts, indirect_ts)
 
     def test_concurrent_free(self):
-        regulation = """
-
-        #! regulation
-        type concurrent-free
-        (r1_S, r2), (r1_T, r2)
-        """
+        regulation = get_model_str("regulation4")
 
         model = self.model_parser.parse(self.model_with_labels + regulation).data
 
@@ -96,12 +63,7 @@ def test_concurrent_free(self):
         self.assertEqual(direct_ts, indirect_ts)
 
     def test_regular(self):
-        regulation = """
-
-        #! regulation
-        type regular
-        (r1_Sr1_Tr2|r1_Tr1_Sr2)
-        """
+        regulation = get_model_str("regulation5")
 
         model = self.model_parser.parse(self.model_with_labels + regulation).data
 
@@ -125,13 +87,7 @@ def test_no_regulation(self):
         self.assertEqual(direct_ts, indirect_ts)
 
     def test_network_free_simulation_regulated(self):
-        regulation = """
-
-        #! regulation
-        type programmed
-        r1_S: {r1_T, r2}
-        r1_T: {r1_S}
-        """
+        regulation = get_model_str("regulation1")
         model = self.model_parser.parse(self.model_with_labels + regulation).data
 
         result = model.network_free_simulation(5)
diff --git a/Testing/test_rule.py b/Testing/test_rule.py
index eea15da..4afa112 100644
--- a/Testing/test_rule.py
+++ b/Testing/test_rule.py
@@ -1,281 +1,104 @@
 import unittest
 
-from eBCSgen.Core.Rate import Rate
-from eBCSgen.Core.Structure import StructureAgent
-from eBCSgen.Core.Atomic import AtomicAgent
-from eBCSgen.Core.Complex import Complex
 from eBCSgen.Core.Rule import Rule
-from eBCSgen.Core.Side import Side
-from eBCSgen.Core.Reaction import Reaction
-from eBCSgen.Parsing.ParseBCSL import Parser
 
+import Testing.objects_testing as objects
 
-class TestRule(unittest.TestCase):
-    def setUp(self):
-        self.a1 = AtomicAgent("S", "u")
-        self.a2 = AtomicAgent("S", "p")
-        self.a3 = AtomicAgent("B", "_")
-        self.a4 = AtomicAgent("B", "-")
-        self.a5 = AtomicAgent("B", "+")
-
-        self.s1 = StructureAgent("K", {self.a1})
-        self.s2 = StructureAgent("B", set())
-        self.s3 = StructureAgent("K", {self.a2})
-        self.s4 = StructureAgent("B", set())
-        self.s5 = StructureAgent("D", {self.a3})
-        self.s6 = StructureAgent("K", {self.a4})
-        self.s7 = StructureAgent("K", {self.a5})
-
-        self.c1 = Complex([self.s1, self.s2], "cyt")
-        self.c2 = Complex([self.s3], "cyt")
-        self.c3 = Complex([self.s2], "cyt")
-        self.c4 = Complex([self.s5], "cell")
-
-        #  rules
-
-        sequence_1 = (self.s1, self.s2, self.s3, self.s4)
-        mid_1 = 2
-        compartments_1 = ["cyt"] * 4
-        complexes_1 = [(0, 1), (2, 2), (3, 3)]
-        pairs_1 = [(0, 2), (1, 3)]
-        rate_1 = Rate("3.0*[K()::cyt]/2.0*v_1")
-
-        self.r1 = Rule(sequence_1, mid_1, compartments_1, complexes_1, pairs_1, rate_1)
-
-        sequence_2 = (self.s1, self.s2, self.s3, self.s4, self.s5)
-        mid_2 = 2
-        compartments_2 = ["cyt"] * 4 + ["cell"]
-        complexes_2 = [(0, 1), (2, 2), (3, 3), (4, 4)]
-        pairs_2 = [(0, 2), (1, 3), (None, 4)]
-        rate_2 = Rate("3.0*[K()::cyt]/2.0*v_1")
-
-        self.r2 = Rule(sequence_2, mid_2, compartments_2, complexes_2, pairs_2, rate_2)
-
-        sequence_3 = (self.s6, self.s2, self.s5, self.s7, self.s4)
-        mid_3 = 3
-        compartments_3 = ["cyt"] * 2 + ["cell"] + ["cyt"] * 2
-        complexes_3 = [(0, 1), (2, 2), (3, 3), (4, 4)]
-        pairs_3 = [(0, 3), (1, 4), (2, None)]
-        rate_3 = Rate("3.0*[K(T{3+})::cyt]/2.0*v_1")
-
-        self.r3 = Rule(sequence_3, mid_3, compartments_3, complexes_3, pairs_3, rate_3)
-
-        # special cases
-
-        self.s1_s = StructureAgent("X", set())
-        self.s2_s = StructureAgent("Y", set())
-        self.s3_s = StructureAgent("Z", set())
-
-        sequence_4 = (self.s1_s, )
-        mid_4 = 1
-        compartments_4 = ["rep"]
-        complexes_4 = [(0, 0)]
-        pairs_4 = [(0, None)]
-        rate_4 = Rate("k1*[X()::rep]")
-
-        self.r4 = Rule(sequence_4, mid_4, compartments_4, complexes_4, pairs_4, rate_4)
-
-        sequence_5 = (self.s2_s, )
-        mid_5 = 0
-        compartments_5 = ["rep"]
-        complexes_5 = [(0, 0)]
-        pairs_5 = [(None, 0)]
-        rate_5 = Rate("1.0/(1.0+([X()::rep])**4.0)")
-
-        self.r5 = Rule(sequence_5, mid_5, compartments_5, complexes_5, pairs_5, rate_5)
-
-        #  reactions
-
-        lhs = Side([self.c1])
-        rhs = Side([self.c2, self.c3, self.c4])
-
-        self.reaction1 = Reaction(lhs, rhs, rate_2)
-
-        #  create
-
-        self.t_i = AtomicAgent("T", "i")
-        self.t_a = AtomicAgent("T", "a")
-
-        self.a4_p = AtomicAgent("C", "p")
-        self.a4_u = AtomicAgent("C", "u")
-
-        self.u2_c1_p = AtomicAgent("U", "p")
-        self.u2_c1_u = AtomicAgent("U", "u")
-
-        self.s6 = StructureAgent("D", set())
-        self.s6_c1_p = StructureAgent("D", {self.a4_p})
-        self.s6_c1_u = StructureAgent("D", {self.a4_u})
-
-        self.s2_c1_p = StructureAgent("B", {self.u2_c1_p})
-        self.s2_c1_u = StructureAgent("B", {self.u2_c1_u})
-
-        self.s1_c1_a = StructureAgent("K", {self.a1, self.t_a})
-        self.s1_c1_i = StructureAgent("K", {self.a1, self.t_i})
-
-        self.s3_c1_a = StructureAgent("K", {self.a2, self.t_a})
-        self.s3_c1_i = StructureAgent("K", {self.a2, self.t_i})
-
-        sequence_c1 = (self.s1, self.s2, self.s3, self.s4, self.s6)
-        mid_c1 = 2
-        compartments_c1 = ["cyt"] * 5
-        complexes_c1 = [(0, 0), (1, 1), (2, 3), (4, 4)]
-        pairs_c1 = [(0, 2), (1, 3), (None, 4)]
-        rate_c1 = Rate("3*[K()::cyt]/2*v_1")
-
-        self.c1_c1 = Complex([self.s2_c1_u], "cyt")  # B(U{u})::cyt
-        self.c1_c2 = Complex([self.s2_c1_p], "cyt")  # B(U{p})::cyt
-        self.c1_c3 = Complex([self.s1_c1_a], "cyt")  # K(S{u},T{a})::cyt
-        self.c1_c4 = Complex([self.s1_c1_i], "cyt")  # K(S{u},T{i})::cyt
-        self.c1_c5 = Complex([self.s3_c1_a, self.s2_c1_u], "cyt")  # K(S{p},T{a}).B(U{u})::c
-        self.c1_c6 = Complex([self.s3_c1_i, self.s2_c1_u], "cyt")  # K(S{p},T{i}).B(U{u})::c
-        self.c1_c7 = Complex([self.s3_c1_i, self.s2_c1_p], "cyt")  # K(S{p},T{i}).B(U{p})::c
-        self.c1_c8 = Complex([self.s3_c1_a, self.s2_c1_p], "cyt")  # K(S{p},T{a}).B(U{p})::c
-        self.c1_c9 = Complex([self.s6_c1_p], "cyt")  # D(C{p})::cyt
-        self.c1_c10 = Complex([self.s6_c1_u], "cyt")  # D(C{u})::cyt
-
-        self.rule_c1 = Rule(sequence_c1, mid_c1, compartments_c1, complexes_c1, pairs_c1, rate_c1)
-
-        self.reaction_c1_1 = Reaction(Side([self.c1_c1, self.c1_c3]),
-                                      Side([self.c1_c5, self.c1_c9]), rate_c1)
-        self.reaction_c1_2 = Reaction(Side([self.c1_c1, self.c1_c3]),
-                                      Side([self.c1_c5, self.c1_c10]), rate_c1)
-        self.reaction_c1_3 = Reaction(Side([self.c1_c2, self.c1_c4]),
-                                      Side([self.c1_c7, self.c1_c10]), rate_c1)
-        self.reaction_c1_4 = Reaction(Side([self.c1_c1, self.c1_c4]),
-                                      Side([self.c1_c6, self.c1_c9]), rate_c1)
-        self.reaction_c1_5 = Reaction(Side([self.c1_c2, self.c1_c3]),
-                                      Side([self.c1_c8, self.c1_c9]), rate_c1)
-        self.reaction_c1_6 = Reaction(Side([self.c1_c2, self.c1_c3]),
-                                      Side([self.c1_c8, self.c1_c10]), rate_c1)
-        self.reaction_c1_7 = Reaction(Side([self.c1_c1, self.c1_c4]),
-                                      Side([self.c1_c6, self.c1_c10]), rate_c1)
-        self.reaction_c1_8 = Reaction(Side([self.c1_c2, self.c1_c4]),
-                                      Side([self.c1_c7, self.c1_c9]), rate_c1)
-
-        self.reactions_c1 = {self.reaction_c1_1, self.reaction_c1_2, self.reaction_c1_3, self.reaction_c1_4,
-                             self.reaction_c1_5, self.reaction_c1_6, self.reaction_c1_7, self.reaction_c1_8}
-
-        # context no change
-
-        sequence_no_change = (self.s1_c1_a, self.s2_c1_u, self.s3_c1_a, self.s2_c1_u, self.s6_c1_p)
-        self.rule_no_change = Rule(sequence_no_change, mid_c1, compartments_c1, complexes_c1, pairs_c1, rate_c1)
-
-        # parsing
-
-        self.parser = Parser("rule")
-
-        self.rule_no_rate = Rule(sequence_1, mid_1, compartments_1, complexes_1, pairs_1, None)
 
+class TestRule(unittest.TestCase):
+    
     def test_eq(self):
-        self.assertEqual(self.r1, self.r1)
+        self.assertEqual(objects.r4, objects.r4)
 
     def test_print(self):
-        self.assertEqual(str(self.r1), "K(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt @ 3.0*[K()::cyt]/2.0*v_1")
-        self.assertEqual(str(self.r2),
+        self.assertEqual(str(objects.r4), "K(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt @ 3.0*[K()::cyt]/2.0*v_1")
+        self.assertEqual(str(objects.r5),
                          "K(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt + D(B{_})::cell @ 3.0*[K()::cyt]/2.0*v_1")
 
     def test_create_complexes(self):
-        lhs = Side([self.c1])
-        rhs = Side([self.c2, self.c3, self.c4])
-        self.assertEqual(self.r2.create_complexes(), (lhs, rhs))
+        self.assertEqual(objects.r5.create_complexes(), (objects.lhs, objects.rhs))
 
     def test_to_reaction(self):
-        self.assertEqual(self.r2.to_reaction(), self.reaction1)
+        self.assertEqual(objects.r5.to_reaction(), objects.reaction1)
 
     def test_create_reactions(self):
         atomic_signature = {"T": {"a", "i"}, "U": {"p", "u"}, "C": {"p", "u"}, "S": {"p", "u"}}
         structure_signature = {"K": {"T", "S"}, "B": {"U"}, "D": {"C"}}
-        self.assertEqual(self.rule_c1.create_reactions(atomic_signature, structure_signature),
-                         self.reactions_c1)
+        self.assertEqual(objects.rule_c1.create_reactions(atomic_signature, structure_signature),
+                         objects.reactions_c1)
 
-        self.assertEqual(self.rule_no_change.create_reactions(atomic_signature, structure_signature),
-                         {self.reaction_c1_1})
+        self.assertEqual(objects.rule_no_change.create_reactions(atomic_signature, structure_signature),
+                         {objects.reaction_c1_1})
 
         rule_exp = "K(T{a}).K().K()::cyt => K(T{i}).K().K()::cyt @ k1*[K(T{a}).K().K()::cyt]"
-        rule = self.parser.parse(rule_exp).data[1]
+        rule = objects.rule_parser.parse(rule_exp).data[1]
         result = rule.create_reactions(atomic_signature, structure_signature)
 
         reactions = set()
         with open("Testing/reactions.txt") as file:
             for complex in file.readlines():
-                rule = self.parser.parse(complex).data[1]
+                complex = complex.strip()
+                rule = objects.rule_parser.parse(complex).data[1]
                 reactions.add(rule.to_reaction())
 
         self.assertEqual(result, reactions)
 
-    def test_parser(self):
-        rule_expr = "K(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt + D(B{_})::cell @ 3*[K()::cyt]/2*v_1"
-        self.assertEqual(self.parser.parse(rule_expr).data[1], self.r2)
-
-        rule_expr = "K(B{-}).B()::cyt + D(B{_})::cell => K(B{+})::cyt + B()::cyt @ 3*[K(T{3+})::cyt]/2*v_1"
-        self.assertEqual(self.parser.parse(rule_expr).data[1], self.r3)
-
-        rule_expr = "X()::rep => @ k1*[X()::rep]"
-        self.assertEqual(self.parser.parse(rule_expr).data[1], self.r4)
-
-        rule_expr = "=> Y()::rep @ 1/(1+([X()::rep])**4)"
-        self.assertEqual(self.parser.parse(rule_expr).data[1], self.r5)
-
-        rule_expr = "K(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt"
-        self.assertEqual(self.parser.parse(rule_expr).data[1], self.rule_no_rate)
-
     def test_compatible(self):
-        self.assertTrue(self.r1.compatible(self.r2))
-        self.assertFalse(self.r2.compatible(self.r1))
+        self.assertTrue(objects.r4.compatible(objects.r5))
+        self.assertFalse(objects.r5.compatible(objects.r4))
 
         rule_expr_1 = "K(S{u}).B()::cyt => K(S{p})::cyt + B()::cyt + D(B{_})::cell @ 3*[K()::cyt]/2*v_1"
         rule_expr_2 = "K().B()::cyt => K()::cyt + B()::cyt + D(B{_})::cell @ 3*[K()::cyt]/2*v_1"
-        rule1 = self.parser.parse(rule_expr_1).data[1]
-        rule2 = self.parser.parse(rule_expr_2).data[1]
+        rule1 = objects.rule_parser.parse(rule_expr_1).data[1]
+        rule2 = objects.rule_parser.parse(rule_expr_2).data[1]
 
         self.assertFalse(rule1.compatible(rule2))
         self.assertTrue(rule2.compatible(rule1))
 
     def test_reduce_context(self):
         rule_expr_1 = "K(S{u}).B{i}::cyt => K(S{p})::cyt + B{a}::cyt + D(B{_})::cell @ 3*[K(S{u}).B{i}::cyt]/2*v_1"
-        rule1 = self.parser.parse(rule_expr_1).data[1]
+        rule1 = objects.rule_parser.parse(rule_expr_1).data[1]
 
         rule_expr_2 = "K().B{_}::cyt => K()::cyt + B{_}::cyt + D()::cell @ 3*[K().B{_}::cyt]/2*v_1"
-        rule2 = self.parser.parse(rule_expr_2).data[1]
+        rule2 = objects.rule_parser.parse(rule_expr_2).data[1]
 
         self.assertEqual(rule1.reduce_context(), rule2)
 
         # next case
 
         rule_expr_1 = "K(S{u})::cyt => K(S{p})::cyt + D(B{_})::cell @ 3*[K(S{u})::cyt]/2*v_1"
-        rule1 = self.parser.parse(rule_expr_1).data[1]
+        rule1 = objects.rule_parser.parse(rule_expr_1).data[1]
 
         rule_expr_2 = "K()::cyt => K()::cyt + D()::cell @ 3*[K()::cyt]/2*v_1"
-        rule2 = self.parser.parse(rule_expr_2).data[1]
+        rule2 = objects.rule_parser.parse(rule_expr_2).data[1]
 
         self.assertEqual(rule1.reduce_context(), rule2)
 
         # next case - covering replication
 
         rule_expr_1 = "K(S{u})::cyt => 2 K(S{u})::cyt @ 3*[K(S{u})::cyt]/2*v_1"
-        rule1 = self.parser.parse(rule_expr_1).data[1]
+        rule1 = objects.rule_parser.parse(rule_expr_1).data[1]
 
         rule_expr_2 = "K()::cyt => 2 K()::cyt @ 3*[K()::cyt]/2*v_1"
-        rule2 = self.parser.parse(rule_expr_2).data[1]
+        rule2 = objects.rule_parser.parse(rule_expr_2).data[1]
 
         self.assertEqual(rule1.reduce_context(), rule2)
 
         # next case - covering replication
 
         rule_expr_1 = "K(S{u})::cyt => 3 K(S{u})::cyt @ 3*[K(S{u})::cyt]/2*v_1"
-        rule1 = self.parser.parse(rule_expr_1).data[1]
+        rule1 = objects.rule_parser.parse(rule_expr_1).data[1]
 
         rule_expr_2 = "K()::cyt => 3 K()::cyt @ 3*[K()::cyt]/2*v_1"
-        rule2 = self.parser.parse(rule_expr_2).data[1]
+        rule2 = objects.rule_parser.parse(rule_expr_2).data[1]
 
         self.assertEqual(rule1.reduce_context(), rule2)
 
     def test_exists_compatible_agent(self):
-        complex_parser = Parser("rate_complex")
         agent = "K(S{a}).A{a}::cyt"
-        complex = complex_parser.parse(agent).data.children[0]
+        complex = objects.rate_complex_parser.parse(agent).data.children[0]
 
         rule_expr = "K().A{i}::cyt => K().A{a}::cyt"
-        rule = self.parser.parse(rule_expr).data[1]
+        rule = objects.rule_parser.parse(rule_expr).data[1]
 
         self.assertTrue(rule.exists_compatible_agent(complex))
diff --git a/Testing/test_side.py b/Testing/test_side.py
index f0c5d77..795ef19 100644
--- a/Testing/test_side.py
+++ b/Testing/test_side.py
@@ -2,83 +2,49 @@
 import collections
 import numpy as np
 
-from eBCSgen.Core.Structure import StructureAgent
-from eBCSgen.Core.Atomic import AtomicAgent
-from eBCSgen.Core.Complex import Complex
-from eBCSgen.Core.Side import Side
 from eBCSgen.TS.State import State, Memory, Vector
-from eBCSgen.Parsing.ParseBCSL import Parser
 
+import Testing.objects_testing as objects
 
-class TestSide(unittest.TestCase):
-    def setUp(self):
-        self.a1 = AtomicAgent("S", "u")
-        self.a2 = AtomicAgent("S", "p")
-        self.a3 = AtomicAgent("T", "_")
-        self.a4 = AtomicAgent("S", "_")
-        self.a5 = AtomicAgent("T", "a")
-
-        self.s1 = StructureAgent("K", {self.a1})
-        self.s2 = StructureAgent("B", set())
-        self.s3 = StructureAgent("K", {self.a2})
-        self.s4 = StructureAgent("B", set())
-        self.s5 = StructureAgent("D", {self.a3})
-        self.s6 = StructureAgent("B", {self.a5, self.a2})
-
-        self.c1 = Complex(collections.Counter({self.s2: 2}), "cyt")
-        self.c2 = Complex(collections.Counter({self.s3: 1}), "cyt")
-        self.c3 = Complex(collections.Counter({self.s2: 1}), "cyt")
-        self.c4 = Complex(collections.Counter({self.s5: 1}), "cell")
-        self.c5 = Complex(collections.Counter({self.s4: 1}), "cell")
-        self.c6 = Complex(collections.Counter({self.s6: 1}), "cell")
 
-        #  rules
-
-        self.side1 = Side([self.c1, self.c2, self.c4])
-        self.side2 = Side([self.c2, self.c3, self.c4])
-        self.side3 = Side([self.c3, self.c4, self.c2])
-        self.side4 = Side([self.c6, self.c1])
-        self.side5 = Side([self.c5, self.c1])
-        self.side6 = Side([self.c5, self.c1, self.c2])
-        self.side7 = Side([self.c6, self.c1, self.c3, self.c4])
+class TestSide(unittest.TestCase):
 
     def test_eq(self):
-        self.assertEqual(self.side2, self.side3)
-        self.assertNotEqual(self.side1, self.side2)
+        self.assertEqual(objects.side7, objects.side8)
+        self.assertNotEqual(objects.side6, objects.side7)
 
     def test_to_counter(self):
-        self.assertEqual(self.side1.to_counter(), collections.Counter({self.c1: 1, self.c2: 1, self.c4: 1}))
+        self.assertEqual(objects.side6.to_counter(), collections.Counter({objects.counter_c1: 1, objects.counter_c2: 1, objects.counter_c4: 1}))
 
     def test_to_vector(self):
-        ordering = (self.c1, self.c2, self.c3, self.c4)
-        self.assertEqual(self.side2.to_vector(ordering), State(Vector(np.array((0, 1, 1, 1))), Memory(0)))
+        ordering = (objects.counter_c1, objects.counter_c2, objects.counter_c3, objects.counter_c4)
+        self.assertEqual(objects.side7.to_vector(ordering), State(Vector(np.array((0, 1, 1, 1))), Memory(0)))
 
     def test_compatible(self):
-        self.assertTrue(self.side5.compatible(self.side4))
-        self.assertTrue(self.side5.compatible(self.side7))
-        self.assertFalse(self.side6.compatible(self.side4))
+        self.assertTrue(objects.side10.compatible(objects.side9))
+        self.assertTrue(objects.side10.compatible(objects.side12))
+        self.assertFalse(objects.side11.compatible(objects.side9))
 
     def test_exists_compatible_agent(self):
-        self.assertTrue(self.side1.exists_compatible_agent(self.c1))
-        self.assertFalse(self.side1.exists_compatible_agent(self.c3))
+        self.assertTrue(objects.side6.exists_compatible_agent(objects.counter_c1))
+        self.assertFalse(objects.side6.exists_compatible_agent(objects.counter_c3))
 
     def test_create_all_compatible(self):
-        complex_parser = Parser("rate_complex")
-        side_parser = Parser("side")
-
         atomic_signature = {"A": {"+", "-"}, "B": {"HA", "HE"}, "S": {"a", "i"}, "T": {"u", "p"}}
         structure_signature = {"KaiB": {"A", "B"}, "KaiC": {"S", "T"}}
 
         results = set()
         with open("Testing/complexes_1.txt") as file:
             for complex in file.readlines():
-                results.add(complex_parser.parse(complex).data.children[0])
+                complex = complex.strip()
+                results.add(objects.rate_complex_parser.parse(complex).data.children[0])
 
         with open("Testing/complexes_2.txt") as file:
             for complex in file.readlines():
-                results.add(complex_parser.parse(complex).data.children[0])
+                complex = complex.strip()
+                results.add(objects.rate_complex_parser.parse(complex).data.children[0])
 
         side = "KaiC().KaiC().KaiC().KaiC().KaiC().KaiC()::cyt + 2 KaiB().KaiB().KaiB()::cyt"
-        side = side_parser.parse(side).data.to_side()
+        side = objects.side_parser.parse(side).data.to_side()
         output_comples = side.create_all_compatible(atomic_signature, structure_signature)
         self.assertEqual(output_comples, results)
diff --git a/Testing/test_structure.py b/Testing/test_structure.py
index b4c5b89..61e157f 100644
--- a/Testing/test_structure.py
+++ b/Testing/test_structure.py
@@ -1,61 +1,44 @@
 import unittest
-from eBCSgen.Core.Structure import StructureAgent
-from eBCSgen.Core.Atomic import AtomicAgent
 
+import Testing.objects_testing as objects
 
-class TestStructure(unittest.TestCase):
-    def setUp(self):
-        self.a1 = AtomicAgent("T", "s")
-        self.a2 = AtomicAgent("S", "a")
-        self.a3 = AtomicAgent("T", "s")
-        self.a4 = AtomicAgent("T", "_")
-        self.a5 = AtomicAgent("T", "p")
-        self.a6 = AtomicAgent("S", "i")
-
-        self.s1 = StructureAgent("strA", {self.a1, self.a2})
-        self.s2 = StructureAgent("strA", {self.a2, self.a4})
-        self.s3 = StructureAgent("strA", {self.a4})
-        self.s4 = StructureAgent("strD", set())
-        self.s5 = StructureAgent("strA", {self.a2, self.a1})
-
-        # context
-
-        self.s6 = StructureAgent("strA", {self.a3})
-        self.s7 = StructureAgent("strA", {self.a5})
-        self.s8 = StructureAgent("strA", {self.a3, self.a2})
-        self.s9 = StructureAgent("strA", {self.a5, self.a2})
-        self.s10 = StructureAgent("strA", {self.a3, self.a6})
-        self.s11 = StructureAgent("strA", {self.a5, self.a6})
-        self.s12 = StructureAgent("strA", set())
 
+class TestStructure(unittest.TestCase):
     def test_eq(self):
-        self.assertEqual(self.s1, self.s5)
-        self.assertNotEqual(self.s1, self.s2)
-        self.assertNotEqual(self.s1, self.a2)  # comparing different classes
+        self.assertEqual(objects.s15, objects.s19)
+        self.assertNotEqual(objects.s15, objects.s16)
+        self.assertNotEqual(objects.s15, objects.a2)  # comparing different classes
 
     def test_print(self):
-        self.assertEqual(str(self.s1), "strA(S{a},T{s})")
-        self.assertEqual(str(self.s4), "strD()")
+        self.assertEqual(str(objects.s15), "strA(S{a},T{s})")
+        self.assertEqual(str(objects.s18), "strD()")
 
     def test_compatibility(self):
-        self.assertTrue(self.s2.compatible(self.s1))
-        self.assertFalse(self.s1.compatible(self.s2))
-        self.assertTrue(self.s3.compatible(self.s1))
-        self.assertFalse(self.s1.compatible(self.a2))  # comparing different classes
+        self.assertTrue(objects.s16.compatible(objects.s15))
+        self.assertFalse(objects.s15.compatible(objects.s16))
+        self.assertFalse(
+            objects.s15.compatible(objects.a2)
+        )  # comparing different classes
 
     def test_add_context(self):
         atomic_signature = {"T": {"s", "p"}, "S": {"i", "a"}}
         structure_signature = {"strA": {"T", "S"}}
-        self.assertEqual(self.s6.add_context(self.s7, atomic_signature, structure_signature),
-                         {(self.s8, self.s9), (self.s10, self.s11)})
-        self.assertEqual(self.s6.add_context(1, atomic_signature, structure_signature),
-                         {(self.s8, None), (self.s10, None)})
-        self.assertEqual(self.s6.add_context(-1, atomic_signature, structure_signature),
-                         {(None, self.s8), (None, self.s10)})
+        self.assertEqual(
+            objects.s20.add_context(objects.s21, atomic_signature, structure_signature),
+            {(objects.s15, objects.s23), (objects.s24, objects.s25)},
+        )
+        self.assertEqual(
+            objects.s20.add_context(1, atomic_signature, structure_signature),
+            {(objects.s15, None), (objects.s24, None)},
+        )
+        self.assertEqual(
+            objects.s20.add_context(-1, atomic_signature, structure_signature),
+            {(None, objects.s15), (None, objects.s24)},
+        )
 
     def test_reduce_context(self):
-        self.assertEqual(self.s11.reduce_context(), self.s12)
+        self.assertEqual(objects.s25.reduce_context(), objects.s26)
 
     def test_replace(self):
-        self.assertEqual(self.s7.replace(self.s8), self.s9)
-        self.assertEqual(self.s12.replace(self.s8), self.s8)
+        self.assertEqual(objects.s21.replace(objects.s15), objects.s23)
+        self.assertEqual(objects.s26.replace(objects.s15), objects.s15)
diff --git a/Testing/test_vector_model.py b/Testing/test_vector_model.py
index 8ca7ec9..7e69b01 100644
--- a/Testing/test_vector_model.py
+++ b/Testing/test_vector_model.py
@@ -13,18 +13,14 @@
 from eBCSgen.TS.VectorModel import VectorModel
 from eBCSgen.TS.VectorReaction import VectorReaction
 
+import Testing.objects_testing as objects
+from Testing.models.get_model_str import get_model_str
+
 
 class TestVectorModel(unittest.TestCase):
     def setUp(self):
-        self.s1 = StructureAgent("X", set())
-        self.s2 = StructureAgent("Y", set())
-        self.s3 = StructureAgent("Z", set())
-
-        self.c1 = Complex([self.s1], "rep")
-        self.c2 = Complex([self.s2], "rep")
-        self.c3 = Complex([self.s3], "rep")
 
-        ordering = (self.c1, self.c2, self.c3)
+        ordering = (objects.c27, objects.c28, objects.c29)
         params = {"k1": 0.05, "k2": 0.1}
 
         self.rate_parser = Parser("rate")
@@ -70,55 +66,14 @@ def setUp(self):
 
         self.model_parser = Parser("model")
 
-        self.model_abstract = \
-            """#! rules
-            => X()::rep @ k2*[T{_}::rep]
-            T{a}::rep => T{i}::rep @ k1*[T{_}::rep]
-
-            #! inits
-            10 T{a}::rep
-
-            #! definitions
-            k1 = 0.05
-            k2 = 0.12
-            """
+        self.model_abstract = get_model_str("model_abstract")
 
         # test transition system generating
 
-        a1 = AtomicAgent("B", "a")
-        a2 = AtomicAgent("S", "u")
-        a3 = AtomicAgent("S", "p")
-        a4 = AtomicAgent("T", "i")
-
-        s1 = StructureAgent("K", {a3, a4})
-        s2 = StructureAgent("K", {a2, a4})
-
-        cx1 = Complex([a1], "cyt")
-        cx2 = Complex([s1], "cyt")
-        cx3 = Complex([s2], "cyt")
-        cx4 = Complex([s1, a1], "cyt")
-        cx5 = Complex([s2, a1], "cyt")
-
-        ordering = (cx5, cx4, cx3, cx2, cx1)
+        ordering = (objects.cx5, objects.cx4, objects.cx3, objects.cx2, objects.cx1)
         # (K(S{u},T{i}).B{a}::cyt, K(S{p},T{i}).B{a}::cyt, K(S{u},T{i})::cyt, K(S{p},T{i})::cyt, B{a}::cyt)
 
-        self.model_TS = \
-            """#! rules
-            => K(S{u},T{i})::cyt @ omega
-            K(S{u})::cyt => K(S{p})::cyt @ alpha*[K(S{u})::cyt]
-            K(S{p})::cyt + B{a}::cyt => K(S{p}).B{a}::cyt @ beta*[K(S{p})::cyt]*[B{a}::cyt]
-            B{_}::cyt => @ gamma*[B{_}::cyt]
-            K(S{u},T{i}).B{a}::cyt => @ 5
-
-            #! inits
-            1 B{a}::cyt
-
-            #! definitions
-            alpha = 10
-            beta = 5
-            gamma = 2
-            omega = 3
-            """
+        self.model_TS = get_model_str("model_TS")
 
         alpha = 10
         beta = 5
@@ -173,76 +128,15 @@ def setUp(self):
 
         # bigger TS
 
-        self.model_bigger_TS = \
-            """#! rules
-            => 2 K(S{u},T{i})::cyt @ omega
-            K(S{u})::cyt => K(S{p})::cyt @ alpha*[K(S{u})::cyt]
-            K(S{p})::cyt + B{a}::cyt => K(S{p}).B{a}::cyt @ beta*[K(S{p})::cyt]*[B{a}::cyt]
-            B{_}::cyt => @ gamma*[B{_}::cyt]
-            K(S{u},T{i}).B{a}::cyt => @ 5
+        self.model_bigger_TS = get_model_str("model_bigger_TS")
 
-            #! inits
-            6 B{a}::cyt
+        # even bigger TS
 
-            #! definitions
-            alpha = 10
-            beta = 5
-            gamma = 2
-            omega = 3
-            """
+        self.model_even_bigger_TS = get_model_str("model_even_bigger_TS")
 
-        # even bigger TS
+        self.model_parametrised = get_model_str("model_parametrised2")
 
-        self.model_even_bigger_TS = \
-            """#! rules
-            => K(S{u},T{i})::cyt @ omega
-            K(S{u})::cyt => K(S{p})::cyt @ alpha*[K(S{u})::cyt]
-            K(S{p})::cyt + B{a}::cyt => K(S{p}).B{a}::cyt @ beta*[K(S{p})::cyt]*[B{a}::cyt]
-            B{_}::cyt => @ gamma*[B{_}::cyt]
-            K(S{u},T{i}).B{a}::cyt => @ 5
-
-            #! inits
-            10 B{a}::cyt
-
-            #! definitions
-            alpha = 10
-            beta = 5
-            gamma = 2
-            omega = 3
-            """
-
-        self.model_parametrised = \
-            """#! rules
-            => K(S{u},T{i})::cyt @ omega
-            K(S{u})::cyt => K(S{p})::cyt @ alpha*[K(S{u})::cyt]
-            K(S{p})::cyt + B{a}::cyt => K(S{p}).B{a}::cyt @ beta*[K(S{p})::cyt]*[B{a}::cyt]
-            B{_}::cyt => @ gamma*[B{_}::cyt]
-            K(S{u},T{i}).B{a}::cyt => @ 5
-
-            #! inits
-            1 B{a}::cyt
-
-            #! definitions
-            alpha = 10
-            beta = 5
-            //gamma = 2
-            omega = 3
-            """
-
-        self.model_with_sinks = \
-            """#! rules
-            K(S{u})::cyt => K(S{p})::cyt @ alpha*[K(S{u})::cyt]
-            K(S{p})::cyt + B{a}::cyt => K(S{p}).B{a}::cyt @ beta*[K(S{p})::cyt]*[B{a}::cyt]
-            B{a}::cyt => B{i}::cyt @ alpha*[B{_}::cyt]
-
-            #! inits
-            1 B{a}::cyt
-            1 K(S{u})::cyt
-
-            #! definitions
-            alpha = 10
-            beta = 5
-            """
+        self.model_with_sinks = get_model_str("model_with_sinks")
 
     def test_compute_bound(self):
         self.assertEqual(self.vm_1.bound, 2)
diff --git a/conda/environment.yml b/conda/environment.yml
index 424a0fc..0850cc4 100644
--- a/conda/environment.yml
+++ b/conda/environment.yml
@@ -12,7 +12,7 @@ dependencies:
   - scipy
   - sortedcontainers
   - sympy
-  - python-libsbml
+  - python-libsbml==5.20.1
   - lark
   - lark-parser
   - pytest
diff --git a/eBCSgen/Core/Formula.py b/eBCSgen/Core/Formula.py
index 853a415..a68e2a9 100644
--- a/eBCSgen/Core/Formula.py
+++ b/eBCSgen/Core/Formula.py
@@ -1,7 +1,7 @@
 from lark import Transformer, Tree
 
 from eBCSgen.Errors.ComplexOutOfScope import ComplexOutOfScope
-from eBCSgen.Core.Rate import tree_to_string
+from eBCSgen.utils import tree_to_string
 
 
 class Formula:
diff --git a/eBCSgen/Core/Rate.py b/eBCSgen/Core/Rate.py
index 648aeae..ba6acfd 100644
--- a/eBCSgen/Core/Rate.py
+++ b/eBCSgen/Core/Rate.py
@@ -4,6 +4,7 @@
 from sortedcontainers import SortedList
 
 from eBCSgen.TS.State import Vector
+from eBCSgen.utils import tree_to_string
 
 STATIC_MATH = """<kineticLaw><math xmlns="http://www.w3.org/1998/Math/MathML"><apply>{}</apply></math></kineticLaw>"""
 
@@ -252,15 +253,3 @@ def fix_operator(self, node, matches):
         operator = self.operators[matches[1].type]
         return Tree(node, [matches[0], Token(matches[1].type, operator), matches[2]])
 
-
-def tree_to_string(tree):
-    """
-    Recursively constructs a list form given lark tree.
-
-    :param tree: given lark tree
-    :return: list of components
-    """
-    if type(tree) == Tree:
-        return sum(list(map(tree_to_string, tree.children)), [])
-    else:
-        return [str(tree)]
diff --git a/eBCSgen/Core/Rule.py b/eBCSgen/Core/Rule.py
index 9faa67f..2bd5932 100644
--- a/eBCSgen/Core/Rule.py
+++ b/eBCSgen/Core/Rule.py
@@ -15,7 +15,16 @@ def column(lst, index):
 
 
 class Rule:
-    def __init__(self, agents: tuple, mid: int, compartments: list, complexes: list, pairs: list, rate: Rate, label=None):
+    def __init__(
+        self,
+        agents: tuple,
+        mid: int,
+        compartments: list,
+        complexes: list,
+        pairs: list,
+        rate: Rate,
+        label=None,
+    ):
         """
         Class to represent BCSL rule
 
@@ -35,9 +44,15 @@ def __init__(self, agents: tuple, mid: int, compartments: list, complexes: list,
         self.label = label
         self.comment = (False, [])
 
-    def __eq__(self, other: 'Rule'):
-        return self.agents == other.agents and self.mid == other.mid and self.compartments == other.compartments and \
-               self.complexes == other.complexes and self.pairs == other.pairs and str(self.rate) == str(other.rate)
+    def __eq__(self, other: "Rule"):
+        return (
+            self.agents == other.agents
+            and self.mid == other.mid
+            and self.compartments == other.compartments
+            and self.complexes == other.complexes
+            and self.pairs == other.pairs
+            and str(self.rate) == str(other.rate)
+        )
 
     def __repr__(self):
         return str(self)
@@ -47,14 +62,23 @@ def __str__(self):
         rate = " @ " + str(self.rate) if self.rate else ""
         pre_comment, post_comment = "", ""
         if self.comment[1]:
-            comment = "// redundant #{" + ", ".join(list(map(str, self.comment[1]))) + "} "
+            comment = (
+                "// redundant #{" + ", ".join(list(map(str, self.comment[1]))) + "} "
+            )
             pre_comment = comment + "// " if self.comment[0] else ""
             post_comment = " " + comment if not self.comment[0] else ""
 
         label = str(self.label) + " ~ " if self.label else ""
 
-        return label + pre_comment + " + ".join(lhs.to_list_of_strings()) + \
-               " => " + " + ".join(rhs.to_list_of_strings()) + rate + post_comment
+        return (
+            label
+            + pre_comment
+            + " + ".join(lhs.to_list_of_strings())
+            + " => "
+            + " + ".join(rhs.to_list_of_strings())
+            + rate
+            + post_comment
+        )
 
     def __lt__(self, other):
         return str(self) < str(other)
@@ -70,10 +94,12 @@ def get_unique_complexes_from_rule(self) -> dict:
         :return: dict of {Complexes:{SBML codes of all isomorphisms in set}}
         """
         unique_complexes_from_rule = dict()
-        for (f, t) in self.complexes:
-            c = Complex(self.agents[f:t + 1], self.compartments[f])
-            double =  (c, c.to_SBML_species_code())
-            unique_complexes_from_rule[c] = unique_complexes_from_rule.get(c, set()) | {double}
+        for f, t in self.complexes:
+            c = Complex(self.agents[f : t + 1], self.compartments[f])
+            double = (c, c.to_SBML_species_code())
+            unique_complexes_from_rule[c] = unique_complexes_from_rule.get(c, set()) | {
+                double
+            }
         return unique_complexes_from_rule
 
     def create_complexes(self):
@@ -83,8 +109,8 @@ def create_complexes(self):
         :return: two multisets of Complexes represented as object Side
         """
         lhs, rhs = [], []
-        for (f, t) in self.complexes:
-            c = Complex(self.agents[f:t + 1], self.compartments[f])
+        for f, t in self.complexes:
+            c = Complex(self.agents[f : t + 1], self.compartments[f])
             lhs.append(c) if t < self.mid else rhs.append(c)
         return Side(lhs), Side(rhs)
 
@@ -108,7 +134,9 @@ def rate_to_vector(self, ordering, definitions: dict):
         if self.rate:
             self.rate.vectorize(ordering, definitions)
 
-    def create_reactions(self, atomic_signature: dict, structure_signature: dict) -> set:
+    def create_reactions(
+        self, atomic_signature: dict, structure_signature: dict
+    ) -> set:
         """
         Create all possible reactions.
         Decide if rule is of replication type and call corresponding lower level method.
@@ -118,13 +146,21 @@ def create_reactions(self, atomic_signature: dict, structure_signature: dict) ->
         :return: set of created reactions
         """
         unique_lhs_indices = set(column(self.pairs, 0))
-        if len(self.pairs) > 1 and len(unique_lhs_indices) == 1 and None not in unique_lhs_indices:
+        if (
+            len(self.pairs) > 1
+            and len(unique_lhs_indices) == 1
+            and None not in unique_lhs_indices
+        ):
             # should be the replication rule
-            return self._create_replication_reactions(atomic_signature, structure_signature)
+            return self._create_replication_reactions(
+                atomic_signature, structure_signature
+            )
         else:
             return self._create_normal_reactions(atomic_signature, structure_signature)
 
-    def _create_replication_reactions(self, atomic_signature: dict, structure_signature: dict) -> set:
+    def _create_replication_reactions(
+        self, atomic_signature: dict, structure_signature: dict
+    ) -> set:
         """
         Create reaction from rule of special form for replication (A -> 2 A)
 
@@ -144,13 +180,22 @@ def _create_replication_reactions(self, atomic_signature: dict, structure_signat
             # replicate RHS agent n times
             for _ in range(len(self.pairs)):
                 new_agents.append(deepcopy(new_agents[-1]))
-            new_rule = Rule(tuple(new_agents), self.mid, self.compartments,
-                            self.complexes, self.pairs, self.rate, self.label)
+            new_rule = Rule(
+                tuple(new_agents),
+                self.mid,
+                self.compartments,
+                self.complexes,
+                self.pairs,
+                self.rate,
+                self.label,
+            )
             reactions.add(new_rule.to_reaction())
 
         return reactions
 
-    def _create_normal_reactions(self, atomic_signature: dict, structure_signature: dict) -> set:
+    def _create_normal_reactions(
+        self, atomic_signature: dict, structure_signature: dict
+    ) -> set:
         """
         Adds context to all agents and generated all possible combinations.
         Then, new rules with these enhances agents are generated and converted to Reactions.
@@ -160,7 +205,7 @@ def _create_normal_reactions(self, atomic_signature: dict, structure_signature:
         :return: set of created reactions
         """
         results = []
-        for (l, r) in self.pairs:
+        for l, r in self.pairs:
             if l is None:
                 right = -1
                 left = self.agents[r]
@@ -170,17 +215,27 @@ def _create_normal_reactions(self, atomic_signature: dict, structure_signature:
             else:
                 left = self.agents[l]
                 right = self.agents[r]
-            results.append(left.add_context(right, atomic_signature, structure_signature))
+            results.append(
+                left.add_context(right, atomic_signature, structure_signature)
+            )
 
         reactions = set()
         for result in itertools.product(*results):
             new_agents = tuple(filter(None, column(result, 0) + column(result, 1)))
-            new_rule = Rule(new_agents, self.mid, self.compartments, self.complexes, self.pairs, self.rate, self.label)
+            new_rule = Rule(
+                new_agents,
+                self.mid,
+                self.compartments,
+                self.complexes,
+                self.pairs,
+                self.rate,
+                self.label,
+            )
             reactions.add(new_rule.to_reaction())
 
         return reactions
 
-    def compatible(self, other: 'Rule') -> bool:
+    def compatible(self, other: "Rule") -> bool:
         """
         Checks whether Rule is compatible (position-wise) with the other Rule.
         Is done by formaly translating to Reactions (just a better object handling).
@@ -201,7 +256,14 @@ def reduce_context(self):
         """
         new_agents = tuple([agent.reduce_context() for agent in self.agents])
         new_rate = self.rate.reduce_context() if self.rate else None
-        return Rule(new_agents, self.mid, self.compartments, self.complexes, self.pairs, new_rate)
+        return Rule(
+            new_agents,
+            self.mid,
+            self.compartments,
+            self.complexes,
+            self.pairs,
+            new_rate,
+        )
 
     def is_meaningful(self) -> bool:
         """
@@ -231,7 +293,9 @@ def create_all_compatible(self, atomic_signature: dict, structure_signature: dic
         :param structure_signature: given structure signature
         :return: set of all created Complexes
         """
-        return self.to_reaction().create_all_compatible(atomic_signature, structure_signature)
+        return self.to_reaction().create_all_compatible(
+            atomic_signature, structure_signature
+        )
 
     def evaluate_rate(self, state, params):
         """
@@ -242,7 +306,7 @@ def evaluate_rate(self, state, params):
         :return: a real number of the rate
         """
         values = dict()
-        for (state_complex, count) in state.content.value.items():
+        for state_complex, count in state.content.value.items():
             for agent in self.rate_agents:
                 if agent.compatible(state_complex):
                     values[agent] = values.get(agent, 0) + count
@@ -277,16 +341,24 @@ def replace(self, aligned_match):
         # replace respective agents
 
         unique_lhs_indices = set(column(self.pairs, 0))
-        if len(self.pairs) > 1 and len(unique_lhs_indices) == 1 and \
-                None not in unique_lhs_indices and len(aligned_match) == 1:
+        if (
+            len(self.pairs) > 1
+            and len(unique_lhs_indices) == 1
+            and None not in unique_lhs_indices
+            and len(aligned_match) == 1
+        ):
             resulting_rhs = self._replace_replicated_rhs(aligned_match[0])
         else:
             resulting_rhs = self._replace_normal_rhs(aligned_match)
 
         # construct resulting complexes
         output_complexes = []
-        for (f, t) in list(filter(lambda item: item[0] >= self.mid, self.complexes)):
-            output_complexes.append(Complex(resulting_rhs[f - self.mid:t - self.mid + 1], self.compartments[f]))
+        for f, t in list(filter(lambda item: item[0] >= self.mid, self.complexes)):
+            output_complexes.append(
+                Complex(
+                    resulting_rhs[f - self.mid : t - self.mid + 1], self.compartments[f]
+                )
+            )
 
         return Multiset(collections.Counter(output_complexes))
 
@@ -298,7 +370,7 @@ def _replace_normal_rhs(self, aligned_match):
         :return: RHS with replaced agents
         """
         resulting_rhs = []
-        for i, rhs_agent in enumerate(self.agents[self.mid:]):
+        for i, rhs_agent in enumerate(self.agents[self.mid :]):
             if len(aligned_match) <= i:
                 resulting_rhs.append(rhs_agent)
             else:
@@ -329,11 +401,11 @@ def reconstruct_complexes_from_match(self, match):
         :return: multiset of constructed agents
         """
         output_complexes = []
-        for (f, t) in list(filter(lambda item: item[1] < self.mid, self.complexes)):
-            output_complexes.append(Complex(match[f:t + 1], self.compartments[f]))
+        for f, t in list(filter(lambda item: item[1] < self.mid, self.complexes)):
+            output_complexes.append(Complex(match[f : t + 1], self.compartments[f]))
         return Multiset(collections.Counter(output_complexes))
 
-    def create_reversible(self):
+    def create_reversible(self, rate: Rate = None):
         """
         Create a reversible version of the rule with _bw label.
 
@@ -343,19 +415,22 @@ def create_reversible(self):
 
         :return: reversed Rule
         """
-        agents = self.agents[self.mid:] + self.agents[:self.mid]
+        agents = self.agents[self.mid :] + self.agents[: self.mid]
         mid = len(self.agents) - self.mid
-        compartments = self.compartments[self.mid:] + self.compartments[:self.mid]
-        complexes = sorted([((f - self.mid) % len(self.agents),
-                             (t - self.mid) % len(self.agents)) for (f, t) in self.complexes])
+        compartments = self.compartments[self.mid :] + self.compartments[: self.mid]
+        complexes = sorted(
+            [
+                ((f - self.mid) % len(self.agents), (t - self.mid) % len(self.agents))
+                for (f, t) in self.complexes
+            ]
+        )
         pairs = []
-        for (l, r) in self.pairs:
+        for l, r in self.pairs:
             if l is None or r is None:
                 pairs.append((r, l))
             else:
                 pairs.append((l, r))
 
-        rate = self.rate
         label = None
         if self.label:
             label = self.label + "_bw"
diff --git a/eBCSgen/Errors/RegulationParsingError.py b/eBCSgen/Errors/RegulationParsingError.py
new file mode 100644
index 0000000..554976b
--- /dev/null
+++ b/eBCSgen/Errors/RegulationParsingError.py
@@ -0,0 +1,6 @@
+class RegulationParsingError(Exception):
+    def __init__(self, error):
+        self.error = error
+
+    def __str__(self):
+        return "Error while parsing the regulation:\n\n{}".format(self.error)
diff --git a/eBCSgen/Parsing/ParseBCSL.py b/eBCSgen/Parsing/ParseBCSL.py
index 69915b4..8368f44 100644
--- a/eBCSgen/Parsing/ParseBCSL.py
+++ b/eBCSgen/Parsing/ParseBCSL.py
@@ -4,7 +4,7 @@
 from numpy import inf
 from copy import deepcopy
 from lark import Lark, Transformer, Tree
-from lark import UnexpectedCharacters, UnexpectedToken
+from lark import UnexpectedCharacters, UnexpectedToken, UnexpectedEOF
 from lark.load_grammar import _TERMINAL_NAMES
 import regex
 from sortedcontainers import SortedList
@@ -24,6 +24,10 @@
 from eBCSgen.TS.Edge import edge_from_dict
 from eBCSgen.Core.Side import Side
 from eBCSgen.Core.Model import Model
+from eBCSgen.Errors.ComplexParsingError import ComplexParsingError
+from eBCSgen.Errors.UnspecifiedParsingError import UnspecifiedParsingError
+from eBCSgen.Errors.RegulationParsingError import RegulationParsingError
+from eBCSgen.utils import tree_to_string
 
 
 def load_TS_from_json(json_file: str) -> TransitionSystem:
@@ -37,19 +41,27 @@ def load_TS_from_json(json_file: str) -> TransitionSystem:
     with open(json_file) as json_file:
         data = json.load(json_file)
 
-        ordering = SortedList(map(lambda agent: complex_parser.parse(agent).data.children[0], data['ordering']))
-        ts = TransitionSystem(ordering, data['bound'])
+        ordering = SortedList(
+            map(
+                lambda agent: complex_parser.parse(agent).data.children[0],
+                data["ordering"],
+            )
+        )
+        ts = TransitionSystem(ordering, data["bound"])
         ts.states_encoding = dict()
-        for node_id in data['nodes']:
-            vector = np.array(eval(data['nodes'][node_id]))
+        for node_id in data["nodes"]:
+            vector = np.array(eval(data["nodes"][node_id]))
             is_hell = True if vector[0] == inf else False
             ts.states_encoding[int(node_id)] = State(Vector(vector), Memory(0), is_hell)
-        ts.edges = {edge_from_dict(edge) for edge in data['edges']}
-        ts.init = data['initial']
-        if 'parameters' in data:
-            ts.params = data['parameters']
-
-        ts.unprocessed = {State(Vector(np.array(eval(state))), Memory(0)) for state in data.get('unprocessed', list())}
+        ts.edges = {edge_from_dict(edge) for edge in data["edges"]}
+        ts.init = data["initial"]
+        if "parameters" in data:
+            ts.params = data["parameters"]
+
+        ts.unprocessed = {
+            State(Vector(np.array(eval(state))), Memory(0))
+            for state in data.get("unprocessed", list())
+        }
         ts.states = set(ts.states_encoding.values()) - ts.unprocessed
         return ts
 
@@ -76,34 +88,43 @@ def __init__(self):
         self.counter = 0
 
     def __str__(self):
-        return " | ".join([str(self.seq), str(self.comp), str(self.complexes), str(self.counter)])
+        return " | ".join(
+            [str(self.seq), str(self.comp), str(self.complexes), str(self.counter)]
+        )
 
     def __repr__(self):
         return str(self)
 
     def to_side(self):
-        return Side([Complex(self.seq[c[0]:c[1] + 1], self.comp[c[0]]) for c in self.complexes])
+        return Side(
+            [
+                Complex(self.seq[c[0] : c[1] + 1], self.comp[c[0]])
+                for c in self.complexes
+            ]
+        )
 
 
 GRAMMAR = r"""
-    model: rules (inits)? (definitions)? (complexes)? (regulation)?
+    model: (sections)* rules (sections | rules)*
+    sections: inits | definitions | complexes | regulation | observables
 
-    rules: RULES_START (rule|COMMENT)+
-    inits: INITS_START (init|COMMENT)+
-    definitions: DEFNS_START (definition|COMMENT)+
-    complexes: COMPLEXES_START (cmplx_dfn|COMMENT)+
-    regulation: REGULATION_START regulation_def
+    rules: RULES_START _NL+ (rule _NL+)* rule _NL*
+    inits: INITS_START _NL+ (init _NL+)* init _NL*
+    definitions: DEFNS_START _NL+ (definition _NL+)* definition _NL*
+    complexes: COMPLEXES_START _NL+ (cmplx_dfn _NL+)* cmplx_dfn _NL*
+    regulation: REGULATION_START _NL+ regulation_def _NL*
+    observables: OBSERVABLES_START _NL+ (observable _NL+)* observable _NL*
 
-    init: const? rate_complex (COMMENT)?
-    definition: def_param "=" number (COMMENT)?
-    rule: (label)? side ARROW side ("@" rate)? (";" variable)? (COMMENT)?
-    cmplx_dfn: cmplx_name "=" value (COMMENT)?
+    init: const? rate_complex 
+    definition: def_param "=" number
+    rule: ((label)? side arrow side ("@" rate)? (";" variable)?) | ((label)? side BI_ARROW side ("@" rate "|" rate )? (";" variable)?)
+    cmplx_dfn: cmplx_name "=" value
 
     side: (const? complex "+")* (const? complex)?
     complex: (abstract_sequence|value|cmplx_name) DOUBLE_COLON compartment
 
-    !rate : fun "/" fun | fun
-    !fun: const | param | rate_agent | fun "+" fun | fun "-" fun | fun "*" fun | fun POW const | "(" fun ")"
+    !rate : fun
+    !fun: const | param | rate_agent | fun "+" fun | fun "-" fun | fun "*" fun | fun POW const | "(" fun ")" | fun "/" fun
 
     !rate_agent: "[" rate_complex "]"
 
@@ -111,12 +132,17 @@ def to_side(self):
 
     COM: "//"
     POW: "**"
-    ARROW: "=>" | "<=>" 
+    arrow: SINGLE_ARROW | REPLICATION_ARROW
+    SINGLE_ARROW: "=>"
+    BI_ARROW: "<=>"
+    REPLICATION_ARROW: "=*>"
     RULES_START: "#! rules"
     INITS_START: "#! inits"
     DEFNS_START: "#! definitions"
     COMPLEXES_START: "#! complexes"
     REGULATION_START: "#! regulation"
+    OBSERVABLES_START: "#! observables"
+    _NL: /(\r?\n[\t ]*)+/
 
     !label: CNAME "~"
 
@@ -130,16 +156,16 @@ def to_side(self):
     %import common.NUMBER
     %import common.INT
     %import common.DECIMAL
-    %import common.WS
-    %ignore WS
+    %import common.WS_INLINE
+    %ignore WS_INLINE
     %ignore COMMENT
 """
 
 EXTENDED_GRAMMAR = """
     abstract_sequence: atomic_complex | atomic_structure_complex | structure_complex
-    atomic_complex: atomic ":" (cmplx_name|VAR)
-    atomic_structure_complex: atomic ":" structure ":" (cmplx_name|VAR)
-    structure_complex: structure ":" (cmplx_name|VAR)
+    atomic_complex: atomic "::" (VAR | value)
+    atomic_structure_complex: atomic "::" structure "::" (VAR | value)
+    structure_complex: structure "::" (VAR | value)
 
     variable: VAR "=" "{" cmplx_name ("," cmplx_name)+ "}"
     VAR: "?"
@@ -169,32 +195,73 @@ def to_side(self):
 REGULATIONS_GRAMMAR = """
     regulation_def: "type" ( regular | programmed | ordered | concurrent_free | conditional ) 
 
-    !regular: "regular" (DIGIT|LETTER| "+" | "*" | "(" | ")" | "[" | "]" | "_" | "|" | "&")+
+    !regular: "regular" _NL+ expression _NL*
 
-    programmed: "programmed" successors+
+    programmed: "programmed" _NL+ (successors _NL+)* successors _NL*
     successors: CNAME ":" "{" CNAME ("," CNAME)* "}"
 
-    ordered: "ordered" order ("," order)*
+    ordered: "ordered" _NL+ order ("," order)* _NL*
     order: ("(" CNAME "," CNAME ")")
 
-    concurrent_free: "concurrent-free" order ("," order)*
+    concurrent_free: "concurrent-free" _NL+ order ("," order)* _NL*
 
-    conditional: "conditional" context+
+    conditional: "conditional" _NL+ context (_NL+ context)* _NL*
     context: CNAME ":" "{" rate_complex ("," rate_complex)* "}"
 """
 
+REGEX_GRAMMAR = r"""
+    !?expression: term ("|" term)*
+
+    ?term: factor+
+
+    ?factor: primary quantifier?
+
+    !quantifier: "??"
+            | "*?"
+            | "+?"
+            | "*+"
+            | "++"
+            | "?+"
+            | "*"
+            | "+"
+            | "?"
+            | "{NUMBER,NUMBER}"
+            | "{NUMBER}"
+
+    !primary: "(" expression ")"
+        | "[" REGEX_CHAR "-" REGEX_CHAR "]"
+        | "[" REGEX_CHAR* "]"
+        | SPECIAL_CHAR
+        | ESCAPED_CHAR
+        | "." 
+        | REGEX_CHAR
+
+    SPECIAL_CHAR: "^" | "$" | "&"
+
+    ESCAPED_CHAR: "\\" ("w"|"W"|"d"|"D"|"s"|"S"|"b"|"B"|"A"|"Z"|"G"|"."|"^"|"["|"]"|"("|")"|"{"|"}"|"?"|"*"|"+"|"|"|"\\")
+
+    REGEX_CHAR: /[^\\^$().*+?{}\[\]|]/
+"""
+
+OBSERVABLES_GRAMMAR = """
+    observable: CNAME ":" observable_pattern
+    !observable_pattern: const | complex | observable_pattern "+" observable_pattern | observable_pattern "-" observable_pattern | observable_pattern "*" observable_pattern | observable_pattern "/" observable_pattern | observable_pattern POW const | "(" observable_pattern ")"
+"""
+
 
 class TransformRegulations(Transformer):
     def regulation(self, matches):
-        return {'regulation': matches[1]}
+        return {"regulation": matches[1]}
 
     def regulation_def(self, matches):
         return matches[0]
 
     def regular(self, matches):
-        re = "".join(matches[1:])
-        # might raise exception
-        regex.compile(re)
+        re = "".join(tree_to_string(matches[1]))
+        try:
+            regex.compile(re)
+        except regex.error as e:
+            raise RegulationParsingError(f"Invalid regular expression: {re}. Error: {e}")
         return Regular(re)
 
     def programmed(self, matches):
@@ -252,27 +319,28 @@ def cmplx_dfn(self, matches):
     def rules(self, matches):
         new_rules = [matches[0]]
         for rule in matches[1:]:
-            if rule.children[-1].data == 'variable':
+            if rule.children[-1].data == "variable":
                 variables = rule.children[-1].children[1:]
                 for variable in variables:
                     replacer = ReplaceVariables(variable)
-                    new_rule = Tree('rule', deepcopy(rule.children[:-1]))
+                    new_rule = Tree("rule", deepcopy(rule.children[:-1]))
                     new_rules.append(replacer.transform(new_rule))
             else:
                 new_rules.append(rule)
-        return Tree('rules', new_rules)
+        return Tree("rules", new_rules)
 
 
 def remove_nested_complex_aliases(complex_defns):
     """
     Removes nested complex aliases from their definitions.
     """
+
     def replace_aliases(complex_defns):
         new_definitions = dict()
         for defn in complex_defns:
             result = []
             for child in complex_defns[defn]:
-                if child.data == 'cmplx_name':
+                if child.data == "cmplx_name":
                     result += complex_defns[str(child.children[0])]
                 else:
                     result.append(child)
@@ -286,7 +354,7 @@ def replace_aliases(complex_defns):
         complex_defns = new_defns
         new_defns = replace_aliases(complex_defns)
 
-    complex_defns = {key: Tree('value', complex_defns[key]) for key in complex_defns}
+    complex_defns = {key: Tree("value", complex_defns[key]) for key in complex_defns}
     return complex_defns
 
 
@@ -303,6 +371,10 @@ def __init__(self, complex_defns):
         self.complex_defns = complex_defns
 
     def cmplx_name(self, matches):
+        if str(matches[0]) not in self.complex_defns:
+            raise ComplexParsingError(
+                f"Complex alias {matches[0]} not found in defined complexes: {list(self.complex_defns.keys())}", matches
+            )
         return deepcopy(self.complex_defns[str(matches[0])])
 
     def abstract_sequence(self, matches):
@@ -332,33 +404,120 @@ def structure_complex(self, matches):
 
     def insert_atomic_to_struct(self, atomic, struct):
         """
-        Adds or replaces atomic subtree in struct tree.
+        Adds an atomic subtree to a struct tree. If a non-empty atomic with the same name
+        already exists in the struct, it raises an error to prevent illegal nesting.
+
+        Args:
+            atomic: The atomic subtree to be added.
+            struct: The struct tree where the atomic will be added.
+
+        Returns:
+            The updated struct tree with the atomic added.
+
+        Raises:
+            ComplexParsingError: If a non-empty atomic with the same name is already present in the struct.
         """
         if len(struct.children) == 2:
+            for i in range(len(struct.children[1].children)):
+                if self.get_name(atomic) == self.get_name(
+                    struct.children[1].children[i]
+                ):
+                    if self.is_empty(struct.children[1].children[i]):
+                        struct.children[1].children[i] = atomic
+                        return struct
+                    raise ComplexParsingError(
+                        f"Illegal nesting sequence: {atomic}:{struct}", struct
+                    )
             struct.children[1].children.append(atomic)
         else:
-            struct.children.append(Tree('composition', [atomic]))
+            struct.children.append(Tree("composition", [atomic]))
         return struct
 
     def insert_struct_to_complex(self, struct, complex):
         """
-        Adds or replaces struct subtree in complex tree.
+        Adds a struct subtree to a complex tree, or merges it with an existing struct subtree. This method first
+        searches for a struct in the complex with the same name as the input struct. If found, it then checks for
+        atomic incompatibility within the structs.
+
+        The method ensures that the struct being added does not contain atomics with names that match any atomics
+        in the corresponding struct in the complex. This step is crucial to maintain the integrity of the complex
+        by avoiding conflicting or duplicate atomic structures.
+
+        Args:
+            struct: The struct subtree to be added or merged.
+            complex: The complex tree where the struct will be added or merged.
+
+        Returns:
+            The updated complex tree with the struct added or merged.
+
+        Raises:
+            ComplexParsingError: If no matching struct is found in the complex.
         """
-        for i in range(len(complex.children)):
-            if self.get_name(struct) == self.get_name(complex.children[i].children[0]):
-                complex.children[i] = Tree('agent', [struct])
-                break
-        return complex
+        if isinstance(complex.children[0].children[0].children[0].children[0], Tree):
+            search = complex.children[0]
+        else:
+            search = complex
+
+        for i in range(len(search.children)):
+            if self.get_name(struct) == self.get_name(search.children[i].children[0]):
+                struct_found = True
+                # search same name structs - if they contain atomics with matching names, they are considered incompatible
+                for j in range(len(struct.children[1].children)):
+                    for k in range(
+                        len(search.children[i].children[0].children[1].children)
+                    ):
+                        if self.get_name(
+                            struct.children[1].children[j]
+                        ) == self.get_name(
+                            search.children[i].children[0].children[1].children[k]
+                        ):
+                            struct_found = False
+                            break
+                    # if no same name atomic found in the struct, we found the suitable complex's struct
+                    if not struct_found:
+                        break
+
+                if struct_found:
+                    # if the complex's struct is empty, replace it with the struct
+                    if self.is_empty(search.children[i]):
+                        search.children[i] = Tree("agent", [struct])
+                    else:
+                        # if the complex's struct is not empty merge the struct's children into the complex's struct
+                        search.children[i].children[0].children[1].children += struct.children[1].children
+                    return complex
+
+        raise ComplexParsingError(
+            f"Illegal struct nesting or duplication: {struct}:{complex}", complex
+        )
 
     def insert_atomic_to_complex(self, atomic, complex):
         """
-        Adds or replaces atomic subtree in complex tree.
+        Adds an atomic subtree to a complex tree. If a non-empty atomic with the same name
+        is already present in the complex, it raises an error to prevent illegal nesting.
+
+        Args:
+            atomic: The atomic subtree to be added.
+            complex: The complex tree where the atomic will be added.
+
+        Returns:
+            The updated complex tree with the atomic added.
+
+        Raises:
+            ComplexParsingError: If an atomic with the same name is already present in the complex.
         """
-        for i in range(len(complex.children)):
-            if self.get_name(atomic) == self.get_name(complex.children[i].children[0]):
-                complex.children[i] = Tree('agent', [atomic])
-                break
-        return complex
+        if isinstance(complex.children[0].children[0].children[0].children[0], Tree):
+            search = complex.children[0]
+        else:
+            search = complex
+
+        for i in range(len(search.children)):
+            if self.get_name(atomic) == self.get_name(search.children[i].children[0]):
+                if self.is_empty(search.children[i].children[0]):
+                    search.children[i] = Tree("agent", [atomic])
+                    return complex
+        raise ComplexParsingError(
+            f"Illegal atomic nesting or duplication: {atomic}:{complex}", complex
+        )
 
     def get_name(self, agent):
         return str(agent.children[0].children[0])
@@ -366,11 +525,21 @@ def get_name(self, agent):
     def complex(self, matches):
         result = []
         for match in matches[0].children:
-            if match.data == 'value':
+            if match.data == "value":
                 result += match.children
             else:
                 result.append(match)
-        return Tree('complex', [Tree('value', result)] + matches[1:])
+        return Tree("complex", [Tree("value", result)] + matches[1:])
+
+    def is_empty(self, agent):
+        """
+        Checks if the agent is empty.
+        """
+        if agent.data == "atomic":
+            return agent.children[1].children[0] == "_"
+        elif agent.data == "agent":
+            return len(agent.children[0].children[1].children) == 0
+        return False
 
 
 class TreeToComplex(Transformer):
@@ -388,11 +557,19 @@ def atomic(self, matches):
 
     def structure(self, matches):
         name = str(matches[0].children[0])
-        if len(matches) > 1:
-            composition = set(matches[1].children)
-            return StructureAgent(name, composition)
-        else:
+        if len(matches) <= 1:
             return StructureAgent(name, set())
+        atomic_names = set()
+        composition = set()
+        for atomic in matches[1].children:
+            if atomic.name in atomic_names:
+                raise ComplexParsingError(
+                    f"Duplicate atomic agent in structure: {atomic.name}", matches
+                )
+            atomic_names.add(atomic.name)
+            composition.add(atomic)
+
+        return StructureAgent(name, composition)
 
     def rate_complex(self, matches):
         sequence = []
@@ -409,6 +586,7 @@ class TreeToObjects(Transformer):
     def __init__(self):
         super(TreeToObjects, self).__init__()
         self.params = set()
+        self.labels = set()
 
     """
     A transformer which is called on a tree in a bottom-up manner and transforms all subtrees/tokens it encounters.
@@ -450,57 +628,84 @@ def side(self, matches):
 
     def rule(self, matches):
         label = None  # TODO create implicit label
-        rate = None
+        rate1 = None
+        rate2 = None
+        if len(matches) == 6:
+            label, lhs, arrow, rhs, rate1, rate2 = matches
         if len(matches) == 5:
-            label, lhs, arrow, rhs, rate = matches
+            if type(matches[0]) == str:
+                label, lhs, arrow, rhs, rate1 = matches
+            else:
+                lhs, arrow, rhs, rate1, rate2 = matches
         elif len(matches) == 4:
             if type(matches[0]) == str:
                 label, lhs, arrow, rhs = matches
             else:
-                lhs, arrow, rhs, rate = matches
+                lhs, arrow, rhs, rate1 = matches
         else:
             lhs, arrow, rhs = matches
+        if label:
+            self.labels.add(label)
         agents = tuple(lhs.seq + rhs.seq)
         mid = lhs.counter
         compartments = lhs.comp + rhs.comp
         complexes = lhs.complexes + list(
-            map(lambda item: (item[0] + lhs.counter, item[1] + lhs.counter), rhs.complexes))
+            map(
+                lambda item: (item[0] + lhs.counter, item[1] + lhs.counter),
+                rhs.complexes,
+            )
+        )
         pairs = [(i, i + lhs.counter) for i in range(min(lhs.counter, rhs.counter))]
-        if lhs.counter > rhs.counter:
+        if type(arrow) is Tree and arrow.children[0].value == "=*>":
+            if lhs.counter >= rhs.counter or lhs.counter != 1 or rhs.counter <= 1:
+                raise UnspecifiedParsingError("Rule does not contain replication")
+            
+            for i in range(lhs.counter, rhs.counter):
+                if lhs.seq[pairs[-1][0]] == rhs.seq[pairs[-1][1] - lhs.counter]:
+                    if rhs.seq[pairs[-1][1] - lhs.counter] == rhs.seq[i]:
+                        pairs += [(pairs[-1][0], i + lhs.counter)]
+                    else:
+                        raise UnspecifiedParsingError("Rule does not contain replication")
+
+        elif lhs.counter > rhs.counter:
             pairs += [(i, None) for i in range(rhs.counter, lhs.counter)]
         elif lhs.counter < rhs.counter:
-            for i in range(lhs.counter, rhs.counter):
-                replication = False
-                if lhs.counter == 1 and rhs.counter > 1:
-                    if lhs.seq[pairs[-1][0]] == rhs.seq[pairs[-1][1] - lhs.counter]:
-                        if rhs.seq[pairs[-1][1] - lhs.counter] == rhs.seq[i]:
-                            pairs += [(pairs[-1][0], i + lhs.counter)]
-                            replication = True
-                if not replication:
-                    pairs += [(None, i + lhs.counter)]
+            pairs += [(None, i + lhs.counter) for i in range(lhs.counter, rhs.counter)]
 
         reversible = False
-        if arrow == '<=>':
+        if arrow == "<=>":
             reversible = True
-        return reversible, Rule(agents, mid, compartments, complexes, pairs, Rate(rate) if rate else None, label)
+        return (
+            reversible,
+            Rule(
+                agents,
+                mid,
+                compartments,
+                complexes,
+                pairs,
+                Rate(rate1) if rate1 else None,
+                label,
+            ),
+            Rate(rate2) if rate2 else None,
+        )
 
     def rules(self, matches):
         rules = []
-        for reversible, rule in matches[1:]:
+        for reversible, rule, new_rate in matches[1:]:
             if reversible:
-                reversible_rule = rule.create_reversible()
+                reversible_rule = rule.create_reversible(new_rate)
                 rules.append(rule)
                 rules.append(reversible_rule)
             else:
                 rules.append(rule)
-        return {'rules': rules}
+        return {"rules": rules}
 
     def definitions(self, matches):
         result = dict()
         for definition in matches[1:]:
             pair = definition.children
             result[pair[0]] = pair[1]
-        return {'definitions': result}
+        return {"definitions": result}
 
     def init(self, matches):
         return matches
@@ -512,53 +717,93 @@ def inits(self, matches):
                 result[init[1].children[0]] = int(init[0])
             else:
                 result[init[0].children[0]] = 1
-        return {'inits': result}
+        return {"inits": result}
+
+    def observable(self, matches):
+        return {str(matches[0]): matches[1].children}
+
+    def observables(self, matches):
+        result = dict()
+        for observable in matches[1:]:
+            result.update(observable)
+        return {"observables": result}
 
     def param(self, matches):
         self.params.add(str(matches[0]))
         return Tree("param", matches)
 
     def model(self, matches):
+        rules = set()
         definitions = dict()
         regulation = None
         inits = collections.Counter()
+        observables = dict()
         for match in matches:
             if type(match) == dict:
                 key, value = list(match.items())[0]
-                if key == 'rules':
-                    rules = set(value)
-                if key == 'inits':
-                    inits = value
-                if key == 'definitions':
-                    definitions = value
-                if key == 'regulation':
+            elif isinstance(match, Tree) and match.data == "sections":
+                if isinstance(match.children[0], Tree):
+                    continue
+                key, value = list(match.children[0].items())[0]
+            else:
+                continue
+
+            if key == "rules":
+                rules.update(value)
+            elif key == "inits":
+                inits.update(value)
+            elif key == "definitions":
+                definitions.update(value)
+            elif key == "observables":
+                observables.update(value)
+            elif key == "regulation":
+                if regulation:
+                    raise UnspecifiedParsingError("Multiple regulations")
+                # check if regulation is in label
+                if value.check_labels(self.labels):
                     regulation = value
+
         params = self.params - set(definitions)
         return Model(rules, inits, definitions, params, regulation)
 
 
 class Parser:
     def __init__(self, start):
-        grammar = "start: " + start + GRAMMAR + COMPLEX_GRAMMAR + EXTENDED_GRAMMAR + REGULATIONS_GRAMMAR
-        self.parser = Lark(grammar, parser='lalr',
-                           propagate_positions=False,
-                           maybe_placeholders=False
-                           )
+        grammar = (
+            "start: "
+            + start
+            + GRAMMAR
+            + COMPLEX_GRAMMAR
+            + EXTENDED_GRAMMAR
+            + REGULATIONS_GRAMMAR
+            + REGEX_GRAMMAR
+            + OBSERVABLES_GRAMMAR
+        )
+        self.parser = Lark(
+            grammar,
+            parser="earley",
+            propagate_positions=False,
+            maybe_placeholders=False,
+        )
 
         self.terminals = dict((v, k) for k, v in _TERMINAL_NAMES.items())
-        self.terminals.update({"COM": "//",
-                               "ARROW": "=>, <=>",
-                               "POW": "**",
-                               "DOUBLE_COLON": "::",
-                               "RULES_START": "#! rules",
-                               "INITS_START": "#! inits",
-                               "DEFNS_START": "#! definitions",
-                               "COMPLEXES_START": "#! complexes",
-                               "REGULATION_START": "#! regulation",
-                               "CNAME": "name",
-                               "NAME": "agent_name",
-                               "VAR": "?"
-                               })
+        self.terminals.update(
+            {
+                "COM": "//",
+                "ARROW": "=>",
+                "BI_ARROW": "<=>",
+                "POW": "**",
+                "DOUBLE_COLON": "::",
+                "RULES_START": "#! rules",
+                "INITS_START": "#! inits",
+                "DEFNS_START": "#! definitions",
+                "COMPLEXES_START": "#! complexes",
+                "REGULATION_START": "#! regulation",
+                "CNAME": "name",
+                "NAME": "agent_name",
+                "VAR": "?",
+            }
+        )
         self.terminals.pop("$END", None)
 
     def replace(self, expected: set) -> set:
@@ -596,7 +841,9 @@ def transform(self, tree: Tree) -> Result:
         try:
             complexer = ExtractComplexNames()
             tree = complexer.transform(tree)
-            complexer.complex_defns = remove_nested_complex_aliases(complexer.complex_defns)
+            complexer.complex_defns = remove_nested_complex_aliases(
+                complexer.complex_defns
+            )
             de_abstracter = TransformAbstractSyntax(complexer.complex_defns)
             tree = de_abstracter.transform(tree)
             tree = TreeToComplex().transform(tree)
@@ -618,11 +865,33 @@ def syntax_check(self, expression: str) -> Result:
         try:
             tree = self.parser.parse(expression)
         except UnexpectedCharacters as u:
-            return Result(False, {"unexpected": expression[u.pos_in_stream],
-                                  "expected": self.replace(u.allowed),
-                                  "line": u.line, "column": u.column})
+            return Result(
+                False,
+                {
+                    "unexpected": expression[u.pos_in_stream],
+                    "expected": self.replace(u.allowed),
+                    "line": u.line,
+                    "column": u.column,
+                },
+            )
         except UnexpectedToken as u:
-            return Result(False, {"unexpected": str(u.token),
-                                  "expected": self.replace(u.expected),
-                                  "line": u.line, "column": u.column})
+            return Result(
+                False,
+                {
+                    "unexpected": str(u.token),
+                    "expected": self.replace(u.expected),
+                    "line": u.line,
+                    "column": u.column,
+                },
+            )
+        except UnexpectedEOF as u:
+            return Result(
+                False,
+                {
+                    "unexpected": str(u.token),
+                    "expected": self.replace(u.expected),
+                    "line": u.line,
+                    "column": u.column,
+                },
+            )
         return Result(True, tree)
diff --git a/eBCSgen/Regulations/ConcurrentFree.py b/eBCSgen/Regulations/ConcurrentFree.py
index ad5022e..282f48f 100644
--- a/eBCSgen/Regulations/ConcurrentFree.py
+++ b/eBCSgen/Regulations/ConcurrentFree.py
@@ -1,3 +1,4 @@
+from eBCSgen.Errors.RegulationParsingError import RegulationParsingError
 from eBCSgen.Regulations.Base import BaseRegulation
 
 
@@ -22,3 +23,10 @@ def filter(self, current_state, candidates):
             if p_rule and non_p_rule:
                 del candidates[non_p_rule.pop()]
         return candidates
+
+    def check_labels(self, model_labels):
+        for tuple in self.regulation:
+            for label in tuple:
+                if label not in model_labels:
+                    raise RegulationParsingError(f"Label {label} in concurrent-free regulation not present in model")
+        return True
\ No newline at end of file
diff --git a/eBCSgen/Regulations/Conditional.py b/eBCSgen/Regulations/Conditional.py
index a11998f..6f6914b 100644
--- a/eBCSgen/Regulations/Conditional.py
+++ b/eBCSgen/Regulations/Conditional.py
@@ -1,3 +1,4 @@
+from eBCSgen.Errors.RegulationParsingError import RegulationParsingError
 from eBCSgen.Regulations.Base import BaseRegulation
 
 
@@ -19,6 +20,12 @@ def __repr__(self):
     def filter(self, current_state, candidates):
         agents = set(current_state.content.value)
         return {rule: values for rule, values in candidates.items() if not self.check_intersection(rule.label, agents)}
+    
+    def check_labels(self, model_labels):
+        for label in self.regulation:
+            if label not in model_labels:
+                raise RegulationParsingError(f"Label {label} in conditional regulation not present in model")
+        return True
 
     def check_intersection(self, label, agents):
         if label not in self.regulation:
diff --git a/eBCSgen/Regulations/Ordered.py b/eBCSgen/Regulations/Ordered.py
index bb37018..5bae8e8 100644
--- a/eBCSgen/Regulations/Ordered.py
+++ b/eBCSgen/Regulations/Ordered.py
@@ -1,3 +1,4 @@
+from eBCSgen.Errors.RegulationParsingError import RegulationParsingError
 from eBCSgen.Regulations.Base import BaseRegulation
 
 
@@ -38,3 +39,10 @@ def filter(self, current_state, candidates):
             return candidates
         last_rule = current_state.memory.history[-1]
         return {rule: values for rule, values in candidates.items() if not (last_rule, rule.label) in self.regulation}
+    
+    def check_labels(self, model_labels):
+        for tuple in self.regulation:
+            for label in tuple:
+                if label not in model_labels:
+                    raise RegulationParsingError(f"Label {label} in programmed regulation not present in model")
+        return True
\ No newline at end of file
diff --git a/eBCSgen/Regulations/Programmed.py b/eBCSgen/Regulations/Programmed.py
index e5580e5..8492053 100644
--- a/eBCSgen/Regulations/Programmed.py
+++ b/eBCSgen/Regulations/Programmed.py
@@ -1,3 +1,4 @@
+from eBCSgen.Errors.RegulationParsingError import RegulationParsingError
 from eBCSgen.Regulations.Base import BaseRegulation
 
 
@@ -24,3 +25,12 @@ def filter(self, current_state, candidates):
         if last_rule in self.regulation:
             return {rule: values for rule, values in candidates.items() if rule.label in self.regulation[last_rule]}
         return candidates
+    
+    def check_labels(self, model_labels):
+        for rule_label, successors_labels in self.regulation.items():
+            if rule_label not in model_labels:
+                raise RegulationParsingError(f"Label {rule_label} in programmed regulation not present in model")
+            if not successors_labels.issubset(model_labels):
+                missing_labels = successors_labels - model_labels
+                raise RegulationParsingError(f"Label(s) {missing_labels} in programmed regulation not present in model")
+        return True
\ No newline at end of file
diff --git a/eBCSgen/Regulations/Regular.py b/eBCSgen/Regulations/Regular.py
index 2b8f859..7b7497b 100644
--- a/eBCSgen/Regulations/Regular.py
+++ b/eBCSgen/Regulations/Regular.py
@@ -1,4 +1,5 @@
 import regex
+from eBCSgen.Errors.RegulationParsingError import RegulationParsingError
 
 from eBCSgen.Regulations.Base import BaseRegulation
 
@@ -23,3 +24,18 @@ def filter(self, current_state, candidates):
         path = "".join(current_state.memory.history)
         return {rule: values for rule, values in candidates.items()
                 if self.regulation.fullmatch(path + rule.label, partial=True) is not None}
+    
+    def check_labels(self, model_labels):
+        patterns = self.regulation.pattern.replace("(", "").replace(")", "").split("|")
+        subpaterns_set = set()
+        for pattern in patterns:
+            subpaterns_set = subpaterns_set.union(set(pattern.split(";")))
+
+        for subpattern in subpaterns_set:
+            subregex = regex.compile(subpattern)
+            if any(subregex.search(label) for label in model_labels):
+                continue
+            raise RegulationParsingError(
+                f"Label in programmed regulation not present in model"
+            )
+        return True
diff --git a/eBCSgen/utils.py b/eBCSgen/utils.py
new file mode 100644
index 0000000..11e6001
--- /dev/null
+++ b/eBCSgen/utils.py
@@ -0,0 +1,14 @@
+from lark import Tree
+
+
+def tree_to_string(tree):
+    """
+    Recursively constructs a list form given lark tree.
+
+    :param tree: given lark tree
+    :return: list of components
+    """
+    if type(tree) == Tree:
+        return sum(list(map(tree_to_string, tree.children)), [])
+    else:
+        return [str(tree)]
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..d86ec0b
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,13 @@
+numpy
+pandas
+pyparsing
+regex
+requests
+scipy
+sortedcontainers
+sympy
+python-libsbml
+lark
+lark-parser
+pytest
+pymodelchecking
\ No newline at end of file