Merge branch 'master' of https://github.com/mschmidt87/multi-area-model

README.md

@@ -179,9 +179,9 @@ pandas, numpy, nested_dict, matplotlib (2.1.2), scipy, NEST 2.14.0
 
 Optional: seaborn, Sumatra
 
-To install the required packages in a conda environment, execute:
+To install the required packages with pip, execute:
 
-`conda env create -f environment.yaml`
+`pip install -r requirements.txt`
 
 Note that NEST needs to be installed separately, see <http://www.nest-simulator.org/installation/>.
 
@@ -191,8 +191,14 @@ In addition, reproducing the figures of [1] requires networkx, python-igraph, py
 
 In addition, Figure 7 of [1] requires installing the `infomap` package to perform the map equation clustering. See <http://www.mapequation.org/code.html> for all necessary information.
 
+Similarly, reproducing the figures of [3] requires statsmodels, networkx, pyx, python-louvain, which can be installed by executing:
+
+`pip install -r figures/Schmidt2018_dyn/additional_requirements.txt`
+
 The SLN fit in `multiarea_model/data_multiarea/VisualCortex_Data.py` and `figures/Schmidt2018/Fig5_cc_laminar_pattern.py` requires an installation of R and the R library `aod` (<http://cran.r-project.org/package=aod>). Without R installation, both scripts will directly use the resulting values of the fit (see Fig. 5 of [1]).
 
+The calculation of BOLD signals from the simulated firing rates for Fig. 8 of [3] requires an installation of R and the R library `neuRosim` (<https://cran.r-project.org/web/packages/neuRosim/index.html>).
+
 ## Contributors
 
 All authors of the publications [1-3] made contributions to the

environment.yaml

-name: multiarea_model
-dependencies:
- - future
- - numpy
- - matplotlib
- - pandas
- - python==3.6
- - scipy
- - pip:
-   - nested_dict
-   - dicthash

figures/Schmidt2018/README.md

@@ -4,12 +4,14 @@ This folder contains the scripts to reproduce all figures of Schmidt M, Bakker R
 
 ![Model overview](../../model_construction.png)
 
-Please note: Figures 2, 5, and 8 show slight deviations from the published figures in the paper. Published Figures 2 and 5 miss a few data points. This error slipped in during the review process. Importantly, the presented fits are identical in the (correct) figures in this repository and in the manuscript. These deviations thus do not affect the scientific conclusions.
-
 Please note that the placement of areas in Figure 7 will deviate from the published figure, because their location depends on the force-directed algorithm implemented in `igraph` and `python-igraph` does not allow manual setting of the random seed for the algorithm. This is a mere visual issue and does not affect the scientific content.
 
 Please note that, since we currently cannot publish the data on Neuronal Densities, Figures 2 and 5 can currently not be produced and executing it throws an error.
 
+Reproducing the figures requires some additional Python packages listed in `additional_requirements.txt`. They can be installed using pip by executing:
+
+`pip install -r additional_requirements.txt`
+
 If snakemake is installed, the figures can be produced by executing
 
 `snakemake`

figures/Schmidt2018_dyn/Fig2_bistability.eps

@@ -2,7 +2,7 @@
 %%BoundingBox: 13 13 515 374
 %%HiResBoundingBox: 13.1732 13.1732 514.071 373.415
 %%Creator: PyX 0.14.1
-%%CreationDate: Wed May 16 12:18:47 2018
+%%CreationDate: Thu Jun  7 09:15:20 2018
 %%EndComments
 %%BeginProlog
 %%BeginResource: BeginEPSF
@@ -42,7 +42,7 @@ BeginEPSF
 %%HiResBoundingBox: 0.000000 0.000000 510.235200 315.342696
 %%Creator: GPL Ghostscript 922 (ps2write)
 %%LanguageLevel: 2
-%%CreationDate: D:20180516121846+09'00'
+%%CreationDate: D:20180607091517+09'00'
 %%EndComments
 %%BeginProlog
 save
@@ -19320,7 +19320,7 @@ EndEPSF
 BeginEPSF
 113.386 240.945 283.465 131.47 rectclip
 [0.554725 0.000000 0.000000 0.554725 113.385827 241.499607] concat
-%%BeginDocument: phasespace_sketch.eps
+%%BeginDocument: Fig2_bistability_phasespace_sketch.eps
 %!PS-Adobe-3.0 EPSF-3.0
 %%Creator: cairo 1.14.6 (http://cairographics.org)
 %%CreationDate: Thu Jul 27 11:53:48 2017

figures/Schmidt2018_dyn/Fig2_bistability.py

@@ -55,9 +55,9 @@ for key, label in zip(keys, data_labels):
 
 labels = ['A', 'B', 'C']
 
-for i, k in enumerate(['LA', 'HA', 'LA_post']):
-    ax = axes[labels[i]]
-    ax2 = axes[labels[i] + '2']
+for k, key in enumerate(['LA', 'HA', 'LA_post']):
+    ax = axes[labels[k]]
+    ax2 = axes[labels[k] + '2']
     print(k)
     matrix = np.zeros((len(M.area_list), 8))
 
@@ -66,7 +66,7 @@ for i, k in enumerate(['LA', 'HA', 'LA_post']):
             if pop not in M.structure[area]:
                 rate = np.nan
             else:
-                rate = data[k][area][pop][0]
+                rate = data[key][area][pop][0]
 
             if rate == 0.0:
                 rate = 1e-5
@@ -74,11 +74,11 @@ for i, k in enumerate(['LA', 'HA', 'LA_post']):
 
     matrix = np.transpose(matrix)
 
-    if i == 0:
+    if k == 0:
         matrix_plot(panel_factory.figure, ax, matrix, position='left')
         rate_histogram_plot(panel_factory.figure, ax2,
                             matrix, position='left')
-    elif i == 1:
+    elif k == 1:
         matrix_plot(panel_factory.figure, ax, matrix, position='center')
         rate_histogram_plot(panel_factory.figure, ax2,
                             matrix, position='center')
@@ -98,7 +98,7 @@ pyx.text.preamble(r"\usepackage{helvet}")
 c = pyx.canvas.canvas()
 c.insert(pyx.epsfile.epsfile(0.5, 0.5, "Fig2_bistability_mpl.eps", width=17.6))
 c.insert(pyx.epsfile.epsfile(
-    4., 8.5, "phasespace_sketch.eps", width=10.))
+    4., 8.5, "Fig2_bistability_phasespace_sketch.eps", width=10.))
 c.insert(pyx.epsfile.epsfile(1., 3.1, "Epop.eps", width=0.75))
 c.insert(pyx.epsfile.epsfile(1., 2., "Ipop.eps", width=0.75))
 

figures/Schmidt2018_dyn/README.md

@@ -22,6 +22,14 @@ and responses to parametric stimuli in macaque V1, available in the crcns.org da
 
 Figure 8 requires the experimental fMRI data (described by Babapoor-Farrokhran et al. (2013), see Methods section of Schmidt et al. 2018 for more details) that are contained in this repository in `Fig8_exp_func_conn.csv`.
 
+### Requirements
+
+Reproducing the figures requires some additional Python packages listed in `additional_requirements.txt`. They can be installed using pip by executing:
+
+`pip install -r additional_requirements.txt`
+
+The calculation of BOLD signals from the simulated firing rates for Fig. 8 requires an installation of R and the R library `neuRosim` (<https://cran.r-project.org/web/packages/neuRosim/index.html>).
+
 ### Snakemake workflow
 
 The entire workflow from raw spike files till the final figures is defined in `Snakefile` and `Snakefile_preprocessing`. If snakemake is installed, the figures can be produced by executing `snakemake`.

figures/Schmidt2018_dyn/Snakefile

@@ -6,7 +6,7 @@ area_list = ['V1', 'V2', 'VP', 'V3', 'V3A', 'MT', 'V4t', 'V4', 'VOT', 'MSTd',
              'STPa', '46', 'AITd', 'TH']
 population_list = ['23E', '23I', '4E',  '4I', '5E', '5I', '6E', '6I']
 
-LOAD_ORIGINAL_DATA = False
+LOAD_ORIGINAL_DATA = True
 
 ORIGINAL_SIM_LABELS = {'all': ['533d73357fbe99f6178029e6054b571b485f40f6',
                                '0adda4a542c3d5d43aebf7c30d876b6c5fd1d63e',
@@ -55,13 +55,17 @@ ORIGINAL_SIM_LABELS = {'all': ['533d73357fbe99f6178029e6054b571b485f40f6',
 
 if LOAD_ORIGINAL_DATA:
     DATA_DIR = original_data_path
-    SIM_LABELS = ORIGINAL_SIM_LABELS_LABELS
+    SIM_LABELS = ORIGINAL_SIM_LABELS
 else:
     from network_simulations import create_label_dict
     from config import data_path
     DATA_DIR = data_path
     SIM_LABELS = create_label_dict()
 
+if DATA_DIR is None:
+    raise TypeError("The path to the data files is None. "
+                    "Please define the data path.")
+
 rule all:
     input:
         'Fig1_model_overview.eps',

figures/Schmidt2018_dyn/additional_requirements.txt

 statsmodels
-R
 networkx
 pyx
 python-louvain
+-e git+https://github.com/mschmidt87/correlation-toolbox.git@enh/add_setup_file#egg=correlation-toolbox

figures/Schmidt2018_dyn/plotstyle.rc

+# set default backend
+backend : GTK3Agg
+
+# resolution of figures in dpi
+# does not influence eps output
+figure.dpi : 150
+savefig.dpi : 300
+
+# set default figuresize
+figure.figsize : 12.0, 7.42
+
+# font
+font.size : 8
+legend.fontsize : 8
+font.family : sans-serif
+
+# size of lines
+lines.linewidth : 1.5
+lines.antialiased : True
+
+# size of markers (points in point plots)
+lines.markersize : 3.0
+patch.linewidth : 1.0
+axes.linewidth : 1.0 # edge linewidth
+
+# ticks distances
+xtick.major.size : 4 # major tick size in points
+xtick.minor.size : 2 # minor tick size in points
+lines.markeredgewidth : 0.5 # line width of ticks
+xtick.major.pad : 4 # distance to major tick label in points
+xtick.minor.pad : 4 # distance to the minor tick label in points
+ytick.major.size : 4  # major tick size in points
+ytick.minor.size : 2 # minor tick size in points
+ytick.major.pad : 4 # distance to major tick label in points
+ytick.minor.pad : 4 # distance to the minor tick label in points
+
+# ticks textsize
+ytick.labelsize : 8
+xtick.labelsize : 8
+
+# axes labelsize
+axes.labelsize : 8
+
+# use latex to generate the labels in plots
+# not needed anymore in newer versions
+# using this, font detection fails on adobe illustrator 2010-07-20
+text.usetex : True
+text.latex.preamble : \usepackage{amsmath}
+ps.useafm : False # use of afm fonts, results in small files
+ps.fonttype : 3 # Output Type 3 (Type3) or Type 42 (TrueType)
+
+# set different default color cycle
+axes.color_cycle : 4c72b0, 55a868, c44e52, 8172b2, ccb974, 64b5cd
\ No newline at end of file

figures/SchueckerSchmidt2017/Fig3_bistability.py

@@ -61,8 +61,8 @@ network_params = {'N_scaling': 1.,
                   'connection_params': conn_params,
                   'neuron_params': neuron_params,
                   'input_params': input_params}
-
-M_LA = MultiAreaModel(network_params, simulation=True,
+M_LA = MultiAreaModel(network_params,
+                      simulation=True,
                       sim_spec=sim_params,
                       analysis=True)
 M_LA.analysis.create_pop_rates()

figures/SchueckerSchmidt2017/Snakefile

+rule all:
+    input:
+        'iteration_1/K_prime.npy',
+        'iteration_2/K_prime.npy',
+        'iteration_3/K_prime.npy',
+        'iteration_4/K_prime.npy',
+        'Fig1_scheme.eps',
+#        'Fig2_EE_example.eps',
+        'Fig3_bistability.eps',
+        'Fig4_meanfield_mam.eps',
+        'Fig5_eigenspace.eps',
+        'Fig6_unstable_FP.eps',
+        'Fig7_stabilization_analysis.eps',
+        'Fig8_conn_analysis.eps',
+        'Fig9_simulation.eps'
+
+rule Stabilization:
+    output:
+        'iteration_1/K_prime.npy',
+        'iteration_2/K_prime.npy',
+        'iteration_3/K_prime.npy',
+        'iteration_4/K_prime.npy'
+    shell:
+        'python3 stabilization.py'
+
+rule Fig2_EE_example:
+    output:
+        'Fig2_EE_example.eps'
+    shell:
+        'python3 Fig2_EE_example.py'
+
+rule Fig3_bistability:
+    output:
+        'Fig3_bistability.eps'
+    shell:
+        'python3 Fig3_bistability.py'
+
+rule Fig4_meanfield_mam:
+    output:
+        'Fig4_meanfield_mam.eps'
+    shell:
+        'python3 Fig4_meanfield_mam.py'
+
+rule Fig5_eigenspace:
+    output:
+        'Fig5_eigenspace.eps'
+    shell:
+        'python3 Fig5_eigenspace.py'
+
+rule Fig6_unstable_FP:
+    output:
+        'Fig6_unstable_FP.eps'
+    shell:
+        'python3 Fig6_unstable_FP.py'
+
+rule Fig7_stabilization_analysis:
+    output:
+        'Fig7_stabilization_analysis.eps'
+    shell:
+        'python3 Fig7_stabilization_analysis.py'
+
+rule Fig8_conn_analysis:
+    output:
+        'Fig8_conn_analysis.eps'
+    shell:
+        'python3 Fig8_conn_analysis.py'
+
+        
+rule Fig9_simulation:
+    output:
+        'Fig9_simulation.eps'
+    shell:
+        'python3 Fig9_simulation.py'
+
+        

figures/SchueckerSchmidt2017/stabilization.py

@@ -5,8 +5,6 @@ from multiarea_model import MultiAreaModel
 from multiarea_model.multiarea_helpers import create_vector_mask
 from multiarea_model.stabilize import stabilize
 import utils
-import seaborn as sns
-cp = sns.color_palette()
 
 """
 Initialization
@@ -34,12 +32,12 @@ M_target = MultiAreaModel(network_params_target, theory=True,
                           theory_spec=theory_params)
 
 THREADS = 4
-load_list = [1, 2, 3, 4]
+load_list = []
 
 # This list defines which of the detected minima of the velocity
 # vector is identified as the unstable fixed point. It has to be
 # created manually.
-ind_list = [1, 1, 0, 1]
+ind_list = []
 
 """
 Main loop
@@ -68,7 +66,8 @@ for iteration in [1, 2, 3, 4, 5]:
         # Scan parameter space to find a good approximation of the
         # critical parameter value where the model crosses the
         # separatrix for the initial condition of zero rates
-        if iteration < 5: # For iteration 5, we just analyze the behavior without performing the stabilization
+        if iteration < 5:
+            # For iteration 5, we just analyze the behavior without performing the stabilization
             data[iteration] = utils.compute_iteration(7, fac_nu_ext_5E_list,
                                                       theory_params, M_base, threads=THREADS)
         else:

figures/SchueckerSchmidt2017/utils.py

+import copy
+import multiprocessing as mp
+import numpy as np
+import os
+import pylab as pl
+from functools import partial
+
+from multiarea_model.multiarea_helpers import create_mask
+from scipy.signal import find_peaks_cwt
+
+
+def integrate(f5, M_base):
+    mask5 = create_mask(M_base.structure, target_pops=['5E'],
+                        source_areas=[], external=True)
+    mask6 = create_mask(M_base.structure, target_pops=['6E'],
+                        source_areas=[], external=True)
+
+    f6 = 10/3.*f5-7/3.
+    conn_params = copy.deepcopy(M_base.params['connection_params'])
+    conn_params.update({'fac_nu_ext_5E': f5,
+                        'fac_nu_ext_6E': f6})
+    M = copy.deepcopy(M_base)
+    M.K_matrix[mask5] *= f5
+    M.K_matrix[mask6] *= f6
+    p, r = M.theory.integrate_siegert()
+    return r
+
+
+def iteration_results(fac_nu_ext_5E_list, M_base, threads=None):
+    if threads is None:
+        results = np.array([integrate(f5, M_base) for f5 in fac_nu_ext_5E_list])
+    else:
+        integrate_part = partial(integrate,
+                                 M_base=M_base)
+        pool = mp.Pool(processes=threads)
+        results = np.array(pool.map(integrate_part, fac_nu_ext_5E_list))
+        pool.close()
+        pool.join()
+    return results
+
+
+def velocity_peaks(time, result, threshold=0.05):
+    d_nu = np.abs(np.diff(np.mean(result, axis=0)) / np.diff(time)[0])
+    ind = np.where(d_nu < threshold)
+    minima = find_peaks_cwt(-1. * np.log(d_nu[ind]), np.array([0.1]))
+    if len(minima) > 0:
+        t_min = time[ind][minima]
+    else:
+        t_min = []
+    min_full = [np.argmin(np.abs(time - t)) for t in t_min]
+    return d_nu, min_full
+
+
+def plot_iteration(results, theory_params, threshold=0.05, full=True):
+    traj = np.mean(results, axis=1)
+    if full:
+        ind = list(range(0, len(traj)))
+    else:
+        i = np.argmax(np.diff(traj[:, -1]))
+        ind = [i, i+1]
+
+    time = np.arange(0., theory_params['T'], theory_params['dt'])
+    fig = pl.figure()
+    ax = fig.add_subplot(121)
+    [ax.plot(time,
+             traj[i]) for i in ind]
+
+    ax.set_yscale('Log')
+    ax.set_xlabel('Time (a.u.)')
+    ax.set_ylabel(r'$\langle \nu \rangle$')
+
+    ax = fig.add_subplot(122)
+    for n, i in enumerate(ind):
+        d_nu, minima = velocity_peaks(time, results[i], threshold=threshold)
+        ax.plot(time[:-1], d_nu)
+        if not full:
+            ax.vlines(time[minima],
+                      1e-6,
+                      1e0,
+                      linestyles='dashed',
+                      color='k')
+    ax.set_yscale('Log')
+    pl.show()
+
+    
+def save_iteration(step, data):
+    data_dir = 'iteration_{}'.format(step)
+    try:
+        os.mkdir(data_dir)
+    except FileExistsError:
+        pass
+    for key in ['parameters', 'K_prime', 'results']:
+        np.save('iteration_{}/{}.npy'.format(step, key), data[key])
+
+
+def load_iteration(step):
+    data_dir = 'iteration_{}'.format(step)
+    data = {}
+    files = os.listdir(data_dir)
+    for f in files:
+        data[os.path.splitext(f)[0]] = np.load(os.path.join(data_dir, f))
+    return data
+
+
+def compute_iteration(max_iter, fac_nu_ext_5E_list, theory_params, M_base, threads=None):
+    par_list = np.zeros(0)
+    results = np.zeros((0, 254, int(theory_params['T'] / theory_params['dt'])))
+
+    i = 0
+    while i < max_iter:
+        print("Iteration: {}".format(i))
+        print(fac_nu_ext_5E_list)
+        r = iteration_results(fac_nu_ext_5E_list, M_base, threads=threads)
+        results = np.vstack((results, r))
+        par_list = np.append(par_list, fac_nu_ext_5E_list)
+        i = np.argmax(np.diff(np.mean(r, axis=1)[:, -1]))
+        i += 1
+        fac_nu_ext_5E_list = np.arange(fac_nu_ext_5E_list[i],
+                                       # to ensure that the array includes the last value, we add a small epsilon
+                                       fac_nu_ext_5E_list[i+1] + 1.e-10,
+                                       10**(-(i+2.)))
+
+    ind = np.argsort(par_list)
+    par_list = par_list[ind]
+    results = results[ind]
+
+    return {'results': results, 'parameters': par_list}
+
+
+def determine_velocity_minima(time, data, threshold=0.05):
+    
+    traj = np.mean(data['results'], axis=1)
+    i = np.argmax(np.diff(traj[:, -1]))
+    r_low = data['results'][i]
+    r_high = data['results'][i+1]
+
+    dnu_low, minima_low = velocity_peaks(time, r_low, threshold=threshold)
+    dnu_high, minima_high = velocity_peaks(time, r_high, threshold=threshold)
+
+    par_transition = data['parameters'][i]
+    return par_transition, r_low, r_high, minima_low, minima_high
+
+    
+

multiarea_model/multiarea_helpers.py

@@ -529,9 +529,7 @@ def extract_area_dict(d, structure, target_area, source_area):
 
 def convert_syn_weight(W, neuron_params):
     """
-    # TODO: this documentation is wrong. this is not the integral.
-    Compute the integral of the PSP of exponential
-    post-synaptic currents evoked by a PSC with amplitude W.
+    Convert the amplitude of the PSC into mV.
 
     Parameters
     ----------

requirements.txt

+future
+numpy
+matplotlib
+pandas
+scipy
+nested_dict
+dicthash

run_example.py

+import numpy as np
 import os
 
 from multiarea_model import MultiAreaModel
@@ -31,12 +32,12 @@ neuron_params = {'V0_mean': -150.,
 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,
-              'input_params': input_params,
               'recording_dict': {'record_vm': False}}
 
 theory_params = {'dt': 0.1}
@@ -46,6 +47,8 @@ M = MultiAreaModel(network_params, simulation=True,
                    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)
 
 
@@ -71,13 +74,13 @@ neuron_params = {'V0_mean': -150.,
 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,
-              'input_params': input_params,
               'recording_dict': {'record_vm': False}}
 
 theory_params = {'dt': 0.1}
@@ -87,4 +90,6 @@ M = MultiAreaModel(network_params, simulation=True,
                    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])))
 M.simulation.simulate()

tests/test_network_scaling.py

@@ -4,79 +4,59 @@ import numpy as np
 import json
 
 
-#def test_network_scaling():
-"""
-Test the downscaling option of the network.
-
-- Test whether indegrees and neuron number are correctly scaled down.
-- Test whether the resulting mean and variance of the input currents
-  as well as the resulting rates are identical, based on mean-field theory.
-"""
-
-network_params = {}
-M0 = MultiAreaModel(network_params, theory=True)
-K0 = M0.K_matrix
-W0 = M0.W_matrix
-N0 = M0.N_vec
-syn0 = M0.syn_matrix
-p, r0 = M0.theory.integrate_siegert()
-
-d = vector_to_dict(r0[:, -1],
-                   M0.area_list,
-                   M0.structure)
-
-with open('mf_rates.json', 'w') as f:
-    json.dump(d, f)
-
-network_params = {'N_scaling': .1,
-                  'K_scaling': .1,
-                  'fullscale_rates': 'mf_rates.json'}
-theory_params = {'initial_rates': r0[:, -1],
-                 'T': 50.}
-M = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
-
-K = M.K_matrix
-W = M.W_matrix
-N = M.N_vec
-syn = M.syn_matrix
-p, r = M.theory.integrate_siegert()
-
-assert(np.allclose(K, network_params['K_scaling'] * K0))
-assert(np.allclose(N, network_params['N_scaling'] * N0))
-assert(np.allclose(syn, network_params['K_scaling'] * network_params['N_scaling'] * syn0))
-assert(np.allclose(W, W0 / np.sqrt(network_params['K_scaling'])))
-
-r0_extend = np.append(r0[:, -1], M0.params['input_params']['rate_ext'])
-tau_m = M.params['neuron_params']['single_neuron_dict']['tau_m']
-C_m = M.params['neuron_params']['single_neuron_dict']['C_m']
-
-mu0 = (1e-3 * tau_m * np.dot(M0.K_matrix * M0.J_matrix, r0_extend)
-       + tau_m / C_m * M0.add_DC_drive)
-mu = 1e-3 * tau_m * np.dot(M.K_matrix * M.J_matrix, r0_extend) + tau_m / C_m * M.add_DC_drive
-
-sigma0 = np.sqrt(1e-3 * tau_m * np.dot(M0.K_matrix * M0.J_matrix**2, r0_extend))
-sigma = np.sqrt(1e-3 * tau_m * np.dot(M.K_matrix * M.J_matrix**2, r0_extend))
-
-
-r0_alex_loaded = np.load('test_network_scaling_r0.npy')
-r_alex_loaded = np.load('test_network_scaling_r.npy')
-
-network_params = {}
-theory_params = {'initial_rates': r0_alex_loaded[:, -1],
-                 'T': 50.}
-M0_alex = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
-p, r0_alex = M0_alex.theory.integrate_siegert()
-
-network_params = {}
-theory_params = {'initial_rates': r_alex_loaded[:, -1],
-                 'T': 50.}
-M0_alex = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
-p, r_alex = M0_alex.theory.integrate_siegert()
-
-
-# p, r0_py = M0.theory.integrate_siegert_python(parallel=False)
-# p, r_py = M.theory.integrate_siegert_python(parallel=False)
-
-assert(np.allclose(mu, mu0))
-assert(np.allclose(sigma, sigma0))
-assert(np.allclose(r[:, -1], r0[:, -1]))
+def test_network_scaling():
+    """
+    Test the downscaling option of the network.
+
+    - Test whether indegrees and neuron number are correctly scaled down.
+    - Test whether the resulting mean and variance of the input currents
+      as well as the resulting rates are identical, based on mean-field theory.
+    """
+
+    network_params = {}
+    M0 = MultiAreaModel(network_params, theory=True)
+    K0 = M0.K_matrix
+    W0 = M0.W_matrix
+    N0 = M0.N_vec
+    syn0 = M0.syn_matrix
+    p, r0 = M0.theory.integrate_siegert()
+
+    d = vector_to_dict(r0[:, -1],
+                       M0.area_list,
+                       M0.structure)
+
+    with open('mf_rates.json', 'w') as f:
+        json.dump(d, f)
+
+    network_params = {'N_scaling': .1,
+                      'K_scaling': .1,
+                      'fullscale_rates': 'mf_rates.json'}
+    theory_params = {'initial_rates': r0[:, -1],
+                     'T': 50.}
+    M = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
+
+    K = M.K_matrix
+    W = M.W_matrix
+    N = M.N_vec
+    syn = M.syn_matrix
+    p, r = M.theory.integrate_siegert()
+
+    assert(np.allclose(K, network_params['K_scaling'] * K0))
+    assert(np.allclose(N, network_params['N_scaling'] * N0))
+    assert(np.allclose(syn, network_params['K_scaling'] * network_params['N_scaling'] * syn0))
+    assert(np.allclose(W, W0 / np.sqrt(network_params['K_scaling'])))
+
+    r0_extend = np.append(r0[:, -1], M0.params['input_params']['rate_ext'])
+    tau_m = M.params['neuron_params']['single_neuron_dict']['tau_m']
+    C_m = M.params['neuron_params']['single_neuron_dict']['C_m']
+
+    mu0 = (1e-3 * tau_m * np.dot(M0.K_matrix * M0.J_matrix, r0_extend)
+           + tau_m / C_m * M0.add_DC_drive)
+    mu = 1e-3 * tau_m * np.dot(M.K_matrix * M.J_matrix, r0_extend) + tau_m / C_m * M.add_DC_drive
+
+    sigma0 = np.sqrt(1e-3 * tau_m * np.dot(M0.K_matrix * M0.J_matrix**2, r0_extend))
+    sigma = np.sqrt(1e-3 * tau_m * np.dot(M.K_matrix * M.J_matrix**2, r0_extend))
+
+    assert(np.allclose(mu, mu0))
+    assert(np.allclose(sigma, sigma0))
+    assert(np.allclose(r[:, -1], r0[:, -1]))