diff --git a/multiarea_model/__init__.py b/multiarea_model/__init__.py
index a5d3ddec2c868de2bc7b0a3ed282930f650981a2..ede7d376c9264ea248928b906ca07314ada7202c 100644
--- a/multiarea_model/__init__.py
+++ b/multiarea_model/__init__.py
@@ -1,4 +1,6 @@
 from . import default_params
+from .default_params import get_default_parameters
+from .default_params import nested_update
 from .analysis import Analysis
 from .simulation import Simulation
 from .theory import Theory
diff --git a/multiarea_model/data_multiarea/Model.py b/multiarea_model/data_multiarea/Model.py
index 6d68fba7906b192ba3631f498599ababfcf5f8b3..661d483ff17cebe0b3a86f19bf533464ef0f8c6b 100644
--- a/multiarea_model/data_multiarea/Model.py
+++ b/multiarea_model/data_multiarea/Model.py
@@ -20,19 +20,16 @@ Sacha van Albada
 import numpy as np
 import json
 import re
-import sys
 import os
 import scipy
 import scipy.integrate
-import pprint
 from copy import deepcopy
 from nested_dict import nested_dict
 from itertools import product
-from multiarea_model.default_params import network_params, nested_update
 from multiarea_model.data_multiarea.VisualCortex_Data import process_raw_data
 
 
-def compute_Model_params(out_label='', mode='default'):
+def compute_Model_params(net_params, out_label='', mode='default'):
     """
     Compute the parameters of the network, in particular the size
     of populations, external inputs to them, and number of synapses
@@ -119,24 +116,10 @@ def compute_Model_params(out_label='', mode='default'):
     overwrite default by custom values for parameters specified
     in the parameter file.
     """
-    net_params = deepcopy(network_params)
     if mode == 'default':
         prefix = 'default'
     elif mode == 'custom':
         prefix = 'custom_data_files/custom'
-        with open(os.path.join(basepath, '.'.join(('_'.join((prefix,
-                                                             out_label,
-                                                             'parameter_dict')),
-                                                   'json'))), 'r') as f:
-            custom_params = json.load(f)
-        nested_update(net_params, custom_params)
-        # print information on overwritten parameters
-        print("\n")
-        print("========================================")
-        print("Customized parameters")
-        print("--------------------")
-        pprint.pprint(custom_params)
-        print("========================================")
 
     """
     Define parameter values
diff --git a/multiarea_model/default_params.py b/multiarea_model/default_params.py
index bbe72df88b54484b920e5030468c77cd076ea747..cef39e5c6a8c17b33b0a53deb29b6f227ebd8162 100644
--- a/multiarea_model/default_params.py
+++ b/multiarea_model/default_params.py
@@ -310,3 +310,10 @@ def check_custom_params(d, def_d):
                 def_val = def_d[key]
             except KeyError:
                 raise KeyError('Unused key in custom parameter dictionary: {}'.format(key))
+
+def get_default_parameters():
+    parameters = {
+            'sim_dict': sim_params,
+            'network_dict': network_params
+            }
+    return parameters
diff --git a/multiarea_model/multiarea_model.py b/multiarea_model/multiarea_model.py
index e6766d8060724b58ef63e1e16ecffec4917c14ef..27c5b33b6896e452e35ab49b290c46dd42de6bf1 100644
--- a/multiarea_model/multiarea_model.py
+++ b/multiarea_model/multiarea_model.py
@@ -30,10 +30,8 @@ import numpy as np
 import os
 import pprint
 import shutil
-from .default_params import complete_area_list, nested_update, network_params
-from .default_params import check_custom_params
+from .default_params import complete_area_list
 from collections import OrderedDict
-from copy import deepcopy
 from .data_multiarea.Model import compute_Model_params
 from .analysis import Analysis
 from config import base_path
@@ -61,7 +59,7 @@ dicthash.FLOOR_SMALL_FLOATS = True
 
 
 class Model:
-    def __init__(self, network_spec, theory=False, simulation=False,
+    def __init__(self, network_spec, theory=False, simulation=True,
                  analysis=False, *args, **keywords):
         """
         Multiarea model class.
@@ -83,27 +81,20 @@ class Model:
             whether to create an instance of the analysis class as member.
 
         """
-        self.params = deepcopy(network_params)
+        self.params = network_spec['network_dict']
+        self.sim_spec = network_spec['sim_dict']
         if isinstance(network_spec, dict):
-            print("Initializing network from dictionary.")
-            check_custom_params(network_spec, self.params)
-            self.custom_params = network_spec
             p_ = 'multiarea_model/data_multiarea/custom_data_files'
             # Draw random integer label for data script to avoid clashes with
             # parallelly created class instances
             rand_data_label = np.random.randint(10000)
             print("RAND_DATA_LABEL", rand_data_label)
-            tmp_parameter_fn = os.path.join(base_path,
-                                            p_,
-                                            'custom_{}_parameter_dict.json'.format(rand_data_label))
             tmp_data_fn = os.path.join(base_path,
                                        p_,
                                        'custom_Data_Model_{}.json'.format(rand_data_label))
 
-            with open(tmp_parameter_fn, 'w') as f:
-                json.dump(self.custom_params, f)
             # Execute Data script
-            compute_Model_params(out_label=str(rand_data_label),
+            compute_Model_params(self.params, out_label=str(rand_data_label),
                                  mode='custom')
         else:
             print("Initializing network from label.")
@@ -115,7 +106,6 @@ class Model:
                                        'custom_Data_Model_{}.json'.format(network_spec))
             with open(parameter_fn, 'r') as f:
                 self.custom_params = json.load(f)
-        nested_update(self.params, self.custom_params)
         with open(tmp_data_fn, 'r') as f:
             dat = json.load(f)
 
@@ -175,8 +165,6 @@ class Model:
                                    'config_files',
                                    'custom_Data_Model_{}.json'.format(self.label))
 
-            shutil.move(tmp_parameter_fn,
-                        parameter_fn)
             shutil.move(tmp_data_fn,
                         data_fn)
 
@@ -191,12 +179,7 @@ class Model:
                 theory_spec = keywords['theory_spec']
             self.init_theory(theory_spec)
 
-        if simulation:
-            if 'sim_spec' not in keywords:
-                sim_spec = {}
-            else:
-                sim_spec = keywords['sim_spec']
-            self.init_simulation(sim_spec)
+        self.init_simulation()
 
         if analysis:
             assert(getattr(self, 'simulation'))
@@ -212,8 +195,8 @@ class Model:
     def connect(self):
         self.simulation.connect()
 
-    def simulate(self):
-        self.simulation.simulate()
+    def simulate(self, t_sim):
+        self.simulation.simulate(t_sim)
 
     def __str__(self):
         s = "Multi-area network {} with custom parameters: \n".format(self.label)
@@ -229,8 +212,8 @@ class Model:
     def init_theory(self, theory_spec):
         self.theory = Theory(self, theory_spec)
 
-    def init_simulation(self, sim_spec):
-        self.simulation = Simulation(self, sim_spec)
+    def init_simulation(self):
+        self.simulation = Simulation(self, self.sim_spec)
 
     def init_analysis(self, ana_spec):
         assert(hasattr(self, 'simulation'))
diff --git a/multiarea_model/simulation.py b/multiarea_model/simulation.py
index 66f719c9b374d5d55123a5eb302eac5c4ee9c653..a0f9d25071dea0c158cd03b9d0dceef546205ccc 100644
--- a/multiarea_model/simulation.py
+++ b/multiarea_model/simulation.py
@@ -82,7 +82,7 @@ class Simulation:
         self.copy_files()
         print("Copied files.")
         d = {'sim_params': self.custom_params,
-             'network_params': self.network.custom_params,
+             'network_params': self.network.params,
              'network_label': self.network.label}
         with open(os.path.join(self.data_dir,
                                '_'.join(('custom_params', self.label))), 'w') as f:
@@ -91,7 +91,6 @@ class Simulation:
 
         self.areas_simulated = self.params['areas_simulated']
         self.areas_recorded = self.params['recording_dict']['areas_recorded']
-        self.T = self.params['t_sim']
 
         self.prepare()
 
@@ -124,8 +123,6 @@ class Simulation:
                               'simulation.py'),
                  os.path.join('multiarea_model',
                               'default_params.py'),
-                 os.path.join('config_files',
-                              ''.join(('custom_Data_Model_', self.network.label, '.json'))),
                  os.path.join('config_files',
                               '_'.join((self.network.label, 'config')))]
         if self.network.params['connection_params']['replace_cc_input_source'] is not None:
@@ -296,11 +293,11 @@ class Simulation:
         self.cortico_cortical_input()
         self.save_network_gids()
 
-    def simulate(self):
+    def simulate(self, t_sim):
         """
         Simulate the network.
         """
-        nest.Simulate(self.T)
+        nest.Simulate(t_sim)
 
     def memory(self):
         """
diff --git a/run_example_downscaled.py b/run_example_downscaled.py
deleted file mode 100644
index c875439494d54634e88a82710595143f224cfacc..0000000000000000000000000000000000000000
--- a/run_example_downscaled.py
+++ /dev/null
@@ -1,46 +0,0 @@
-import numpy as np
-import os
-
-import multiarea_model as model
-from config import base_path
-
-"""
-Down-scaled model.
-Neurons and indegrees are both scaled down to 10 %.
-Can usually be simulated on a local machine.
-
-Warning: This will not yield reasonable dynamical results from the
-network and is only meant to demonstrate the simulation workflow.
-"""
-d = {}
-conn_params = {'replace_non_simulated_areas': 'het_poisson_stat',
-               'g': -11.,
-               'K_stable': 'K_stable.npy',
-               'fac_nu_ext_TH': 1.2,
-               'fac_nu_ext_5E': 1.125,
-               'fac_nu_ext_6E': 1.41666667,
-               'av_indegree_V1': 3950.}
-input_params = {'rate_ext': 10.}
-neuron_params = {'V0_mean': -150.,
-                 'V0_sd': 50.}
-network_params = {'N_scaling': 0.01,
-                  'K_scaling': 0.01,
-                  'fullscale_rates': os.path.join(base_path, 'tests/fullscale_rates.json'),
-                  'input_params': input_params,
-                  'connection_params': conn_params,
-                  'neuron_params': neuron_params}
-
-sim_params = {'t_sim': 2000.,
-              'num_processes': 1,
-              'local_num_threads': 1,
-              'recording_dict': {'record_vm': False}}
-
-theory_params = {'dt': 0.1}
-
-M = model.Model(network_params, simulation=True,
-                   sim_spec=sim_params,
-                   theory=True,
-                   theory_spec=theory_params)
-M.create()
-M.connect()
-M.simulate()
diff --git a/run_example_fullscale.py b/run_example_fullscale.py
deleted file mode 100644
index 2b70f9d03bcf6e29f8126465fbe93a5807a5b7bb..0000000000000000000000000000000000000000
--- a/run_example_fullscale.py
+++ /dev/null
@@ -1,52 +0,0 @@
-import numpy as np
-import os
-
-from multiarea_model import Model
-from start_jobs import start_job
-from config import submit_cmd, jobscript_template
-from config import base_path
-
-"""
-Example script showing how to simulate the multi-area model
-on a cluster.
-
-We choose the same configuration as in
-Fig. 3 of Schmidt et al. (2018).
-
-"""
-
-"""
-Full model. Needs to be simulated with sufficient
-resources, for instance on a compute cluster.
-"""
-d = {}
-conn_params = {'g': -11.,
-               'K_stable': os.path.join(base_path, 'K_stable.npy'),
-               'fac_nu_ext_TH': 1.2,
-               'fac_nu_ext_5E': 1.125,
-               'fac_nu_ext_6E': 1.41666667,
-               'av_indegree_V1': 3950.}
-input_params = {'rate_ext': 10.}
-neuron_params = {'V0_mean': -150.,
-                 'V0_sd': 50.}
-network_params = {'N_scaling': 1.,
-                  'K_scaling': 1.,
-                  'connection_params': conn_params,
-                  'input_params': input_params,
-                  'neuron_params': neuron_params}
-
-sim_params = {'t_sim': 2000.,
-              'num_processes': 720,
-              'local_num_threads': 1,
-              'recording_dict': {'record_vm': False}}
-
-theory_params = {'dt': 0.1}
-
-M = Model(network_params, simulation=True,
-                   sim_spec=sim_params,
-                   theory=True,
-                   theory_spec=theory_params)
-p, r = M.theory.integrate_siegert()
-print("Mean-field theory predicts an average "
-      "rate of {0:.3f} spikes/s across all populations.".format(np.mean(r[:, -1])))
-start_job(M.simulation.label, submit_cmd, jobscript_template)
diff --git a/run_groundstate.py b/run_groundstate.py
new file mode 100644
index 0000000000000000000000000000000000000000..97037c219feacebe5cbe1b03e00cb758aea2f7cc
--- /dev/null
+++ b/run_groundstate.py
@@ -0,0 +1,76 @@
+import os
+
+import multiarea_model as model
+from config import base_path
+import time
+import nest
+
+
+time_start = time.time()
+
+conn_params = {'replace_non_simulated_areas': 'het_poisson_stat',
+               'g': -11.,
+               'K_stable': 'K_stable.npy',
+               'fac_nu_ext_TH': 1.2,
+               'fac_nu_ext_5E': 1.125,
+               'fac_nu_ext_6E': 1.41666667,
+               'av_indegree_V1': 3950.}
+input_params = {'rate_ext': 10.}
+neuron_params = {'V0_mean': -150.,
+                 'V0_sd': 50.}
+network_params = {'N_scaling': 0.01,
+                  'K_scaling': 0.01,
+                  'fullscale_rates': os.path.join(
+                      base_path,
+                      'tests/fullscale_rates.json'
+                      ),
+                  'input_params': input_params,
+                  'connection_params': conn_params,
+                  'neuron_params': neuron_params}
+
+sim_params = {'t_sim': 2000.,
+              'num_processes': 1,
+              'local_num_threads': 1,
+              'recording_dict': {'record_vm': False}}
+
+parameters = model.get_default_parameters()
+new_params = {
+        'sim_dict': sim_params,
+        'network_dict': network_params
+        }
+model.nested_update(parameters, new_params)
+
+M = model.Model(parameters)
+
+time_model = time.time()
+
+M.create()
+
+time_create = time.time()
+
+M.connect()
+
+time_connect = time.time()
+
+M.simulate(parameters['sim_dict']['t_sim'])
+
+time_simulate = time.time()
+
+print(
+    '\nTimes of Rank {}:\n'.format(
+        nest.Rank()) +
+    '  Total time:          {:.3f} s\n'.format(
+        time_simulate -
+        time_start) +
+    '  Time to initialize:  {:.3f} s\n'.format(
+        time_model -
+        time_start) +
+    '  Time to create:      {:.3f} s\n'.format(
+        time_create -
+        time_model) +
+    '  Time to connect:     {:.3f} s\n'.format(
+        time_connect -
+        time_create) +
+    '  Time to simulate:    {:.3f} s\n'.format(
+        time_simulate -
+        time_create))
diff --git a/run_metastable.py b/run_metastable.py
new file mode 100644
index 0000000000000000000000000000000000000000..0f4e13579c95c57167fedda02c8778b9a342c16e
--- /dev/null
+++ b/run_metastable.py
@@ -0,0 +1,68 @@
+import multiarea_model as model
+import time
+import nest
+
+
+time_start = time.time()
+
+conn_params = {'g': -11.0,
+               'fac_nu_ext_TH': 1.2,
+               'fac_nu_ext_5E': 1.125,
+               'fac_nu_ext_6E': 1.41666667,
+               'av_indegree_V1': 3950.0,
+               'K_stable': 'K_stable.npy',
+               'cc_weights_factor': 1.9,
+               'cc_weights_I_factor': 2.0}
+neuron_params = {'V0_mean': -150.,
+                 'V0_sd': 50.}
+network_params = {'N_scaling': .01,
+                  'K_scaling': .01,
+                  'connection_params': conn_params,
+                  'neuron_params': neuron_params}
+
+sim_params = {'t_sim': 2000.,
+              'num_processes': 1,
+              'local_num_threads': 1,
+              'recording_dict': {'record_vm': False}}
+
+parameters = model.get_default_parameters()
+new_params = {
+        'sim_dict': sim_params,
+        'network_dict': network_params
+        }
+model.nested_update(parameters, new_params)
+
+M = model.Model(parameters)
+
+time_model = time.time()
+
+M.create()
+
+time_create = time.time()
+
+M.connect()
+
+time_connect = time.time()
+
+M.simulate(parameters['sim_dict']['t_sim'])
+
+time_simulate = time.time()
+
+print(
+    '\nTimes of Rank {}:\n'.format(
+        nest.Rank()) +
+    '  Total time:          {:.3f} s\n'.format(
+        time_simulate -
+        time_start) +
+    '  Time to initialize:  {:.3f} s\n'.format(
+        time_model -
+        time_start) +
+    '  Time to create:      {:.3f} s\n'.format(
+        time_create -
+        time_model) +
+    '  Time to connect:     {:.3f} s\n'.format(
+        time_connect -
+        time_create) +
+    '  Time to simulate:    {:.3f} s\n'.format(
+        time_simulate -
+        time_create))