diff --git a/CITATION b/CITATION
index 89bd8489e45e1576db9690fc803b09cb05ed0c15..035ed60dc2d6c7d3f1b9c0f25658bf385fdb83cc 100644
--- a/CITATION
+++ b/CITATION
@@ -6,25 +6,29 @@
     year = {2017},
     month = {02},
     volume = {13},
-    pages = {1-25},
+    pages = {e1005179},
     number = {2},
     doi = {10.1371/jour.pcbi.1005179}
 }
 
 @article{Schmidt18_1,
-    title={Multi-scale account of the network structure of macaque visual cortex},
-    author={Schmidt, Maximilian and Bakker, Rembrandt and Hilgetag, Claus Christian and Diesmann, Markus and van Albada, Sacha Jennifer},
-    journal={Brain Structure and Function},
-    year={2018},
-    volume={223},
-    number={3},
-    pages={1409--1435},
-    doi={https://doi.org/10.1007/s00429-017-1554-4}
+    title = {Multi-scale account of the network structure of macaque visual cortex},
+    author = {Schmidt, Maximilian and Bakker, Rembrandt and Hilgetag, Claus Christian and Diesmann, Markus and van Albada, Sacha Jennifer},
+    journal = {Brain Structure and Function},
+    year = {2018},
+    volume = {223},
+    number = {3},
+    pages = {1409--1435},
+    doi = {https://doi.org/10.1007/s00429-017-1554-4}
 }
 
 @article{Schmidt18_2,
-  title={A multi-scale layer-resolved spiking network model of resting-state dynamics in macaque visual cortical areas}
-  author={Schmidt, Maximilian and Bakker, Rembrandt and Shen, Kelly and Bezgin, Gleb and Diesmann, Markus and Sacha Jennifer van Albada},
-  year={2018},
-  journal={under review}
- }
+  title = {A multi-scale layer-resolved spiking network model of resting-state dynamics in macaque visual cortical areas}
+  author = {Schmidt, Maximilian and Bakker, Rembrandt and Shen, Kelly and Bezgin, Gleb and Diesmann, Markus and Sacha Jennifer van Albada},
+  year = {in press},
+  journal = {PLOS Computational Biology},
+  volume = {14},
+  number = {10},
+  pages = {e1006359},
+  doi = {10.1371/journal.pcbi.1006359}
+}
diff --git a/README.md b/README.md
index 552ed0221736d22dafa59f21534793ac6cb517a1..0bc6574d3db58fa024f10fb62be3c75f8c13c156 100644
--- a/README.md
+++ b/README.md
@@ -4,20 +4,20 @@
 ![Model overview](model_construction.png)
 
 This code implements the spiking network model of macaque visual cortex developed
-at the Institute of Neuroscience and Medicine (INM-6), Research Center Jülich. 
+at the Institute of Neuroscience and Medicine (INM-6), Research Center Jülich.
 The model has been documented in the following publications:
 
 1. Schmidt M, Bakker R, Hilgetag CC, Diesmann M & van Albada SJ
    Multi-scale account of the network structure of macaque visual cortex
-   Brain Structure and Function (2018), 223:1409 [https://doi.org/10.1007/s00429-017-1554-4](https://doi.org/10.1007/s00429-017-1554-4)
+   Brain Structure and Function (2018), 223: 1409 [https://doi.org/10.1007/s00429-017-1554-4](https://doi.org/10.1007/s00429-017-1554-4)
 
 2. Schuecker J, Schmidt M, van Albada SJ, Diesmann M & Helias M (2017)
    Fundamental Activity Constraints Lead to Specific Interpretations of the Connectome.
-   PLOS Computational Biology, 13(2). [https://doi.org/10.1371/journal.pcbi.1005179](https://doi.org/10.1371/journal.pcbi.1005179)
+   PLOS Computational Biology, 13(2): e1005179. [https://doi.org/10.1371/journal.pcbi.1005179](https://doi.org/10.1371/journal.pcbi.1005179)
 
-3. Schmidt M, Bakker R, Shen K, Bezgin B, Diesmann M & van Albada SJ (2018) 
+3. Schmidt M, Bakker R, Shen K, Bezgin B, Diesmann M & van Albada SJ (accepted)
    A multi-scale layer-resolved spiking network model of
-   resting-state dynamics in macaque cortex. PLOS Computational Biology (accepted)
+   resting-state dynamics in macaque cortex. PLOS Computational Biology, 14(9): e1006359.
 
 The code in this repository is self-contained and allows one to
 reproduce the results of all three papers.
@@ -49,7 +49,7 @@ Furthermore, please add the path to the repository to your PYTHONPATH:
 `MultiAreaModel`
 
 The central class that initializes the network and contains all
-information about population sizes and network connectivity. This 
+information about population sizes and network connectivity. This
 enables reproducing all figures in [1]. Network parameters only
 refer to the structure of the network and ignore any information on
 its dynamical simulation or description via analytical theory.
@@ -82,9 +82,9 @@ basic analysis and plotting.
 
 The `figures` folder contains a subfolder with all scripts necessary to produce
 the figures from [1]. The scripts for [2] and [3] will follow soon.
-If snakemake is installed, the figures can be produced by executing
+If Snakemake (Köster J & Rahmann S, Bioinformatics (2012) 28(19): 2520-2522) is installed, the figures can be produced by executing
 `snakemake` in the respective folder, e.g.:
-	
+
 	cd figures/Schmidt2018/
 	snakemake
 
@@ -97,9 +97,9 @@ A simple simulation can be run in the following way:
        custom_params = ...
        custom_simulation_params = ...
 2. Instantiate the model class together with a simulation class instance.
-   
+
        M = MultiAreaModel(custom_params, simulation=True, sim_spec=custom_simulation_params)
-	   
+
 3. Start the simulation.
 
        M.simulation.simulate()
@@ -113,13 +113,13 @@ The procedure is similar to a simple simulation:
        custom_params = ...
        custom_simulation_params = ...
 2. Instantiate the model class together with a simulation class instance.
-   
+
        M = MultiAreaModel(custom_params, simulation=True, sim_spec=custom_simulation_params)
 3. Start the simulation.
    Call `start_job` to create a job file using the `jobscript_template` from the configuration file
    and submit it to the queue with the user-defined `submit_cmd`.
-   
-The file `run_example.py` provides an example.
+
+The file `run_example_fullscale.py` provides an example.
 
 Be aware that, depending on the chosen parameters and initial conditions, the network can enter a high-activity state, which slows down the simulation drastically and can cost a significant amount of computing resources.
 
@@ -150,7 +150,7 @@ The multi-area model can be run in different modes.
    - `hom_poisson_stat`: all non-simulated areas are replaced by Poissonian spike trains with the
      same rate as the stationary background input (`rate_ext` in `input_params`).
    - `het_poisson_stat`: all non-simulated areas are replaced by Poissonian spike trains with
-      population-specific stationary rate stored in an external file. 
+      population-specific stationary rate stored in an external file.
    - `current_nonstat`: all non-simulated areas are replaced by stepwise constant currents with
      population-specific, time-varying time series defined in an external file.
 
@@ -210,7 +210,7 @@ Meegen.
 
 ## Citation
 
-If you use this code, we ask you to cite the appropriate papers in your publication. For the multi-area model itself, please cite [1] and [3]. If you use the mean-field theory or the stabilization method, please cite [2] in addition. We provide bibtex entries in `CITATION`.
+If you use this code, we ask you to cite the appropriate papers in your publication. For the multi-area model itself, please cite [1] and [3]. If you use the mean-field theory or the stabilization method, please cite [2] in addition. We provide bibtex entries in the file called `CITATION`.
 
 If you have questions regarding the code or scientific content, please create an issue on github.
 
@@ -225,7 +225,7 @@ This work was supported by the Helmholtz Portfolio Supercomputing and
 Modeling for the Human Brain (SMHB), the European Union 7th Framework
 Program (Grant 269921, BrainScaleS and 604102, Human Brain Project,
 Ramp up phase) and European Unions Horizon 2020 research and
-innovation program (Grant 720270, Human Brain Project, SGA1), the
+innovation program (Grants 720270 and 737691, Human Brain Project, SGA1 and SGA2), the
 Jülich Aachen Research Alliance (JARA), the Helmholtz young
 investigator group VH-NG-1028,and the German Research Council (DFG
 Grants SFB936/A1,Z1 and TRR169/A2) and computing time granted by the
diff --git a/_config.yml b/_config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..c7da01a893f5f9aa118959607a3f2c2b9cff047f
--- /dev/null
+++ b/_config.yml
@@ -0,0 +1,3 @@
+theme: jekyll-theme-dinky
+description: A large-scale spiking model of the vision-related areas of macaque cortex. 
+future: true
diff --git a/run_example.py b/run_example_downscaled.py
similarity index 51%
rename from run_example.py
rename to run_example_downscaled.py
index 9a16c376ae9b3c8ab8f38989b5c7d75135b34889..d5f278175460e345009752bf4e995ea2fc1e19e8 100644
--- a/run_example.py
+++ b/run_example_downscaled.py
@@ -2,56 +2,8 @@ import numpy as np
 import os
 
 from multiarea_model import MultiAreaModel
-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 = MultiAreaModel(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)
-
-
 """
 Down-scaled model.
 Neurons and indegrees are both scaled down to 10 %.
diff --git a/run_example_fullscale.py b/run_example_fullscale.py
new file mode 100644
index 0000000000000000000000000000000000000000..30e0c8bcbf9877a52a65a8ca699c889c5011b0bc
--- /dev/null
+++ b/run_example_fullscale.py
@@ -0,0 +1,52 @@
+import numpy as np
+import os
+
+from multiarea_model import MultiAreaModel
+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 = MultiAreaModel(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)