diff --git a/punit_responses.py b/punit_responses.py
index 68ec34b..f360bb6 100644
--- a/punit_responses.py
+++ b/punit_responses.py
@@ -1,5 +1,6 @@
 import numpy as np
 import nixio as nix
+import argparse
 import os
 from numpy.core.fromnumeric import repeat
 from traitlets.traitlets import Instance
@@ -183,11 +184,11 @@ def simulate_responses(stimulus_params, model_params, repeats=10, deltaf=20):
     print("\n")
 
 
-def simulate_cell(cell_id, models):
-    deltafs = [-200, -100, -50, -20, -10, -5, 5, 10, 20, 50, 100, 200]  # Hz, difference frequency between self and other
-    chirp_sizes = [40, 60, 100]
+def simulate_cell(cell_id, models, args):
+    deltafs =  args.deltafs  # Hz, difference frequency between self and other
+    chirp_sizes = args.chirpsizes  # Hz, the chirp size, i.e. the frequency excursion 
     stimulus_params = { "eodfs": {"self": 0.0, "other": 0.0}, # eod frequency in Hz, to be overwritten
-                        "contrasts": [20, 10, 5, 2.5, 1.25, 0.625, 0.3125],
+                        "contrasts": args.contrasts,
                         "chirp_size": 100,  # Hz, frequency excursion
                         "chirp_duration": 0.015,  # s, chirp duration
                         "chirp_amplitude_dip": 0.05,  # %, amplitude drop during chirp
@@ -212,14 +213,32 @@ def simulate_cell(cell_id, models):
                                     stimulus_params["duration"] - stimulus_params["chirp_duration"], 
                                     1./stimulus_params["chirp_frequency"])
             stimulus_params["chirp_times"] = chirp_times     
-            simulate_responses(stimulus_params, model_params, repeats=25, deltaf=deltaf)
+            simulate_responses(stimulus_params, model_params, repeats=args.trials, deltaf=deltaf)
 
 
 def main():
+    parser = argparse.ArgumentParser(description="Simulate P-unit responses using the model parameters from the models.csv file. Calling it without any arguments works with the defaults, may need some time.")
+    parser.add_argument("-n", "--number", type=int, default=20, help="Number of simulated neurons. Randomly chosen from model list. Defaults to 20")
+    parser.add_argument("-t", "--trials", type=int, default=25, help="Number of stimulus repetitions, trials. Defaults to 25")
+    parser.add_argument("-dfs", "--deltafs", type=float, nargs="+", default=[-200, -100, -50, -20, -10, -5, 5, 10, 20, 50, 100, 200],
+                        help="List of difference frequencies. Defaults to [-200, -100, -50, -20, -10, -5, 5, 10, 20, 50, 100, 200]")
+    parser.add_argument("-cs", "--chirpsizes", type=float, nargs="+", default=[40, 60, 100],
+                        help="List of chirp sizes. Defaults to [40, 60, 100]")
+    parser.add_argument("-ct", "--contrasts", type=float, nargs="+", default=[20, 10, 5, 2.5, 1.25, 0.625, 0.3125],
+                        help="List of foreign fish contrasts. Defaults to [20, 10, 5, 2.5, 1.25, 0.625, 0.3125].")
+    parser.add_argument("-o", "--output_folder", type=str, default=data_folder, help="Where to store the data. Defaults to %s"%os.path.join(".", data_folder))
+    parser.add_argument("-j", "--jobs", type=int, default=max(1, int(np.floor(multiprocessing.cpu_count() * 0.5))), help="Number of parallel processes (simulations) defaults to half of the available cores.")
+    args = parser.parse_args()
+
     models = load_models("models.csv")
-    num_cores = multiprocessing.cpu_count() - 6
+    num_models = args.number
+    if args.number > len(models):
+        print("INFO: number of cells larger than number of available models. Reset to max number of models.")
+        num_models = len(models)        
+    indices = list(range(len(models)))
+    np.random.shuffle(indices)
     
-    Parallel(n_jobs=num_cores)(delayed(simulate_cell)(cell_id, models) for cell_id in range(len(models[:20])))       
+    Parallel(n_jobs=args.jobs)(delayed(simulate_cell)(cell_id, models, args) for cell_id in indices)    
 
 
 if __name__ == "__main__":