diff --git a/ASTIM.py b/ASTIM.py
index c4229fd..4b9150a 100644
--- a/ASTIM.py
+++ b/ASTIM.py
@@ -1,136 +1,141 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # @Author: Theo Lemaire
 # @Date:   2017-02-13 18:16:09
 # @Email: theo.lemaire@epfl.ch
 # @Last Modified by:   Theo Lemaire
-# @Last Modified time: 2018-09-23 15:27:47
+# @Last Modified time: 2018-09-23 16:11:08
 
 ''' Run A-STIM simulations of a specific point-neuron. '''
 
 import os
 import logging
 import numpy as np
 from argparse import ArgumentParser
 
 from PySONIC.core import NeuronalBilayerSonophore
 from PySONIC.utils import logger, checkBatchLog, selectDirDialog
 from PySONIC.neurons import getNeuronsDict
 from PySONIC.batches import createSimQueue, runBatch
 from PySONIC.plt import plotBatch
 
 # Default parameters
 defaults = dict(
     neuron='RS',
     diams=[32.0],  # nm
     freqs=[500.0],  # kHz
     amps=[100.0],  # kPa
     durations=[100.0],  # ms
     PRFs=[100.0],  # Hz
     DCs=[100.0],  # %
     offsets=[50.],  # ms
     method='sonic'
 )
 
 
 def runAStimBatch(outdir, logpath, nbls, stim_params, method, mpi=False):
     ''' Run batch A-STIM simulations of the system for various neuron types and
         stimulation parameters.
 
         :param outdir: full path to output directory
         :param logpath: full path log file
         :param stim_params: dictionary containing sweeps for all stimulation parameters
         :param method: numerical integration method ("classic", "hybrid" or "sonic")
         :param mpi: boolean statting wether or not to use multiprocessing
         :return: list of full paths to the output files
     '''
 
     mandatory_params = ['freqs', 'durations', 'offsets', 'PRFs', 'DCs']
     for mparam in mandatory_params:
         if mparam not in stim_params:
             raise ValueError('Missing stimulation parameter field: "{}"'.format(mparam))
 
     logger.info("Starting A-STIM simulation batch")
 
     # Generate queue
     nofreq_queue = createSimQueue(stim_params.get('amps', [None]), stim_params['durations'],
                                   stim_params['offsets'], stim_params['PRFs'], stim_params['DCs'])
 
     # Repeat queue for each US frequency
     nofreq_queue = np.array(nofreq_queue)
     freqs = stim_params['freqs']
     nf = len(freqs)
     nqueue = nofreq_queue.shape[0]
     queue = np.tile(nofreq_queue, (nf, 1))
     freqs_col = np.vstack([np.ones(nqueue) * f for f in freqs]).reshape(nf * nqueue, 1)
     queue = np.hstack((freqs_col, queue)).tolist()
 
     # Add method to queue items
     for item in queue:
         item.append(method)
 
     # Run batch
     return runBatch(nbls, 'runAndSave', queue, extra_params=[outdir, logpath], mpi=mpi)
 
 
-if __name__ == '__main__':
+def main():
 
     ap = ArgumentParser()
 
     # Runtime options
     ap.add_argument('--mpi', default=False, action='store_true', help='Use multiprocessing')
     ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity')
     ap.add_argument('-p', '--plot', default=False, action='store_true', help='Plot results')
     ap.add_argument('-o', '--outputdir', type=str, default=None, help='Output directory')
     ap.add_argument('-t', '--titrate', default=False, action='store_true', help='Perform titration')
     ap.add_argument('-m', '--method', type=str, default=defaults['method'],
                     help='Numerical integration method ("classic", "hybrid" or "sonic")')
 
     # Stimulation parameters
     ap.add_argument('-n', '--neuron', type=str, default=defaults['neuron'],
                     help='Neuron name (string)')
     ap.add_argument('-a', '--diams', nargs='+', type=float, help='Sonophore diameter (nm)')
     ap.add_argument('-f', '--freqs', nargs='+', type=float, help='US frequency (kHz)')
     ap.add_argument('-A', '--amps', nargs='+', type=float, help='Acoustic pressure amplitude (kPa)')
     ap.add_argument('-d', '--durations', nargs='+', type=float, help='Stimulus duration (ms)')
     ap.add_argument('--offset', nargs='+', type=float, help='Offset duration (ms)')
     ap.add_argument('--PRF', nargs='+', type=float, help='PRF (Hz)')
     ap.add_argument('--DC', nargs='+', type=float, help='Duty cycle (%%)')
 
     # Parse arguments
     args = {key: value for key, value in vars(ap.parse_args()).items() if value is not None}
     loglevel = logging.DEBUG if args['verbose'] is True else logging.INFO
     logger.setLevel(loglevel)
     outdir = args['outputdir'] if 'outputdir' in args else selectDirDialog()
     mpi = args['mpi']
     plot = args['plot']
     titrate = args['titrate']
     method = args['method']
     neuron_str = args['neuron']
     diams = np.array(args.get('diams', defaults['diams'])) * 1e-9  # m
     stim_params = dict(
         freqs=np.array(args.get('freqs', defaults['freqs'])) * 1e3,  # Hz
         amps=np.array(args.get('amps', defaults['amps'])) * 1e3,  # Pa
         durations=np.array(args.get('durations', defaults['durations'])) * 1e-3,  # s
         PRFs=np.array(args.get('PRFs', defaults['PRFs'])),  # Hz
         DCs=np.array(args.get('DCs', defaults['DCs'])) * 1e-2,  # (-)
         offsets=np.array(args.get('offsets', defaults['offsets'])) * 1e-3  # s
     )
     if titrate:
         stim_params['amps'] = [None]
 
     # Run A-STIM batch
     logpath, _ = checkBatchLog(outdir, 'A-STIM')
     if neuron_str not in getNeuronsDict():
-        raise ValueError('Unknown neuron type: "{}"'.format(neuron_str))
+        logger.error('Unknown neuron type: "%s"', neuron_str)
+        return
     neuron = getNeuronsDict()[neuron_str]()
     pkl_filepaths = []
     for a in diams:
         nbls = NeuronalBilayerSonophore(a, neuron)
         pkl_filepaths += runAStimBatch(outdir, logpath, nbls, stim_params, method, mpi=mpi)
     pkl_dir, _ = os.path.split(pkl_filepaths[0])
 
     # Plot resulting profiles
     if plot:
         yvars = {'Q_m': ['Qm']}
         plotBatch(pkl_dir, pkl_filepaths, yvars)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/ESTIM.py b/ESTIM.py
index bd69d3d..f397bc9 100644
--- a/ESTIM.py
+++ b/ESTIM.py
@@ -1,106 +1,110 @@
 # -*- coding: utf-8 -*-
 # @Author: Theo Lemaire
 # @Date:   2017-08-24 11:55:07
 # @Last Modified by:   Theo Lemaire
-# @Last Modified time: 2018-09-23 14:32:39
+# @Last Modified time: 2018-09-23 16:11:57
 
 ''' Run E-STIM simulations of a specific point-neuron. '''
 
 import os
 import logging
 import numpy as np
 from argparse import ArgumentParser
 
 from PySONIC.utils import logger, checkBatchLog, selectDirDialog
 from PySONIC.neurons import *
 from PySONIC.batches import createSimQueue, runBatch
 from PySONIC.plt import plotBatch
 
 # Default parameters
 defaults = dict(
     neuron='RS',
     amps=[10.0],  # mA/m2
     durations=[100.0],  # ms
     PRFs=[100.0],  # Hz
     DCs=[100.0],  # %
     offsets=[50.],  # ms
     method='sonic'
 )
 
 
 def runEStimBatch(outdir, logpath, neuron, stim_params, mpi=False):
     ''' Run batch E-STIM simulations of the system for various neuron types and
         stimulation parameters.
 
         :param outdir: full path to output directory
         :param logpath: full path log file
         :param stim_params: dictionary containing sweeps for all stimulation parameters
         :param mpi: boolean statting wether or not to use multiprocessing
         :return: list of full paths to the output files
     '''
     mandatory_params = ['durations', 'offsets', 'PRFs', 'DCs']
     for mparam in mandatory_params:
         if mparam not in stim_params:
             raise ValueError('Missing stimulation parameter field: "{}"'.format(mparam))
 
     logger.info("Starting E-STIM simulation batch")
 
     # Generate simulations queue
     queue = createSimQueue(stim_params.get('amps', [None]), stim_params['durations'],
                            stim_params['offsets'], stim_params['PRFs'], stim_params['DCs'])
 
     # Run batch
     return runBatch(neuron, 'runAndSave', queue, extra_params=[outdir, logpath], mpi=mpi)
 
 
-if __name__ == '__main__':
-
+def main():
     ap = ArgumentParser()
 
     # Runtime options
     ap.add_argument('--mpi', default=False, action='store_true', help='Use multiprocessing')
     ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity')
     ap.add_argument('-p', '--plot', default=False, action='store_true', help='Plot results')
     ap.add_argument('-o', '--outputdir', type=str, default=None, help='Output directory')
     ap.add_argument('-t', '--titrate', default=False, action='store_true', help='Perform titration')
 
     # Stimulation parameters
     ap.add_argument('-n', '--neuron', type=str, default=defaults['neuron'],
                     help='Neuron name (string)')
     ap.add_argument('-A', '--amps', nargs='+', type=float, help='Injected current density (mA/m2)')
     ap.add_argument('-d', '--durations', nargs='+', type=float, help='Stimulus duration (ms)')
     ap.add_argument('--offset', nargs='+', type=float, help='Offset duration (ms)')
     ap.add_argument('--PRF', nargs='+', type=float, help='PRF (Hz)')
     ap.add_argument('--DC', nargs='+', type=float, help='Duty cycle (%%)')
 
     # Parse arguments
     args = {key: value for key, value in vars(ap.parse_args()).items() if value is not None}
     loglevel = logging.DEBUG if args['verbose'] is True else logging.INFO
     logger.setLevel(loglevel)
     outdir = args['outputdir'] if 'outputdir' in args else selectDirDialog()
     mpi = args['mpi']
     plot = args['plot']
     titrate = args['titrate']
     neuron_str = args['neuron']
     stim_params = dict(
         amps=np.array(args.get('amps', defaults['amps'])),  # mA/m2
         durations=np.array(args.get('durations', defaults['durations'])) * 1e-3,  # s
         PRFs=np.array(args.get('PRFs', defaults['PRFs'])),  # Hz
         DCs=np.array(args.get('DCs', defaults['DCs'])) * 1e-2,  # (-)
         offsets=np.array(args.get('offsets', defaults['offsets'])) * 1e-3  # s
     )
     if titrate:
         stim_params['amps'] = [None]
 
     # Run E-STIM batch
     logpath, _ = checkBatchLog(outdir, 'E-STIM')
     if neuron_str not in getNeuronsDict():
-        raise ValueError('Unknown neuron type: "{}"'.format(neuron_str))
+        logger.error('Unknown neuron type: "%s"', neuron_str)
+        return
     neuron = getNeuronsDict()[neuron_str]()
     pkl_filepaths = runEStimBatch(outdir, logpath, neuron, stim_params, mpi=mpi)
     pkl_dir, _ = os.path.split(pkl_filepaths[0])
 
     # Plot resulting profiles
     if plot:
         yvars = {'V_m': ['Vm']}
         plotBatch(pkl_dir, pkl_filepaths, yvars)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/MECH.py b/MECH.py
index bcbaba0..c008214 100644
--- a/MECH.py
+++ b/MECH.py
@@ -1,112 +1,115 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # @Author: Theo Lemaire
 # @Date:   2016-11-21 10:46:56
 # @Email: theo.lemaire@epfl.ch
 # @Last Modified by:   Theo Lemaire
-# @Last Modified time: 2018-09-23 13:03:02
+# @Last Modified time: 2018-09-23 16:12:32
 
 ''' Run simulations of the NICE mechanical model. '''
 
 import os
 import logging
 import numpy as np
 from argparse import ArgumentParser
 
 from PySONIC.core import BilayerSonophore
 from PySONIC.utils import logger, selectDirDialog, checkBatchLog
 from PySONIC.neurons import CorticalRS
 from PySONIC.batches import createQueue, runBatch
 from PySONIC.plt import plotBatch
 
 # Default parameters
 defaults = dict(
     Cm0=CorticalRS().Cm0 * 1e2,  # uF/m2
     Qm0=CorticalRS().Vm0,  # nC/m2
     diams=[32.0],  # nm
     embeddings=[0.],  # um
     freqs=[500.0],  # kHz
     amps=[100.0],  # kPa
     charges=[0.]  # nC/cm2
 )
 
 
 def runMechBatch(outdir, logpath, bls, stim_params, mpi=False):
     ''' Run batch simulations of the mechanical system with imposed values of charge density.
 
         :param outdir: full path to output directory
         :param logpath: full path log file
         :param bls: BilayerSonophore object
         :param stim_params: dictionary containing sweeps for all stimulation parameters
         :param mpi: boolean stating whether or not to use multiprocessing
         :return: list of full paths to the output files
     '''
 
     # Checking validity of stimulation parameters
     mandatory_params = ['freqs', 'amps', 'charges']
     for mparam in mandatory_params:
         if mparam not in stim_params:
             raise ValueError('Missing stimulation parameter field: "{}"'.format(mparam))
 
     logger.info("Starting mechanical simulation batch")
 
     # Unpack stimulation parameters
     freqs = np.array(stim_params['freqs'])
     amps = np.array(stim_params['amps'])
     charges = np.array(stim_params['charges'])
 
     # Generate simulations queue and run batch
     queue = createQueue((freqs, amps, charges))
     return runBatch(bls, 'runAndSave', queue, extra_params=[outdir, logpath], mpi=mpi)
 
 
-if __name__ == '__main__':
-
+def main():
     ap = ArgumentParser()
 
     # Runtime options
     ap.add_argument('--mpi', default=False, action='store_true', help='Use multiprocessing')
     ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity')
     ap.add_argument('-p', '--plot', default=False, action='store_true', help='Plot results')
     ap.add_argument('-o', '--outputdir', type=str, default=None, help='Output directory')
 
     # Stimulation parameters
     ap.add_argument('-a', '--diams', nargs='+', type=float, help='Sonophore diameter (nm)')
     ap.add_argument('--Cm0', type=float, default=defaults['Cm0'],
                     help='Resting membrane capacitance (uF/cm2)')
     ap.add_argument('--Qm0', type=float, default=defaults['Qm0'],
                     help='Resting membrane charge density (nC/cm2)')
     ap.add_argument('-d', '--embeddings', nargs='+', type=float, help='Embedding depth (um)')
     ap.add_argument('-f', '--freqs', nargs='+', type=float, help='US frequency (kHz)')
     ap.add_argument('-A', '--amps', nargs='+', type=float, help='Acoustic pressure amplitude (kPa)')
     ap.add_argument('-Q', '--charges', nargs='+', type=float, help='Membrane charge density (nC/cm2)')
 
     # Parse arguments
     args = {key: value for key, value in vars(ap.parse_args()).items() if value is not None}
     loglevel = logging.DEBUG if args['verbose'] is True else logging.INFO
     logger.setLevel(loglevel)
     outdir = args['outputdir'] if 'outputdir' in args else selectDirDialog()
     mpi = args['mpi']
     plot = args['plot']
     Cm0 = args['Cm0'] * 1e-2  # F/m2
     Qm0 = args['Qm0'] * 1e-5  # C/m2
     diams = np.array(args.get('diams', defaults['diams'])) * 1e-9  # m
     embeddings = np.array(args.get('embeddings', defaults['embeddings'])) * 1e-6  # m
     stim_params = dict(
         freqs=np.array(args.get('freqs', defaults['freqs'])) * 1e3,  # Hz
         amps=np.array(args.get('amps', defaults['amps'])) * 1e3,  # Pa
         charges=np.array(args.get('charges', defaults['charges'])) * 1e-5  # C/m2
     )
 
     # Run MECH batch
     logpath, _ = checkBatchLog(outdir, 'MECH')
     pkl_filepaths = []
     for a in diams:
         for d in embeddings:
             bls = BilayerSonophore(a, Cm0, Qm0, embedding_depth=d)
             pkl_filepaths += runMechBatch(outdir, logpath, bls, stim_params, mpi=mpi)
     pkl_dir, _ = os.path.split(pkl_filepaths[0])
 
     # Plot resulting profiles
     if plot:
         plotBatch(pkl_dir, pkl_filepaths)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/README.md b/README.md
index a4c4a3c..689d1f2 100644
--- a/README.md
+++ b/README.md
@@ -1,68 +1,68 @@
 ## Description
 
 Python implementation of the **multi-Scale Optimized Neuronal Intramembrane Cavitation** (SONIC) model to compute individual neural responses to acoustic stimuli, as predicted by the *intramembrane cavitation* hypothesis.
 
 This package contains several core modules:
 - `bls` defines the underlying biomechanical model of intramembrane cavitation (`BilayerSonophore` class), and provides an integration method to predict compute the mechanical oscillations of the plasma membrane subject to a periodic acoustic perturbation.
 - `pneuron` defines a generic `PointNeuron` class that contains methods to simulate a *Hodgkin-Huxley* point-neuron model.
 - `neurons` contains classes that inherit from `PointNeuron` and define intrinsic channels mechanisms of several neurons types.
 - `nbls` defines a `NeuronalBilayerSonophore` class that inherits from `BilayerSonophore` and receives a specific `PointNeuron` child instance at initialization, and implements the bi-directional link between the mechanical and electrical parts of the NICE model. It also provides several integration methods (detailed below) to simulate the full electro-mechanical model behavior upon sonication:
 	- `runFull` solves all variables for the entire duration of the simulation. This method uses very small time steps and is computationally expensive (simulation time: several hours)
 	- `runHybrid` integrates the system by consecutive “slices” of time, during which the full system is solved until mechanical stabilization, at which point the electrical system is solely solved with predicted mechanical variables until the end of the slice. This method is more efficient (simulation time: several minutes) and provides accurate results.
 	- `runSONIC` integrates a coarse-grained, effective differential system to grasp the net effect of the acoustic stimulus on the electrical system. This method allows to run simulations of the electrical system in only a few seconds, with very accurate results of the net membrane charge density evolution, but requires the pre-computation of lookup tables.
 
 As well as some additional modules:
 - `plt` defines graphing utilities to load results of several simulations and display/compare temporal profiles of multiple variables of interest across simulations.
 - `utils` defines generic utilities used across the different modules
 
 
 ## Installation
 
 Install Python 3 if not already done.
 
 Open a terminal.
 
 Activate a Python3 environment if needed, e.g. on the tnesrv5 machine:
 
 ```$ source /opt/apps/anaconda3/bin activate```
 
 Check that the appropriate version of pip is activated:
 
 ```$ pip --version```
 
 Go to the package directory (where the setup.py file is located) and install it:
 
 ```
 $ cd <path_to_directory>
 $ pip install -e .
 ```
 
 *PySONIC* and all its dependencies will be installed.
 
 ## Usage
 
 ### Command line scripts
 
 Open a terminal at the package root directory.
 
-Use `ESTIM.py` to simulate a point-neuron model upon electrical stimulation, e.g. for a *regular-spiking neuron* injected with 10 mA/m2 intracellular current for 30 ms:
+Use `estim.py` to simulate a point-neuron model upon electrical stimulation, e.g. for a *regular-spiking neuron* injected with 10 mA/m2 intracellular current for 30 ms:
 
-```$ python ESTIM.py -n=RS -A=10 -d=30```
+```$ python estim.py -n=RS -A=10 -d=30```
 
-Use `MECH.py` to simulate mechanical model upon sonication (until periodic stabilization), e.g. for a 32 nm diameter sonophore sonicated at 500 kHz and 100 kPa:
+Use `mech.py` to simulate mechanical model upon sonication (until periodic stabilization), e.g. for a 32 nm diameter sonophore sonicated at 500 kHz and 100 kPa:
 
-```$ python MECH.py -a=32 -f=500 -A=100```
+```$ python mech.py -a=32 -f=500 -A=100```
 
-Use `ASTIM.py` to simulate the full electro-mechanical model of a given neuron type upon sonication, e.g. for a 32 nm diameter sonophore within a *regular-spiking neuron* sonicated at 500 kHz and 100 kPa for 150 ms:
+Use `astim.py` to simulate the full electro-mechanical model of a given neuron type upon sonication, e.g. for a 32 nm diameter sonophore within a *regular-spiking neuron* sonicated at 500 kHz and 100 kPa for 150 ms:
 
-```$ python ASTIM.py -n=RS -a=32 -f=500 -A=100 -d=150```
+```$ python astim.py -n=RS -a=32 -f=500 -A=100 -d=150```
 
 If several values are defined for a given parameter, a batch of simulations is run (for every value of the parameter sweep).
 You can also specify these values from within the script (*defaults* dictionary)
 
 The simulation results will be save in an output PKL file. To view these results, you can use the `-p` option
 
 Several more options are available. To view them, type in
 
 ```$ python <script_name> -h```
 
diff --git a/lookups.py b/lookups.py
index f0a2713..6333826 100644
--- a/lookups.py
+++ b/lookups.py
@@ -1,145 +1,139 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # @Author: Theo Lemaire
 # @Date:   2017-06-02 17:50:10
 # @Email: theo.lemaire@epfl.ch
 # @Last Modified by:   Theo Lemaire
-# @Last Modified time: 2018-09-23 15:27:50
+# @Last Modified time: 2018-09-23 16:16:28
 
 """ Create lookup table for specific neuron. """
 
 import os
 import pickle
 import logging
 import numpy as np
 from argparse import ArgumentParser
 
 from PySONIC.utils import logger, getLookupDir
 from PySONIC.batches import createQueue, runBatch
 from PySONIC.neurons import getNeuronsDict
 from PySONIC.core import NeuronalBilayerSonophore
 
 
 # Default parameters
-default = {
-    'neuron': 'RS',
-    'diams': np.array([16.0, 32.0, 64.0]),  # nm
-    'freqs': np.array([20., 100., 500., 1e3, 2e3, 3e3, 4e3]),  # kHz
-    'amps': np.insert(np.logspace(np.log10(0.1), np.log10(600), num=50), 0, 0.0),  # kPa
-}
+defaults = dict(
+    neuron='RS',
+    diams=np.array([16.0, 32.0, 64.0]),  # nm
+    freqs=np.array([20., 100., 500., 1e3, 2e3, 3e3, 4e3]),  # kHz
+    amps=np.insert(np.logspace(np.log10(0.1), np.log10(600), num=50), 0, 0.0)  # kPa
+)
 
 
 def computeAStimLookups(neuron, aref, fref, Aref, phi=np.pi, mpi=False, loglevel=logging.INFO):
     ''' Run simulations of the mechanical system for a multiple combinations of
         imposed US frequencies, acoustic amplitudes and charge densities, compute
         effective coefficients and store them in a dictionary of 3D arrays.
 
         :param neuron: neuron object
         :param aref: array of sonophore diameters (m)
         :param fref: array of acoustic drive frequencies (Hz)
         :param Aref: array of acoustic drive amplitudes (Pa)
         :param phi: acoustic drive phase (rad)
         :param mpi: boolean statting wether or not to use multiprocessing
         :return: lookups dictionary
     '''
 
     # Check validity of input parameters
     for key, values in {'diameters': aref, 'frequencies': fref, 'amplitudes': Aref}.items():
         if not (isinstance(values, list) or isinstance(values, np.ndarray)):
             raise TypeError('Invalid {} (must be provided as list or numpy array)'.format(key))
         if not all(isinstance(x, float) for x in values):
             raise TypeError('Invalid {} (must all be float typed)'.format(key))
         if len(values) == 0:
             raise ValueError('Empty {} array'.format(key))
         if key in ('diameters', 'frequencies') and min(values) <= 0:
             raise ValueError('Invalid {} (must all be strictly positive)'.format(key))
         if key is 'amplitudes' and min(values) < 0:
             raise ValueError('Invalid {} (must all be positive or null)'.format(key))
 
     # create simulation queue
     Qref = np.arange(neuron.Qbounds[0], neuron.Qbounds[1] + 1e-5, 1e-5)  # C/m2
     na, nf, nA, nQ = len(aref), len(fref), len(Aref), len(Qref)
     queue = createQueue((fref, Aref, Qref))
 
     # run simulations and populate outputs (list of lists)
     logger.info('Starting simulation batch for %s neuron', neuron.name)
     outputs = []
     for a in aref:
         nbls = NeuronalBilayerSonophore(a, neuron)
         outputs += runBatch(nbls, 'computeEffVars', queue, mpi=mpi, loglevel=loglevel)
     outputs = np.array(outputs).T
 
     # populate lookups dictionary with input vectors
     lookups = dict(
         a=aref,  # nm
         f=fref,  # Hz
         A=Aref,  # Pa
         Q=Qref  # C/m2
     )
 
     # reshape outputs into 4D arrays and add them to lookups dictionary
     logger.info('Reshaping output into lookup tables')
     for cname, output in zip(neuron.coeff_names, outputs):
         lookups[cname] = output.reshape(na, nf, nA, nQ)
 
     return lookups
 
 
 def main():
-
-    # Define argument parser
     ap = ArgumentParser()
 
+    # Runtime options
+    ap.add_argument('--mpi', default=False, action='store_true', help='Use multiprocessing')
+    ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity')
+
     # Stimulation parameters
-    ap.add_argument('-n', '--neuron', type=str, default=default['neuron'],
+    ap.add_argument('-n', '--neuron', type=str, default=defaults['neuron'],
                     help='Neuron name (string)')
-    ap.add_argument('-a', '--diameters', nargs='+', type=float, default=None,
-                    help='Sonophore diameters (nm)')
-    ap.add_argument('-f', '--frequencies', nargs='+', type=float, default=None,
-                    help='Acoustic drive frequencies (kHz)')
-    ap.add_argument('-A', '--amplitudes', nargs='+', type=float, default=None,
-                    help='Acoustic pressure amplitudes (kPa)')
-
-    # Boolean parameters
-    ap.add_argument('-m', '--multiprocessing', default=False, action='store_true',
-                    help='Use multiprocessing')
-    ap.add_argument('-v', '--verbose', default=False, action='store_true',
-                    help='Increase verbosity')
+    ap.add_argument('-a', '--diams', nargs='+', type=float, help='Sonophore diameter (nm)')
+    ap.add_argument('-f', '--freqs', nargs='+', type=float, help='US frequency (kHz)')
+    ap.add_argument('-A', '--amps', nargs='+', type=float, help='Acoustic pressure amplitude (kPa)')
 
     # Parse arguments
-    args = ap.parse_args()
-    neuron_str = args.neuron
-    diams = (default['diams'] if args.diameters is None else np.array(args.diameters)) * 1e-9  # m
-    freqs = (default['freqs'] if args.frequencies is None else np.array(args.frequencies)) * 1e3  # Hz
-    amps = (default['amps'] if args.amplitudes is None else np.array(args.amplitudes)) * 1e3  # Pa
-
-    loglevel = logging.DEBUG if args.verbose else logging.INFO
+    args = {key: value for key, value in vars(ap.parse_args()).items() if value is not None}
+    loglevel = logging.DEBUG if args['verbose'] is True else logging.INFO
     logger.setLevel(loglevel)
+    mpi = args['mpi']
+    neuron_str = args['neuron']
+    diams = np.array(args.get('diams', defaults['diams'])) * 1e-9  # m
+    freqs = np.array(args.get('freqs', defaults['freqs'])) * 1e3  # Hz
+    amps = np.array(args.get('amps', defaults['amps'])) * 1e3  # Pa
 
     # Check neuron name validity
     if neuron_str not in getNeuronsDict():
         logger.error('Unknown neuron type: "%s"', neuron_str)
         return
     neuron = getNeuronsDict()[neuron_str]()
 
     # Check if lookup file already exists
     lookup_file = 'SONIC_{}.pkl'.format(neuron.name)
     lookup_filepath = '{0}/{1}'.format(getLookupDir(), lookup_file)
     if os.path.isfile(lookup_filepath):
         logger.warning('"%s" file already exists and will be overwritten. ' +
                        'Continue? (y/n)', lookup_file)
         user_str = input()
         if user_str not in ['y', 'Y']:
             logger.error('%s Lookup creation canceled', neuron.name)
             return
 
-    lookup_dict = computeAStimLookups(neuron, diams, freqs, amps,
-                                      mpi=args.multiprocessing, loglevel=loglevel)
+    # compute lookups
+    lookup_dict = computeAStimLookups(neuron, diams, freqs, amps, mpi=mpi, loglevel=loglevel)
+
     # Save dictionary in lookup file
     logger.info('Saving %s neuron lookup table in file: "%s"', neuron.name, lookup_file)
     with open(lookup_filepath, 'wb') as fh:
         pickle.dump(lookup_dict, fh)
 
 
 if __name__ == '__main__':
     main()
diff --git a/setup.py b/setup.py
index a414186..4a9e0c6 100644
--- a/setup.py
+++ b/setup.py
@@ -1,51 +1,51 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # @Author: Theo Lemaire
 # @Date:   2017-06-13 09:40:02
 # @Email: theo.lemaire@epfl.ch
 # @Last Modified by:   Theo Lemaire
-# @Last Modified time: 2018-09-23 13:51:29
+# @Last Modified time: 2018-09-23 16:14:35
 
 from setuptools import setup
 
 
 def readme():
     with open('README.md', encoding="utf8") as f:
         return f.read()
 
 
 setup(
     name='PySONIC',
     version='1.0',
     description='Python implementation of the **multi-Scale Optimized Neuronal Intramembrane \
                Cavitation** (SONIC) model to compute individual neural responses to acoustic \
                stimuli, as predicted by the *intramembrane cavitation* hypothesis.',
     long_description=readme(),
     url='???',
     classifiers=[
         'Development Status :: 4 - Beta',
         'Intended Audience :: Science/Research',
         'Programming Language :: Python :: 3',
         'Topic :: Scientific/Engineering :: Physics'
     ],
     keywords=('SONIC NICE acoustic ultrasound ultrasonic neuromodulation neurostimulation excitation\
              computational model intramembrane cavitation'),
     author='Théo Lemaire',
     author_email='theo.lemaire@epfl.ch',
     license='MIT',
     packages=['PySONIC'],
-    scripts=['scripts/MECH.py', 'scripts/ESTIM.py', 'scripts/ASTIM.py'],
+    scripts=['mech.py', 'estim.py', 'astim.py'],
     install_requires=[
         'numpy>=1.10',
         'scipy>=0.17',
         'matplotlib>=2'
         'pandas>=0.22.0',
         'openpyxl>=2.4',
         'pyyaml>=3.11',
         'pycallgraph>=1.0.1',
         'colorlog>=3.0.1',
         'progressbar2>=3.18.1',
         'lockfile>=0.1.2'
     ],
     zip_safe=False
 )