diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
deleted file mode 100644
index 551a37799b..0000000000
--- a/.github/workflows/main.yml
+++ /dev/null
@@ -1,17 +0,0 @@
-on: [push, pull_request]
-
-jobs:
-  run_tests:
-    runs-on: ubuntu-20.04
-    name: StochSS Continuous Testing
-    steps:
-      # Checkout
-    - name: Checkout
-      id: checkout
-      uses: actions/checkout@v2.1.0
-      # Setup container & execute Python unit tests
-    - name: run_tests.py
-      id: run_tests
-      uses: ./
-      with:
-        testin: '/stochss/stochss/tests/run_tests.py'
diff --git a/.github/workflows/run-integration-tests.yaml b/.github/workflows/run-integration-tests.yaml
new file mode 100644
index 0000000000..e911d3a908
--- /dev/null
+++ b/.github/workflows/run-integration-tests.yaml
@@ -0,0 +1,23 @@
+name: Run StochSS Integration Tests
+
+on:
+  push:
+    branches: [main, develop]
+
+jobs:
+  run-tests:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ ubuntu-20.04 ]
+
+    steps:
+    - name: Initialize environment
+      id: checkout
+      uses: actions/checkout@v2
+
+    - name: Run Integration Tests
+      id: run_integration_tests
+      uses: ./
+      with:
+        testin: '/stochss/stochss/tests/run_integration_tests.py'
diff --git a/.github/workflows/run-unit-tests.yml b/.github/workflows/run-unit-tests.yml
new file mode 100644
index 0000000000..6168fce283
--- /dev/null
+++ b/.github/workflows/run-unit-tests.yml
@@ -0,0 +1,21 @@
+name: Run StochSS Unit Tests
+
+on: [push]
+
+jobs:
+  run-tests:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ ubuntu-20.04 ]
+
+    steps:
+    - name: Initialize environment
+      id: checkout
+      uses: actions/checkout@v2
+
+    - name: Run Unit Tests
+      id: run_unit_tests
+      uses: ./
+      with:
+        testin: '/stochss/stochss/tests/run_unit_tests.py'
diff --git a/Makefile b/Makefile
index a4f706ddc8..3df4c94335 100644
--- a/Makefile
+++ b/Makefile
@@ -148,7 +148,7 @@ test: create_working_dir
 		-v $(DOCKER_WORKING_DIR):/home/jovyan/ \
 		-p 8888:8888 \
 		$(DOCKER_STOCHSS_IMAGE):latest \
-		/stochss/stochss/tests/run_tests.py
+		/stochss/stochss/tests/run_all_tests.py
 
 build_and_test: build test
 
diff --git a/__version__.py b/__version__.py
index d3476703d8..445149aff1 100644
--- a/__version__.py
+++ b/__version__.py
@@ -5,7 +5,7 @@
 # @website https://github.com/stochss/stochss
 # =============================================================================
 
-__version__ = '2.5.0'
+__version__ = '2.5.1'
 __title__ = 'StochSS'
 __description__ = 'StochSS is an integrated development environment (IDE) \
                     for simulation of biochemical networks.'
diff --git a/jupyterhub/Dockerfile.jupyterhub b/jupyterhub/Dockerfile.jupyterhub
index fd28f24821..d959eaaea7 100644
--- a/jupyterhub/Dockerfile.jupyterhub
+++ b/jupyterhub/Dockerfile.jupyterhub
@@ -40,14 +40,14 @@ EXPOSE 8000
 
 # Install dockerspawner, oauth, postgres
 RUN python3 -m pip install --no-cache-dir \
-        jupyterhub \
+        jupyterhub==3.1.1 \
         oauthenticator==0.11.* \
         dockerspawner==0.11.* \
         psycopg2-binary==2.9.* \
         nbviewer==1.0.1 \
         notebook \
-        gillespy2==1.7.2 \
-        spatialpy==1.1.2 \
+        gillespy2==1.8.1 \
+        spatialpy==1.2.1 \
         plotly
 
 #  build the PUG pages
diff --git a/requirements.txt b/requirements.txt
index 2488d89577..671f6e35f7 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -3,7 +3,7 @@ python-libsedml==2.0.9
 python-libcombine==0.2.7
 pyabc==0.12.6
 escapism==1.0.1
-gillespy2==1.8.0
+gillespy2==1.8.1
 spatialpy==1.2.1
 stochss-compute[AWS]==1.0.1
 git+https://github.com/StochSS/sciope.git@master
diff --git a/stochss/tests/__init__.py b/stochss/tests/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/stochss/tests/example_models.py b/stochss/tests/example_models.py
index 267a3c2db7..7d2d9556c5 100644
--- a/stochss/tests/example_models.py
+++ b/stochss/tests/example_models.py
@@ -22,329 +22,437 @@
     Model,
     Species,
     Reaction,
-    Parameter
+    Parameter,
+    TimeSpan
 )
 
 # pylint: disable=line-too-long
 # pylint: disable=missing-class-docstring
 # pylint: disable=too-few-public-methods
-class Brusselator(Model):
-    def __init__(self):
-        Model.__init__(self, name="Brusselator")
-        self.volume = 1000
-
-        # Parameters
-        self.add_parameter(Parameter(name="rate1", expression="5000"))
-        self.add_parameter(Parameter(name="rate2", expression="50"))
-        self.add_parameter(Parameter(name="rate3", expression="5e-05"))
-        self.add_parameter(Parameter(name="rate4", expression="5"))
-
-        # Variables
-        self.add_species(Species(name="A", initial_value=100000, mode="discrete"))
-        self.add_species(Species(name="B", initial_value=100000, mode="discrete"))
-        self.add_species(Species(name="C", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="D", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="X", initial_value=2000, mode="discrete"))
-        self.add_species(Species(name="Y", initial_value=1000, mode="discrete"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="reaction1", reactants={'A': 1}, products={'X': 1, 'A': 1}, propensity_function="rate1"))
-        self.add_reaction(Reaction(name="reaction2", reactants={'B': 1, 'X': 1}, products={'Y': 1, 'C': 1, 'B': 1}, propensity_function="rate2*X"))
-        self.add_reaction(Reaction(name="reaction3", reactants={'X': 2, 'Y': 1}, products={'X': 3}, propensity_function="rate3*Y*X*(X-1)/2"))
-        self.add_reaction(Reaction(name="reaction4", reactants={'X': 1}, products={'D': 1}, propensity_function="rate4*X"))
-
-        # Timespan
-        self.timespan(np.arange(0, 30, 0.01))
-
-
-class Degradation(Model):
-    def __init__(self):
-        Model.__init__(self, name="Degradation")
-        self.volume = 1
-
-        # Parameters
-        self.add_parameter(Parameter(name="decayrate", expression="0.05"))
-
-        # Variables
-        self.add_species(Species(name="protein", initial_value=50, mode="discrete"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="reaction", reactants={'protein': 1}, products={}, rate=self.listOfParameters["decayrate"]))
-
-        # Timespan
-        self.timespan(np.arange(0, 100, 1))
-
-
-class Dimerization(Model):
-    def __init__(self):
-        Model.__init__(self, name="Dimerization")
-        self.volume = 1
-
-        # Parameters
-        self.add_parameter(Parameter(name="k_c", expression="0.005"))
-        self.add_parameter(Parameter(name="k_d", expression="0.08"))
-
-        # Variables
-        self.add_species(Species(name="Monomer", initial_value=30, mode="discrete"))
-        self.add_species(Species(name="Dimer", initial_value=0, mode="discrete"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="r_creation", reactants={'Monomer': 2}, products={'Dimer': 1}, rate=self.listOfParameters["k_c"]))
-        self.add_reaction(Reaction(name="r_dissociation", reactants={'Dimer': 1}, products={'Monomer': 2}, rate=self.listOfParameters["k_d"]))
-
-        # Timespan
-        self.timespan(np.arange(0, 100, 1))
-
-
-class LotkavolterraOscillator(Model):
-    def __init__(self):
-        Model.__init__(self, name="Lotkavolterra_Oscillator")
-        self.volume = 1
-
-        # Parameters
-        self.add_parameter(Parameter(name="k1", expression="0.004"))
-        self.add_parameter(Parameter(name="k2", expression="0.5"))
-        self.add_parameter(Parameter(name="k3", expression="0.0045"))
-
-        # Variables
-        self.add_species(Species(name="R", initial_value=0.0095, mode="continuous"))
-        self.add_species(Species(name="W", initial_value=0.008, mode="continuous"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="r1", reactants={'R': 1}, products={'R': 2}, rate=self.listOfParameters["k1"]))
-        self.add_reaction(Reaction(name="r2", reactants={'R': 1, 'W': 1}, products={'W': 2}, rate=self.listOfParameters["k2"]))
-        self.add_reaction(Reaction(name="r3", reactants={'W': 1}, products={}, rate=self.listOfParameters["k3"]))
-
-        # Timespan
-        self.timespan(np.arange(0, 9000, 1))
-
-
-class MichaelisMenten(Model):
-    def __init__(self):
-        Model.__init__(self, name="Michaelis_Menten")
-        self.volume = 1
-
-        # Parameters
-        self.add_parameter(Parameter(name="rate1", expression="0.0017"))
-        self.add_parameter(Parameter(name="rate2", expression="0.5"))
-        self.add_parameter(Parameter(name="rate3", expression="0.1"))
-
-        # Variables
-        self.add_species(Species(name="A", initial_value=301, mode="discrete"))
-        self.add_species(Species(name="B", initial_value=120, mode="discrete"))
-        self.add_species(Species(name="C", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="D", initial_value=0, mode="discrete"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="r1", reactants={'A': 1, 'B': 1}, products={'C': 1}, rate=self.listOfParameters["rate1"]))
-        self.add_reaction(Reaction(name="r2", reactants={'C': 1}, products={'A': 1, 'B': 1}, rate=self.listOfParameters["rate2"]))
-        self.add_reaction(Reaction(name="r3", reactants={'C': 1}, products={'B': 1, 'D': 1}, rate=self.listOfParameters["rate3"]))
-
-        # Timespan
-        self.timespan(np.arange(0, 100, 1))
-
-
-class Opioid(Model):
-    def __init__(self):
-        Model.__init__(self, name="Opioid")
-        self.volume = 1
-
-        # Parameters
-        self.add_parameter(Parameter(name="alpha", expression="0.15"))
-        self.add_parameter(Parameter(name="epsilon", expression="0.8"))
-        self.add_parameter(Parameter(name="beta_p", expression="0.00266"))
-        self.add_parameter(Parameter(name="beta_a", expression="0.00094"))
-        self.add_parameter(Parameter(name="gamma", expression="0.00744"))
-        self.add_parameter(Parameter(name="zeta", expression="0.2"))
-        self.add_parameter(Parameter(name="delta", expression="0.1"))
-        self.add_parameter(Parameter(name="sigma", expression="0.9"))
-        self.add_parameter(Parameter(name="mu", expression="0.00729"))
-        self.add_parameter(Parameter(name="mu_prime", expression="0.01159"))
-
-        # Variables
-        self.add_species(Species(name="Susceptible", initial_value=200, mode="discrete"))
-        self.add_species(Species(name="Prescribed_User", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="Addicted", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="Rehab", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="Natural_Deaths", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="Addicted_Deaths", initial_value=0, mode="discrete"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="SP", reactants={'Susceptible': 1}, products={'Prescribed_User': 1}, rate=self.listOfParameters["alpha"]))
-        self.add_reaction(Reaction(name="SA_a", reactants={'Susceptible': 1}, products={'Addicted': 1}, rate=self.listOfParameters["beta_a"]))
-        self.add_reaction(Reaction(name="SA_p", reactants={'Susceptible': 1}, products={'Addicted': 1}, rate=self.listOfParameters["beta_p"]))
-        self.add_reaction(Reaction(name="PA", reactants={'Prescribed_User': 1}, products={'Addicted': 1}, rate=self.listOfParameters["gamma"]))
-        self.add_reaction(Reaction(name="PS", reactants={'Prescribed_User': 1}, products={'Susceptible': 1}, rate=self.listOfParameters["epsilon"]))
-        self.add_reaction(Reaction(name="AR", reactants={'Addicted': 1}, products={'Rehab': 1}, rate=self.listOfParameters["zeta"]))
-        self.add_reaction(Reaction(name="RA", reactants={'Rehab': 1}, products={'Addicted': 1}, rate=self.listOfParameters["delta"]))
-        self.add_reaction(Reaction(name="RS", reactants={'Rehab': 1}, products={'Susceptible': 1}, rate=self.listOfParameters["sigma"]))
-        self.add_reaction(Reaction(name="mu_S", reactants={'Susceptible': 1}, products={'Susceptible': 1, 'Natural_Deaths': 1}, rate=self.listOfParameters["mu"]))
-        self.add_reaction(Reaction(name="mu_P", reactants={'Prescribed_User': 1}, products={'Susceptible': 1, 'Natural_Deaths': 1}, rate=self.listOfParameters["mu"]))
-        self.add_reaction(Reaction(name="mu_R", reactants={'Rehab': 1}, products={'Susceptible': 1, 'Natural_Deaths': 1}, rate=self.listOfParameters["mu"]))
-        self.add_reaction(Reaction(name="mu_prime_A", reactants={'Addicted': 1}, products={'Susceptible': 1, 'Addicted_Deaths': 1}, rate=self.listOfParameters["mu_prime"]))
-
-        # Timespan
-        self.timespan(np.arange(0, 200, 1))
-
-
-class Schlogl(Model):
-    def __init__(self):
-        Model.__init__(self, name="Schlogl")
-        self.volume = 1
-
-        # Parameters
-        self.add_parameter(Parameter(name="k1", expression="1"))
-        self.add_parameter(Parameter(name="k2", expression="1"))
-
-        # Variables
-        self.add_species(Species(name="s1", initial_value=300, mode="discrete"))
-        self.add_species(Species(name="s2", initial_value=300, mode="discrete"))
-        self.add_species(Species(name="s3", initial_value=300, mode="discrete"))
-        self.add_species(Species(name="s4", initial_value=300, mode="discrete"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="r1", reactants={'s1': 1, 's4': 1}, products={'s4': 2}, rate=self.listOfParameters["k1"]))
-        self.add_reaction(Reaction(name="r2", reactants={'s2': 1, 's4': 1}, products={'s3': 1}, rate=self.listOfParameters["k2"]))
-
-        # Timespan
-        self.timespan(np.arange(0, 100000, 1000))
-
-
-class ToggleSwitch(Model):
-    def __init__(self):
-        Model.__init__(self, name="Toggle_Switch")
-        self.volume = 1
-
-        # Parameters
-        self.add_parameter(Parameter(name="alpha1", expression="10"))
-        self.add_parameter(Parameter(name="alpha2", expression="10"))
-        self.add_parameter(Parameter(name="beta", expression="2"))
-        self.add_parameter(Parameter(name="gamma", expression="2"))
-        self.add_parameter(Parameter(name="mu", expression="1"))
-
-        # Variables
-        self.add_species(Species(name="A", initial_value=2, mode="discrete"))
-        self.add_species(Species(name="B", initial_value=2, mode="discrete"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="cu", reactants={}, products={'A': 1}, propensity_function="alpha1/(1+pow(B, beta))"))
-        self.add_reaction(Reaction(name="cv", reactants={}, products={'B': 1}, propensity_function="alpha2/(1+pow(A, gamma))"))
-        self.add_reaction(Reaction(name="du", reactants={'A': 1}, products={}, rate=self.listOfParameters["mu"]))
-        self.add_reaction(Reaction(name="dv", reactants={'B': 1}, products={}, rate=self.listOfParameters["mu"]))
-
-        # Timespan
-        self.timespan(np.arange(0, 250, 1))
-
-
-class VilarOscillator(Model):
-    def __init__(self):
-        Model.__init__(self, name="Vilar_Oscillator")
-        self.volume = 1
-
-        # Parameters
-        self.add_parameter(Parameter(name="alphaA", expression="50"))
-        self.add_parameter(Parameter(name="alphaA_prime", expression="500"))
-        self.add_parameter(Parameter(name="alphaR", expression="0.01"))
-        self.add_parameter(Parameter(name="alphaR_prime", expression="50"))
-        self.add_parameter(Parameter(name="betaA", expression="50"))
-        self.add_parameter(Parameter(name="betaR", expression="5"))
-        self.add_parameter(Parameter(name="deltaMA", expression="10"))
-        self.add_parameter(Parameter(name="deltaMR", expression="0.5"))
-        self.add_parameter(Parameter(name="deltaA", expression="1"))
-        self.add_parameter(Parameter(name="deltaR", expression="0.2"))
-        self.add_parameter(Parameter(name="gammaA", expression="1"))
-        self.add_parameter(Parameter(name="gammaR", expression="1"))
-        self.add_parameter(Parameter(name="gammaC", expression="2"))
-        self.add_parameter(Parameter(name="thetaA", expression="50"))
-        self.add_parameter(Parameter(name="thetaR", expression="100"))
-
-        # Variables
-        self.add_species(Species(name="Da", initial_value=1, mode="discrete"))
-        self.add_species(Species(name="Da_prime", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="Ma", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="Dr", initial_value=1, mode="discrete"))
-        self.add_species(Species(name="Dr_prime", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="Mr", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="C", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="A", initial_value=0, mode="discrete"))
-        self.add_species(Species(name="R", initial_value=0, mode="discrete"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="r1", reactants={'A': 1, 'R': 1}, products={'C': 1}, rate="gammaC"))
-        self.add_reaction(Reaction(name="r2", reactants={'A': 1}, products={}, rate="deltaA"))
-        self.add_reaction(Reaction(name="r3", reactants={'C': 1}, products={'R': 1}, rate="deltaA"))
-        self.add_reaction(Reaction(name="r4", reactants={'R': 1}, products={}, rate="deltaR"))
-        self.add_reaction(Reaction(name="r5", reactants={'A': 1, 'Da': 1}, products={'Da_prime': 1}, rate="gammaA"))
-        self.add_reaction(Reaction(name="r6", reactants={'Da_prime': 1}, products={'A': 1, 'Da': 1}, rate="thetaA"))
-        self.add_reaction(Reaction(name="r7", reactants={'Da': 1}, products={'Da': 1, 'Ma': 1}, rate="alphaA"))
-        self.add_reaction(Reaction(name="r8", reactants={'Da_prime': 1}, products={'Da_prime': 1, 'Ma': 1}, rate="alphaA_prime"))
-        self.add_reaction(Reaction(name="r9", reactants={'Ma': 1}, products={}, rate="deltaMA"))
-        self.add_reaction(Reaction(name="r10", reactants={'Ma': 1}, products={'A': 1, 'Ma': 1}, rate="betaA"))
-        self.add_reaction(Reaction(name="r11", reactants={'A': 1, 'Dr': 1}, products={'Dr_prime': 1}, rate="gammaR"))
-        self.add_reaction(Reaction(name="r12", reactants={'Dr_prime': 1}, products={'A': 1, 'Dr': 1}, rate="thetaR"))
-        self.add_reaction(Reaction(name="r13", reactants={'Dr': 1}, products={'Dr': 1, 'Mr': 1}, rate="alphaR"))
-        self.add_reaction(Reaction(name="r14", reactants={'Dr_prime': 1}, products={'Dr_prime': 1, 'Mr': 1}, rate="alphaR_prime"))
-        self.add_reaction(Reaction(name="r15", reactants={'Mr': 1}, products={}, rate="deltaMR"))
-        self.add_reaction(Reaction(name="r16", reactants={'Mr': 1}, products={'Mr': 1, 'R': 1}, rate="betaR"))
-
-        # Timespan
-        self.timespan(np.arange(0, 201, 1))
-
-
-class Oregonator(Model):
-    def __init__(self):
-        Model.__init__(self, name="Oregonator")
-        self.volume = 1
-
-        # Parameters
-        self.add_parameter(Parameter(name="k1", expression="2"))
-        self.add_parameter(Parameter(name="k2", expression="0.1"))
-        self.add_parameter(Parameter(name="k3", expression="104"))
-        self.add_parameter(Parameter(name="k4", expression="4e-07"))
-        self.add_parameter(Parameter(name="k5", expression="26"))
-
-        # Variables
-        self.add_species(Species(name="F", initial_value=2, mode="continuous"))
-        self.add_species(Species(name="A", initial_value=250, mode="continuous"))
-        self.add_species(Species(name="B", initial_value=500, mode="continuous"))
-        self.add_species(Species(name="C", initial_value=1000, mode="continuous"))
-        self.add_species(Species(name="P", initial_value=0, mode="continuous"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="reaction1", reactants={'B': 1, 'F': 1}, products={'A': 1, 'F': 1}, rate=self.listOfParameters["k1"]))
-        self.add_reaction(Reaction(name="reaction2", reactants={'A': 1, 'B': 1}, products={'P': 1}, rate=self.listOfParameters["k2"]))
-        self.add_reaction(Reaction(name="reaction3", reactants={'A': 1, 'F': 1}, products={'A': 2, 'C': 1, 'F': 1}, rate=self.listOfParameters["k3"]))
-        self.add_reaction(Reaction(name="reaction4", reactants={'A': 2}, products={'P': 1}, rate=self.listOfParameters["k4"]))
-        self.add_reaction(Reaction(name="reaction5", reactants={'C': 1, 'F': 1}, products={'B': 1, 'F': 1}, rate=self.listOfParameters["k5"]))
-
-        # Timespan
-        self.timespan(np.arange(0, 5, 0.1))
-
-
-class TysonOscillator(Model):
-    def __init__(self):
-        Model.__init__(self, name="Tyson_Oscillator")
-        self.volume = 300
-
-        # Parameters
-        self.add_parameter(Parameter(name="P", expression="2"))
-        self.add_parameter(Parameter(name="kt", expression="20"))
-        self.add_parameter(Parameter(name="kd", expression="1"))
-        self.add_parameter(Parameter(name="a0", expression="0.005"))
-        self.add_parameter(Parameter(name="a1", expression="0.05"))
-        self.add_parameter(Parameter(name="a2", expression="0.1"))
-        self.add_parameter(Parameter(name="kdx", expression="1"))
-
-        # Variables
-        self.add_species(Species(name="X", initial_value=197, mode="discrete"))
-        self.add_species(Species(name="Y", initial_value=255, mode="discrete"))
-
-        # Reactions
-        self.add_reaction(Reaction(name="rxn1", reactants={}, products={'X': 1}, propensity_function="vol*1/(1+(Y*Y/((vol*vol))))"))
-        self.add_reaction(Reaction(name="rxn2", reactants={'X': 1}, products={}, rate=self.listOfParameters["kdx"]))
-        self.add_reaction(Reaction(name="rxn3", reactants={'X': 1}, products={'X': 1, 'Y': 1}, rate=self.listOfParameters["kt"]))
-        self.add_reaction(Reaction(name="rxn4", reactants={'Y': 1}, products={}, rate=self.listOfParameters["kd"]))
-        self.add_reaction(Reaction(name="rxn5", reactants={'Y': 1}, products={}, propensity_function="Y/(a0 + a1*(Y/vol)+a2*Y*Y/(vol*vol))"))
-
-        # Timespan
-        self.timespan(np.arange(0, 100, 1))
+def create_vilar_oscillator(parameter_values=None):
+    model = Model(name="VilarOscillator")
+    model.volume = 1
+
+    # Parameters
+    alphaA = Parameter(name="alphaA", expression=50)
+    alphaA_prime = Parameter(name="alphaA_prime", expression=500)
+    alphaR = Parameter(name="alphaR", expression=0.01)
+    alphaR_prime = Parameter(name="alphaR_prime", expression=50)
+    betaA = Parameter(name="betaA", expression=50)
+    betaR = Parameter(name="betaR", expression=5)
+    deltaMA = Parameter(name="deltaMA", expression=10)
+    deltaMR = Parameter(name="deltaMR", expression=0.5)
+    deltaA = Parameter(name="deltaA", expression=1)
+    deltaR = Parameter(name="deltaR", expression=0.2)
+    gammaA = Parameter(name="gammaA", expression=1)
+    gammaR = Parameter(name="gammaR", expression=1)
+    gammaC = Parameter(name="gammaC", expression=2)
+    thetaA = Parameter(name="thetaA", expression=50)
+    thetaR = Parameter(name="thetaR", expression=100)
+    model.add_parameter([
+        alphaA, alphaA_prime, alphaR, alphaR_prime, betaA, betaR, deltaMA, deltaMR, deltaA, deltaR, gammaA, gammaR, gammaC, thetaA, thetaR
+    ])
+
+    # Species
+    Da = Species(name="Da", initial_value=1, mode="discrete")
+    Da_prime = Species(name="Da_prime", initial_value=0, mode="discrete")
+    Ma = Species(name="Ma", initial_value=0, mode="discrete")
+    Dr = Species(name="Dr", initial_value=1, mode="discrete")
+    Dr_prime = Species(name="Dr_prime", initial_value=0, mode="discrete")
+    Mr = Species(name="Mr", initial_value=0, mode="discrete")
+    C = Species(name="C", initial_value=0, mode="discrete")
+    A = Species(name="A", initial_value=0, mode="discrete")
+    R = Species(name="R", initial_value=0, mode="discrete")
+    model.add_species([Da, Da_prime, Ma, Dr, Dr_prime, Mr, C, A, R])
+
+    # Reactions
+    r01 = Reaction(name="r1", reactants={'A': 1, 'R': 1}, products={'C': 1}, rate="gammaC")
+    r02 = Reaction(name="r2", reactants={'A': 1}, products={}, rate="deltaA")
+    r03 = Reaction(name="r3", reactants={'C': 1}, products={'R': 1}, rate="deltaA")
+    r04 = Reaction(name="r4", reactants={'R': 1}, products={}, rate="deltaR")
+    r05 = Reaction(name="r5", reactants={'A': 1, 'Da': 1}, products={'Da_prime': 1}, rate="gammaA")
+    r06 = Reaction(name="r6", reactants={'Da_prime': 1}, products={'A': 1, 'Da': 1}, rate="thetaA")
+    r07 = Reaction(name="r7", reactants={'Da': 1}, products={'Da': 1, 'Ma': 1}, rate="alphaA")
+    r08 = Reaction(name="r8", reactants={'Da_prime': 1}, products={'Da_prime': 1, 'Ma': 1}, rate="alphaA_prime")
+    r09 = Reaction(name="r9", reactants={'Ma': 1}, products={}, rate="deltaMA")
+    r10 = Reaction(name="r10", reactants={'Ma': 1}, products={'A': 1, 'Ma': 1}, rate="betaA")
+    r11 = Reaction(name="r11", reactants={'A': 1, 'Dr': 1}, products={'Dr_prime': 1}, rate="gammaR")
+    r12 = Reaction(name="r12", reactants={'Dr_prime': 1}, products={'A': 1, 'Dr': 1}, rate="thetaR")
+    r13 = Reaction(name="r13", reactants={'Dr': 1}, products={'Dr': 1, 'Mr': 1}, rate="alphaR")
+    r14 = Reaction(name="r14", reactants={'Dr_prime': 1}, products={'Dr_prime': 1, 'Mr': 1}, rate="alphaR_prime")
+    r15 = Reaction(name="r15", reactants={'Mr': 1}, products={}, rate="deltaMR")
+    r16 = Reaction(name="r16", reactants={'Mr': 1}, products={'Mr': 1, 'R': 1}, rate="betaR")
+    model.add_reaction([r01, r02, r03, r04, r05, r06, r07, r08, r09, r10, r11, r12, r13, r14, r15, r16])
+
+    # Timespan
+    model.timespan(np.linspace(0,200,201))
+    return model
+
+def create_dimerization(parameter_values=None):
+    model = Model(name="Dimerization")
+
+    # Parameters
+    k_c = Parameter(name='k_c', expression=0.005)
+    k_d = Parameter(name='k_d', expression=0.08)
+    model.add_parameter([k_c, k_d])
+
+    # Species
+    m = Species(name='monomer', initial_value=30)
+    d = Species(name='dimer',   initial_value=0)
+    model.add_species([m, d])
+
+    # Reactions
+    r_creation = Reaction(name="r_creation", rate=k_c, reactants={m:2}, products={d:1})
+    r_dissociation = Reaction(name="r_dissociation", rate=k_d, reactants={d:1}, products={m:2})
+    model.add_reaction([r_creation, r_dissociation])
+
+    # Timespan
+    model.timespan(np.linspace(0, 100, 101))
+    return model
+
+def create_trichloroethylene(parameter_values=None):
+    model = Model(name="Trichloroethylene")
+
+    # Species
+    A = Species(name='TCE', initial_value=300)
+    B = Species(name='Epoxide', initial_value=120)
+    C = Species(name='Dichloracatate', initial_value=0)
+    D = Species(name='LossOfOneCL', initial_value=0)
+    E = Species(name='Glyoxylate', initial_value=0)
+    F = Species(name='Output', initial_value=0)
+    model.add_species([A, B, C, D, E, F])
+
+    # Parameters
+    K1 = Parameter(name='K1', expression=0.00045 * 0.000025)
+    K2 = Parameter(name='K2', expression=14)
+    K3 = Parameter(name='K3', expression=0.033)
+    K4 = Parameter(name='K4', expression=500 * 0.0001)
+    K5 = Parameter(name='K5', expression=500 * 0.0001)
+    model.add_parameter([K1, K2, K3, K4, K5])
+
+    # Reactions
+    J1 = Reaction(name="J1", reactants={A: 1}, products={B: 1}, rate=K2)
+    J2 = Reaction(name="J2", reactants={B: 1}, products={C: 1}, rate=K3)
+    J3 = Reaction(name="J3", reactants={C: 1}, products={D: 1}, rate=K4)
+    J4 = Reaction(name="J4", reactants={D: 1}, products={E: 1}, rate=K4)
+    J5 = Reaction(name="J5", reactants={E: 1}, products={F: 1}, rate=K5)
+    model.add_reaction([J1, J2, J3, J4, J5])
+
+    model.timespan(np.linspace(0, 10000, 100))
+    return model
+
+def create_lac_operon(parameter_values=None):
+    model = Model(name="LacOperon")
+
+    # Species
+    s1 = Species(name='MR', initial_value=0)
+    s2 = Species(name='R', initial_value=0)
+    s3 = Species(name='R2', initial_value=0)
+    s4 = Species(name='O', initial_value=1)
+    s5 = Species(name='R2O', initial_value=0)
+    s6 = Species(name='I', initial_value=0)
+    s7 = Species(name='Iex', initial_value=0)
+    s8 = Species(name='I2R2', initial_value=0)
+    s9 = Species(name='MY', initial_value=0)
+    s10 = Species(name='Y', initial_value=0)
+    s11 = Species(name='YIex', initial_value=0)
+    s12 = Species(name='Ytot', initial_value=0)
+    model.add_species([s1, s2, s3, s4, s5, s6, s7, s8, s9, s10, s11, s12])
+
+    # Parameters
+    k1 = Parameter(name='k1', expression=0.111)
+    k2 = Parameter(name='k2', expression=15.0)
+    k3 = Parameter(name='k3', expression=103.8)
+    k4 = Parameter(name='k4', expression=0.001)
+    k5 = Parameter(name='k5', expression=1992.7)
+    k6 = Parameter(name='k6', expression=2.40)
+    k7 = Parameter(name='k7', expression=1.293e-6)
+    k8 = Parameter(name='k8', expression=12.0)
+    k9 = Parameter(name='k9', expression=1.293e-6)
+    k10 = Parameter(name='k10', expression=9963.2)
+    k11 = Parameter(name='k11', expression=0.50)
+    k12 = Parameter(name='k12', expression=0.010)
+    k13 = Parameter(name='k13', expression=30)
+    k14 = Parameter(name='k14', expression=0.249)
+    k15 = Parameter(name='k15', expression=0.10)
+    k16 = Parameter(name='k16', expression=60000)
+    k17 = Parameter(name='k17', expression=0.920)
+    k18 = Parameter(name='k18', expression=0.920)
+    k19 = Parameter(name='k19', expression=0.462)
+    k20 = Parameter(name='k20', expression=0.462)
+    k21 = Parameter(name='k21', expression=0.20)
+    k22 = Parameter(name='k22', expression=0.20)
+    k23 = Parameter(name='k23', expression=0.20)
+    k24 = Parameter(name='k24', expression=0.20)
+    k25 = Parameter(name='k25', expression=0.20)
+    model.add_parameter(
+        [k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15, k16, k17, k18, k19, k20, k21, k22, k23, k24, k25
+    ])
+
+    # Reactions
+    j1 = Reaction(name="j1", reactants={s12: 1}, products={s1: 1}, rate=k1)
+    j2 = Reaction(name="j2", reactants={s1: 1}, products={s1: 1, s2: 1}, rate=k2)
+    j3 = Reaction(name="j3", reactants={s2: 2}, products={s3: 1}, rate=k3)
+    j4 = Reaction(name="j4", reactants={s3: 1}, products={s2: 2}, rate=k4)
+    j5 = Reaction(name="j5", reactants={s3: 1, s4: 1}, products={s5: 1}, rate=k5)
+    j6 = Reaction(name="j6", reactants={s5: 1}, products={s3: 1, s4: 1}, rate=k6)
+    j7 = Reaction(name="j7", reactants={s6: 2, s3: 1}, products={s8: 1}, propensity_function='((3e-7)/((8e-16)*(6.0221367e14))**(2))*(R2)*(I)*(I-1)')
+    j8 = Reaction(name="j8", reactants={s8: 1}, products={s6: 2, s3: 1}, rate=k8)
+    j9 = Reaction(name="j9", reactants={s6: 2, s5: 1}, products={s8: 1, s4: 1}, propensity_function='((3e-7)/((8e-16)*(6.0221367e14))**(2))*(R2O)*(I)*(I-1)')
+    j10 = Reaction(name="j10", reactants={s8: 1, s4: 1}, products={s6: 2, s5: 1}, rate=k10)
+    j11 = Reaction(name="j11", reactants={s4: 1}, products={s4: 1, s9: 1}, rate=k11)
+    j12 = Reaction(name="j12", reactants={s5: 1}, products={s5: 1, s9: 1}, rate=k12)
+    j13 = Reaction(name="j13", reactants={s9: 1}, products={s9: 1, s10: 1}, rate=k13)
+    j14 = Reaction(name="j14", reactants={s10: 1, s7: 1}, products={s11: 1}, rate=k14)
+    j15 = Reaction(name="j15", reactants={s11: 1}, products={s10: 1, s7: 1}, rate=k15)
+    j16 = Reaction(name="j16", reactants={s11: 1}, products={s10: 1, s6: 1}, rate=k16)
+    j17 = Reaction(name="j17", reactants={s7: 1}, products={s6: 1}, rate=k17)
+    j18 = Reaction(name="j18", reactants={s6: 1}, products={s7: 1}, rate=k18)
+    j19 = Reaction(name="j19", reactants={s1: 1}, products={s12: 1}, rate=k19)
+    j20 = Reaction(name="j20", reactants={s9: 1}, products={s12: 1}, rate=k20)
+    j21 = Reaction(name="j21", reactants={s2: 1}, products={s12: 1}, rate=k21)
+    j22 = Reaction(name="j22", reactants={s3: 1}, products={s12: 1}, rate=k22)
+    j23 = Reaction(name="j23", reactants={s10: 1}, products={s12: 1}, rate=k23)
+    j24 = Reaction(name="j24", reactants={s11: 1}, products={s6: 1}, rate=k24)
+    j25 = Reaction(name="j25", reactants={s8: 1}, products={s6: 2}, rate=k25)
+    model.add_reaction(
+        [j1, j2, j3, j4, j5, j6, j7, j8, j9, j10, j11, j12, j13, j14, j15, j16, j17, j18, j19, j20, j21, j22, j23, j24, j25
+    ])
+
+    model.timespan(np.linspace(0, 100, 10))
+    return model
+
+def create_schlogl(parameter_values=None):
+    model = Model(name="Schlogl")
+
+    # Species
+    s1 = Species(name='A', initial_value=300)
+    s2 = Species(name='B', initial_value=300)
+    s3 = Species(name='C', initial_value=300)
+    s4 = Species(name='X', initial_value=300)
+    model.add_species([s1, s2, s3, s4])
+
+    k1 = Parameter(name='k1', expression=1)
+    k2 = Parameter(name='k2', expression=1)
+    model.add_parameter([k1, k2])
+
+    j1 = Reaction(name="j1", reactants={s1: 1, s4: 1}, products={s4: 2}, rate=k1)
+    j2 = Reaction(name="j2", reactants={s2: 1, s4: 1}, products={s3: 1}, rate=k2)
+    model.add_reaction([j1, j2])
+
+    model.timespan(np.linspace(0, 100000, 100))
+    return model
+
+def create_michaelis_menten(parameter_values=None):
+    model = Model(name="Michaelis_Menten")
+
+    # parameters
+    rate1 = Parameter(name='rate1', expression=0.0017)
+    rate2 = Parameter(name='rate2', expression=0.5)
+    rate3 = Parameter(name='rate3', expression=0.1)
+    model.add_parameter([rate1, rate2, rate3])
+
+    # Species
+    A = Species(name='A', initial_value=301)
+    B = Species(name='B', initial_value=120)
+    C = Species(name='C', initial_value=0)
+    D = Species(name='D', initial_value=0)
+    model.add_species([A, B, C, D])
+
+    # reactions
+    r1 = Reaction(name="r1", reactants={A: 1, B: 1}, products={C: 1}, rate=rate1)
+    r2 = Reaction(name="r2", reactants={C: 1}, products={A: 1, B: 1}, rate=rate2)
+    r3 = Reaction(name="r3", reactants={C: 1}, products={B: 1, D: 1}, rate=rate3)
+    model.add_reaction([r1, r2, r3])
+
+    model.timespan(np.linspace(0, 100, 101))
+    return model
+
+def create_toggle_switch(parameter_values=None):
+    model = Model(name="Toggle_Switch")
+
+    # Species
+    A = Species(name='A', initial_value=10)
+    B = Species(name='B', initial_value=10)
+    model.add_species([A, B])
+
+    # Parameters
+    alpha1 = Parameter(name='alpha1', expression=10)
+    alpha2 = Parameter(name='alpha2', expression=10)
+    beta = Parameter(name='beta', expression=2)
+    gamma = Parameter(name='gamma', expression=2)
+    mu = Parameter(name='mu', expression=1)
+    model.add_parameter([alpha1, alpha2, beta, gamma, mu])
+
+    # Reactions
+    model.add_reaction(Reaction(name="cu", reactants={}, products={'A': 1}, propensity_function="alpha1/(1+pow(B, beta))"))
+    model.add_reaction(Reaction(name="cv", reactants={}, products={'B': 1}, propensity_function="alpha2/(1+pow(A, gamma))"))
+    model.add_reaction(Reaction(name="du", reactants={'A': 1}, products={}, rate=model.listOfParameters["mu"]))
+    model.add_reaction(Reaction(name="dv", reactants={'B': 1}, products={}, rate=model.listOfParameters["mu"]))
+
+    model.timespan(np.linspace(0, 250, 251))
+    return model
+
+def create_decay(parameter_values=None):
+    model = Model(name="Example")
+
+    # Species
+    S = Species(name='Sp', initial_value=100)
+    model.add_species([S])
+
+    # Parameters
+    k1 = Parameter(name='k1', expression=3.0)
+    model.add_parameter([k1])
+
+    # Reactions
+    rxn1 = Reaction(name='S degradation', reactants={S: 1}, products={}, rate=k1)
+    model.add_reaction([rxn1])
+
+    model.timespan(np.linspace(0, 20, 101))
+    return model
+
+def create_tyson_2_state_oscillator(parameter_values=None):
+    model = Model(name="tyson-2-state", volume=300.0)
+
+    # Species
+    X = Species(name='X', initial_value=int(0.65609071 * 300.0))
+    Y = Species(name='Y', initial_value=int(0.85088331 * 300.0))
+    model.add_species([X, Y])
+
+    # Parameters
+    P = Parameter(name='p', expression=2.0)
+    kt = Parameter(name='kt', expression=20.0)
+    kd = Parameter(name='kd', expression=1.0)
+    a0 = Parameter(name='a0', expression=0.005)
+    a1 = Parameter(name='a1', expression=0.05)
+    a2 = Parameter(name='a2', expression=0.1)
+    kdx = Parameter(name='kdx', expression=1.0)
+    model.add_parameter([P, kt, kd, a0, a1, a2, kdx])
+
+    # Reactions
+    rxn1 = Reaction(
+        name='X production', reactants={}, products={X: 1}, propensity_function='300 * 1.0 / (1.0 + (Y * Y / (300 * 300)))'
+    )
+    rxn2 = Reaction(name='X degradation', reactants={X: 1}, products={}, rate=kdx)
+    rxn3 = Reaction(name='Y production', reactants={X: 1}, products={X: 1, Y: 1}, rate=kt)
+    rxn4 = Reaction(name='Y degradation', reactants={Y: 1}, products={}, rate=kd)
+    rxn5 = Reaction(
+        name='Y nonlin', reactants={Y: 1}, products={}, propensity_function='Y / a0 + a1 * (Y / 300) + a2 * Y * Y / (300 * 300)'
+    )
+    model.add_reaction([rxn1, rxn2, rxn3, rxn4, rxn5])
+
+    model.timespan(np.linspace(0, 100, 101))
+    return model
+
+def create_oregonator(parameter_values = None):
+    model = Model(name="Oregonator")
+
+    # Species
+    F = Species(name="F", initial_value=2)
+    A = Species(name="A", initial_value=250)
+    B = Species(name="B", initial_value=500)
+    C = Species(name="C", initial_value=1000)
+    P = Species(name="P", initial_value=0)
+    model.add_species([F, A, B, C, P])
+
+    # Parameters (rates)
+    k1 = Parameter(name="k1", expression=2.0)
+    k2 = Parameter(name="k2", expression=0.1)
+    k3 = Parameter(name="k3", expression=104)
+    k4 = Parameter(name="k4", expression=4e-7)
+    k5 = Parameter(name="k5", expression=26.0)
+    model.add_parameter([k1, k2, k3, k4, k5])
+
+    # Reactions
+    reaction1 = Reaction(name="reaction1", reactants={B: 1, F: 1}, products={A: 1, F: 1}, rate=k1)
+    reaction2 = Reaction(name="reaction2", reactants={A: 1, B: 1}, products={P: 1}, rate=k2)
+    reaction3 = Reaction(name="reaction3", reactants={A: 1, F: 1}, products={A: 2, C: 1, F: 1}, rate=k3)
+    reaction4 = Reaction(name="reaction4", reactants={A: 2}, products={P: 1}, rate=k4)
+    reaction5 = Reaction(name="reaction5", reactants={C: 1, F: 1}, products={B: 1, F: 1}, rate=k5)
+    model.add_reaction([reaction1, reaction2, reaction3, reaction4, reaction5])
+
+    # Set timespan of model
+    model.timespan(np.linspace(0, 5, 501))
+    return model
+
+def create_opioid():
+    model = Model(name="Opioid")
+
+    # Species
+    S = Species(name='Susceptibles', initial_value=1000)
+    P = Species(name='Prescribed_Users', initial_value=0)
+    A = Species(name='Addicted', initial_value=0)
+    R = Species(name='Rehab', initial_value=0)
+    Natural_Deaths = Species(name='Natural_Deaths', initial_value=0)
+    Addiction_Deaths = Species(name='Addiction_Deaths', initial_value=0)
+    model.add_species([S,P,A,R,Natural_Deaths,Addiction_Deaths])
+
+    # Parameters
+    alpha = Parameter(name='alpha', expression= 0.15)
+    epsilon = Parameter(name='epsilon', expression= 0.8)
+    beta_p = Parameter(name='beta_p', expression= 0.00266)
+    beta_a = Parameter(name='beta_a', expression= 0.00094)
+    gamma = Parameter(name='gamma', expression= 0.00744)
+    zeta = Parameter(name='zeta', expression= 0.2)
+    delta = Parameter(name='delta', expression= 0.1)
+    sigma = Parameter(name='sigma', expression= 0.9)
+    mu = Parameter(name='mu', expression= 0.00729)
+    mu_prime = Parameter(name='mu_prime', expression= 0.01159)
+    model.add_parameter([alpha, epsilon, beta_p, beta_a, gamma, zeta, delta, sigma, mu, mu_prime])
+
+    # Reactions
+    SP = Reaction(
+        name="SP", reactants={'Susceptibles': 1}, products={'Prescribed_Users': 1}, rate='alpha'
+    )
+    SA_a = Reaction(name="SA_a", reactants={'Susceptibles': 1}, products={'Addicted': 1}, rate='beta_a')
+    SA_p = Reaction(name="SA_p", reactants={'Susceptibles': 1}, products={'Addicted': 1}, rate='beta_p')
+    mu_S = Reaction(
+        name="mu_S", reactants={'Susceptibles': 1}, products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu'
+    )
+    PA = Reaction(name="PA", reactants={'Prescribed_Users': 1}, products={'Addicted': 1}, rate='gamma')
+    PS = Reaction(
+        name="PS", reactants={'Prescribed_Users': 1}, products={'Susceptibles': 1}, rate='epsilon'
+    )
+    AR = Reaction(name="AR", reactants={'Addicted': 1}, products={'Rehab': 1}, rate='zeta')
+    RA = Reaction(name="RA", reactants={'Rehab': 1}, products={'Addicted': 1}, rate='delta')
+    RS = Reaction(name="RS", reactants={'Rehab': 1}, products={'Susceptibles': 1}, rate='sigma')
+    mu_P = Reaction(
+        name="mu_P", reactants={'Prescribed_Users': 1},
+        products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu'
+    )
+    mu_R = Reaction(
+        name="m_R", reactants={'Rehab': 1}, products={'Susceptibles': 1, 'Natural_Deaths': 1}, rate='mu'
+    )
+    mu_prime_A = Reaction(
+        name="mu_prime_A", reactants={'Addicted': 1},
+        products={'Susceptibles': 1, 'Addiction_Deaths': 1}, rate='mu_prime'
+    )
+    model.add_reaction([SP, PS, SA_a, SA_p, PA, AR, RA, RS, mu_P, mu_R, mu_prime_A, mu_S])
+
+    # Timespan
+    tspan = TimeSpan.linspace(t=10, num_points=11)
+    model.timespan(tspan)
+    return model
+
+def create_telegraph_model():
+    model = Model(name="on_off")
+
+    # Species
+    ON = Species(name="ON", initial_value=0)
+    OFF = Species(name="OFF", initial_value=1)
+    model.add_species([ON, OFF])
+
+    # Parameters
+    kon  = Parameter(name="kon",  expression=0.1)
+    koff = Parameter(name="koff", expression=1.0)
+    model.add_parameter([kon, koff])
+
+    # Reactions
+    r1 = Reaction(name="on", reactants={'OFF': 1}, products= {'ON': 1}, rate='kon')
+    r2 = Reaction(name="off", reactants={'ON': 1}, products= {'OFF': 1}, rate='koff')
+    model.add_reaction([r1, r2])
+
+    # Timespan
+    tspan = TimeSpan.linspace(t=100, num_points=100)
+    model.timespan(tspan)
+    return model
diff --git a/stochss/tests/integration_tests/__init__.py b/stochss/tests/integration_tests/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/stochss/tests/integration_tests/test_gillespy2.py b/stochss/tests/integration_tests/test_gillespy2.py
new file mode 100644
index 0000000000..7c2112a18e
--- /dev/null
+++ b/stochss/tests/integration_tests/test_gillespy2.py
@@ -0,0 +1,158 @@
+'''
+StochSS is a platform for simulating biochemical systems
+Copyright (C) 2019-2023 StochSS developers.
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+'''
+import os
+import unittest
+
+import gillespy2
+
+from example_models import (
+    create_vilar_oscillator, create_dimerization, create_trichloroethylene, create_lac_operon, create_schlogl,
+    create_michaelis_menten, create_toggle_switch, create_decay, create_tyson_2_state_oscillator,
+    create_oregonator, create_opioid, create_telegraph_model
+)
+
+os.chdir('/stochss')
+
+# pylint: disable=import-outside-toplevel
+class TestGillesPy2Dependency(unittest.TestCase):
+    '''
+    ####################################################################################################################
+    Unit tests for GillesPy2 dependency.
+    ####################################################################################################################
+    '''
+    def setUp(self):
+        ''' Create a list of common example paths for each test. '''
+        self.test_ode_models = [
+            create_vilar_oscillator, create_dimerization, create_trichloroethylene, create_lac_operon, create_schlogl,
+            create_michaelis_menten, create_toggle_switch, create_decay, create_tyson_2_state_oscillator,
+            create_oregonator, create_opioid, create_telegraph_model
+        ]
+        self.test_ssa_models = [
+            create_vilar_oscillator, create_dimerization, create_trichloroethylene, create_schlogl,
+            create_michaelis_menten, create_toggle_switch, create_decay, create_tyson_2_state_oscillator,
+            create_opioid, create_telegraph_model
+        ]
+        self.test_hybrid_models = [
+            create_dimerization, create_trichloroethylene, create_schlogl, create_michaelis_menten,
+            create_toggle_switch, create_decay, create_opioid, create_telegraph_model
+        ]
+
+    ####################################################################################################################
+    # Unit tests for GillesPy2 dependency check_cpp_support.
+    ####################################################################################################################
+
+    def test_check_cpp_support(self):
+        ''' Check if the check cpp support functions works in StochSS. '''
+        from gillespy2.solvers.utilities.cpp_support_test import check_cpp_support
+
+        self.assertIsInstance(check_cpp_support(), bool)
+
+    ####################################################################################################################
+    # Unit tests for GillesPy2 dependency get_best_solver.
+    ####################################################################################################################
+
+    def test_get_best_solver(self):
+        ''' Check if the get best solver function works in StochSS. '''
+        test_model = self.test_ode_models[0]()
+        test_solver = test_model.get_best_solver()
+        self.assertIsInstance(test_solver(model=test_model), gillespy2.GillesPySolver)
+
+    ####################################################################################################################
+    # Unit tests for GillesPy2 dependency get_best_solver_algo.
+    ####################################################################################################################
+
+    def test_get_best_solver_algo(self):
+        ''' Check if the get best solver algo function works in StochSS. '''
+        test_algos = ["ODE", "SSA", "CLE", "Tau-Leaping", "Tau-Hybrid"]
+        test_model = self.test_ode_models[0]()
+        for test_algo in test_algos:
+            with self.subTest(test_algo=test_algo):
+                test_solver = test_model.get_best_solver_algo(algorithm=test_algo)
+                self.assertIsInstance(test_solver(model=test_model), gillespy2.GillesPySolver)
+
+    ####################################################################################################################
+    # Unit tests for GillesPy2 dependency solvers.
+    ####################################################################################################################
+
+    def test_ode_solver(self):
+        ''' Check if the test_models run with the ODESolver. '''
+        for model in self.test_ode_models:
+            test_model = model()
+            msg = f"Running {test_model.name} using {gillespy2.ODESolver.name} failed!"
+            with self.subTest(msg=msg):
+                test_model.run(solver=gillespy2.ODESolver, timeout=30)
+
+    def test_ode_c_solver(self):
+        ''' Check if the test_models run with the ODECSolver. '''
+        for model in self.test_ode_models:
+            test_model = model()
+            test_solver = gillespy2.ODECSolver(model=test_model)
+            msg = f"Running {test_model.name} using {gillespy2.ODECSolver.name} failed!"
+            with self.subTest(msg=msg):
+                test_model.run(solver=test_solver, timeout=30)
+
+    def test_numpy_ssa_solver(self):
+        ''' Check if the test_models run with the NumPySSASolver. '''
+        for model in self.test_ssa_models:
+            test_model = model()
+            msg = f"Running {test_model.name} using {gillespy2.NumPySSASolver.name} failed!"
+            with self.subTest(msg=msg):
+                test_model.run(solver=gillespy2.NumPySSASolver, timeout=30)
+
+    def test_ssa_c_solver(self):
+        ''' Check if the test_models run with the SSACSolver. '''
+        for model in self.test_ssa_models:
+            test_model = model()
+            test_solver = gillespy2.SSACSolver(model=test_model)
+            msg = f"Running {test_model.name} using {gillespy2.SSACSolver.name} failed!"
+            with self.subTest(msg=msg):
+                test_model.run(solver=test_solver, timeout=30)
+
+    def test_tau_leaping_solver(self):
+        ''' Check if the test_models run with the TauLeapingSolver. '''
+        for model in self.test_ssa_models:
+            test_model = model()
+            msg = f"Running {test_model.name} using {gillespy2.TauLeapingSolver.name} failed!"
+            with self.subTest(msg=msg):
+                test_model.run(solver=gillespy2.TauLeapingSolver, timeout=30)
+
+    def test_tau_leaping_c_solver(self):
+        ''' Check if the test_models run with the TauLeapingCSolver. '''
+        for model in self.test_ssa_models:
+            test_model = model()
+            test_solver = gillespy2.TauLeapingCSolver(model=test_model)
+            msg = f"Running {test_model.name} using {gillespy2.TauLeapingCSolver.name} failed!"
+            with self.subTest(msg=msg):
+                test_model.run(solver=test_solver, timeout=30)
+
+    def test_tau_hybrid_solver(self):
+        ''' Check if the test_models run with the TauHybridSolver. '''
+        for model in self.test_hybrid_models:
+            test_model = model()
+            msg = f"Running {test_model.name} using {gillespy2.TauHybridSolver.name} failed!"
+            with self.subTest(msg=msg):
+                test_model.run(solver=gillespy2.TauHybridSolver, timeout=30)
+
+    def test_tau_hybrid_c_solver(self):
+        ''' Check if the test_models run with the TauHybridCSolver. '''
+        for model in self.test_hybrid_models:
+            test_model = model()
+            test_solver = gillespy2.TauHybridCSolver(model=test_model)
+            msg = f"Running {test_model.name} using {gillespy2.TauHybridCSolver.name} failed!"
+            with self.subTest(msg=msg):
+                test_model.run(solver=test_solver, timeout=30)
diff --git a/stochss/tests/run_all_tests.py b/stochss/tests/run_all_tests.py
new file mode 100755
index 0000000000..a4154ada73
--- /dev/null
+++ b/stochss/tests/run_all_tests.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python3
+
+'''
+StochSS is a platform for simulating biochemical systems
+Copyright (C) 2019-2023 StochSS developers.
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+'''
+import os
+import argparse
+
+parser = argparse.ArgumentParser()
+parser.add_argument(
+    '-m', '--mode', default='staging', choices=['staging', 'release'],
+    help='Run tests in staging mode or release mode.'
+)
+
+if __name__ == "__main__":
+    os.chdir('/stochss')
+    args = parser.parse_args()
+    print(os.path.dirname(__file__))
+
+    from run_unit_tests import run as run_unit
+    from run_integration_tests import run as run_integration
+
+    run_unit(mode=args.mode)
+    run_integration(mode=args.mode)
diff --git a/stochss/tests/run_integration_tests.py b/stochss/tests/run_integration_tests.py
new file mode 100755
index 0000000000..2e4989268b
--- /dev/null
+++ b/stochss/tests/run_integration_tests.py
@@ -0,0 +1,54 @@
+#!/usr/bin/env python3
+
+'''
+StochSS is a platform for simulating biochemical systems
+Copyright (C) 2019-2023 StochSS developers.
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+'''
+import os
+import sys
+import unittest
+import argparse
+
+parser = argparse.ArgumentParser()
+parser.add_argument(
+    '-m', '--mode', default='staging', choices=['staging', 'release'],
+    help='Run integration tests in staging mode or release mode.'
+)
+
+def run(mode=None):
+    ''' Run the integration tests. '''
+    from integration_tests import test_gillespy2
+
+    modules = [
+        test_gillespy2
+    ]
+
+    for module in modules:
+        suite = unittest.TestLoader().loadTestsFromModule(module)
+        runner = unittest.TextTestRunner(failfast=mode == 'staging', verbosity=1)
+
+        print("Executing: {}".format(module))
+        result = runner.run(suite)
+        print('=' * 70)
+        if not result.wasSuccessful():
+            sys.exit(not result.wasSuccessful())
+
+if __name__ == "__main__":
+    os.chdir('/stochss')
+    args = parser.parse_args()
+    print(os.path.dirname(__file__))
+
+    run(mode=args.mode)
diff --git a/stochss/tests/run_tests.py b/stochss/tests/run_unit_tests.py
similarity index 73%
rename from stochss/tests/run_tests.py
rename to stochss/tests/run_unit_tests.py
index b794f7c240..0187a4367a 100755
--- a/stochss/tests/run_tests.py
+++ b/stochss/tests/run_unit_tests.py
@@ -17,43 +17,44 @@
 You should have received a copy of the GNU General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 '''
-
-
-
-import unittest
-import sys
 import os
+import sys
+import unittest
 import argparse
 
 parser = argparse.ArgumentParser()
-parser.add_argument('-m', '--mode', default='staging', choices=['staging', 'release'],
-                    help='Run tests in staging mode or release mode.')
-
-if __name__ == "__main__":
-    os.chdir('/stochss')
-    args = parser.parse_args()
-    print(os.path.dirname(__file__))
-
-    import test_model_template
-    import test_settings_template
-    import test_stochss_base
-    import test_stochss_file
-    import test_gillespy2
-
+parser.add_argument(
+    '-m', '--mode', default='staging', choices=['staging', 'release'],
+    help='Run unit tests in staging mode or release mode.'
+)
+
+def run(mode=None):
+    ''' Run the unit tests. '''
+    from unit_tests import test_model_template
+    from unit_tests import test_settings_template
+    from unit_tests import test_stochss_base
+    from unit_tests import test_stochss_file
+    
     modules = [
         test_model_template,
         test_settings_template,
-        test_gillespy2,
         test_stochss_base,
         test_stochss_file
     ]
 
     for module in modules:
         suite = unittest.TestLoader().loadTestsFromModule(module)
-        runner = unittest.TextTestRunner(failfast=args.mode == 'staging', verbosity=1)
+        runner = unittest.TextTestRunner(failfast=mode == 'staging', verbosity=1)
 
         print("Executing: {}".format(module))
         result = runner.run(suite)
         print('=' * 70)
         if not result.wasSuccessful():
             sys.exit(not result.wasSuccessful())
+
+if __name__ == "__main__":
+    os.chdir('/stochss')
+    args = parser.parse_args()
+    print(os.path.dirname(__file__))
+
+    run(mode=args.mode)
diff --git a/stochss/tests/test_gillespy2.py b/stochss/tests/test_gillespy2.py
deleted file mode 100644
index 9741595518..0000000000
--- a/stochss/tests/test_gillespy2.py
+++ /dev/null
@@ -1,182 +0,0 @@
-'''
-StochSS is a platform for simulating biochemical systems
-Copyright (C) 2019-2023 StochSS developers.
-
-This program is free software: you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation, either version 3 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with this program.  If not, see <http://www.gnu.org/licenses/>.
-'''
-
-import os
-import unittest
-
-from gillespy2 import GillesPySolver
-
-from example_models import (
-    Brusselator,
-    Degradation,
-    Dimerization,
-    LotkavolterraOscillator,
-    MichaelisMenten,
-    Opioid,
-    Schlogl,
-    ToggleSwitch,
-    VilarOscillator,
-    Oregonator,
-    TysonOscillator
-)
-
-os.chdir('/stochss')
-
-# pylint: disable=import-outside-toplevel
-class TestGillesPy2Dependency(unittest.TestCase):
-    '''
-    ################################################################################################
-    Unit tests for GillesPy2 dependency.
-    ################################################################################################
-    '''
-    def setUp(self):
-        ''' Create a list of common example paths for each test. '''
-        self.test_models = [Brusselator, Degradation, Dimerization, LotkavolterraOscillator,
-                            MichaelisMenten, Opioid, Schlogl, ToggleSwitch, VilarOscillator]
-
-    ################################################################################################
-    # Unit tests for GillesPy2 dependency check_cpp_support.
-    ################################################################################################
-
-    def test_check_cpp_support(self):
-        ''' Check if the check cpp support functions works in StochSS. '''
-        from gillespy2.solvers.utilities.cpp_support_test import check_cpp_support
-
-        self.assertIsInstance(check_cpp_support(), bool)
-
-    ################################################################################################
-    # Unit tests for GillesPy2 dependency get_best_solver.
-    ################################################################################################
-
-    def test_get_best_solver(self):
-        ''' Check if the get best solver function works in StochSS. '''
-        test_model = self.test_models[0]()
-        test_solver = test_model.get_best_solver()
-        self.assertIsInstance(test_solver(model=test_model), GillesPySolver)
-
-    ################################################################################################
-    # Unit tests for GillesPy2 dependency get_best_solver_algo.
-    ################################################################################################
-
-    def test_get_best_solver_algo(self):
-        ''' Check if the get best solver algo function works in StochSS. '''
-        test_algos = ["ODE", "SSA", "CLE", "Tau-Leaping", "Tau-Hybrid"]
-        test_model = self.test_models[0]()
-        for test_algo in test_algos:
-            with self.subTest(test_algo=test_algo):
-                test_solver = test_model.get_best_solver_algo(algorithm=test_algo)
-                self.assertIsInstance(test_solver(model=test_model), GillesPySolver)
-
-    ################################################################################################
-    # Unit tests for GillesPy2 dependency solvers.
-    ################################################################################################
-
-    def test_ode_solver(self):
-        ''' Check if the test_models run with the ODESolver. '''
-        from gillespy2 import ODESolver
-
-        self.test_models.append(Oregonator)
-        for model in self.test_models:
-            test_model = model()
-            with self.subTest(model=test_model.name):
-                test_model.run(solver=ODESolver)
-
-
-    def test_ode_c_solver(self):
-        ''' Check if the test_models run with the ODECSolver. '''
-        from gillespy2 import ODECSolver
-
-        self.test_models.append(Oregonator)
-        for model in self.test_models:
-            test_model = model()
-            with self.subTest(model=test_model.name):
-                test_solver = ODECSolver(model=test_model)
-                test_model.run(solver=test_solver)
-
-
-    def test_numpy_ssa_solver(self):
-        ''' Check if the test_models run with the NumPySSASolver. '''
-        from gillespy2 import NumPySSASolver
-
-        self.test_models.append(TysonOscillator)
-        for model in self.test_models:
-            test_model = model()
-            with self.subTest(model=test_model.name):
-                test_model.run(solver=NumPySSASolver)
-
-
-    def test_ssa_c_solver(self):
-        ''' Check if the test_models run with the SSACSolver. '''
-        from gillespy2 import SSACSolver
-
-        self.test_models.append(TysonOscillator)
-        for model in self.test_models:
-            test_model = model()
-            with self.subTest(model=test_model.name):
-                test_solver = SSACSolver(model=test_model)
-                test_model.run(solver=test_solver)
-
-
-    def test_tau_leaping_solver(self):
-        ''' Check if the test_models run with the TauLeapingSolver. '''
-        from gillespy2 import TauLeapingSolver
-
-        self.test_models.append(TysonOscillator)
-        for model in self.test_models:
-            test_model = model()
-            with self.subTest(model=test_model.name):
-                test_model.run(solver=TauLeapingSolver)
-
-
-    def test_tau_leaping_c_solver(self):
-        ''' Check if the test_models run with the TauLeapingCSolver. '''
-        from gillespy2 import TauLeapingCSolver
-
-        self.test_models.append(TysonOscillator)
-        for model in self.test_models:
-            test_model = model()
-            with self.subTest(model=test_model.name):
-                test_solver = TauLeapingCSolver(model=test_model)
-                test_model.run(solver=test_solver)
-
-
-    def test_tau_hybrid_solver(self):
-        ''' Check if the test_models run with the TauHybridSolver. '''
-        from gillespy2 import TauHybridSolver
-
-        self.test_models.append(Oregonator)
-        self.test_models.append(TysonOscillator)
-        self.test_models.remove(VilarOscillator)
-        for model in self.test_models:
-            test_model = model()
-            with self.subTest(model=test_model.name):
-                test_model.run(solver=TauHybridSolver)
-
-
-    def test_tau_hybrid_c_solver(self):
-        ''' Check if the test_models run with the TauHybridCSolver. '''
-        from gillespy2 import TauHybridCSolver
-
-        self.test_models.append(Oregonator)
-        self.test_models.append(TysonOscillator)
-        self.test_models.remove(VilarOscillator)
-        self.test_models.remove(ToggleSwitch)
-        for model in self.test_models:
-            test_model = model()
-            with self.subTest(model=test_model.name):
-                test_model.run(solver=TauHybridCSolver)
diff --git a/stochss/tests/unit_tests/__init__.py b/stochss/tests/unit_tests/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/stochss/tests/test_model_template.py b/stochss/tests/unit_tests/test_model_template.py
similarity index 100%
rename from stochss/tests/test_model_template.py
rename to stochss/tests/unit_tests/test_model_template.py
diff --git a/stochss/tests/test_settings_template.py b/stochss/tests/unit_tests/test_settings_template.py
similarity index 100%
rename from stochss/tests/test_settings_template.py
rename to stochss/tests/unit_tests/test_settings_template.py
diff --git a/stochss/tests/test_stochss_base.py b/stochss/tests/unit_tests/test_stochss_base.py
similarity index 100%
rename from stochss/tests/test_stochss_base.py
rename to stochss/tests/unit_tests/test_stochss_base.py
diff --git a/stochss/tests/test_stochss_file.py b/stochss/tests/unit_tests/test_stochss_file.py
similarity index 100%
rename from stochss/tests/test_stochss_file.py
rename to stochss/tests/unit_tests/test_stochss_file.py