diff --git a/examples/plot_simulate_evoked.py b/examples/plot_simulate_evoked.py
index 9819dfb57..15a42a723 100644
--- a/examples/plot_simulate_evoked.py
+++ b/examples/plot_simulate_evoked.py
@@ -35,7 +35,7 @@
 # Now let's simulate the dipole
 # You can simulate multiple trials in parallel by using n_jobs > 1
 net = Network(params)
-dpls = simulate_dipole(net, n_jobs=1, n_trials=2)
+dpls = simulate_dipole(net, n_jobs=3, n_trials=3)
 
 ###############################################################################
 # and then plot it
diff --git a/mne_neuron/dipole.py b/mne_neuron/dipole.py
index 99e654cf2..b9137535f 100644
--- a/mne_neuron/dipole.py
+++ b/mne_neuron/dipole.py
@@ -16,13 +16,16 @@ def _hammfilt(x, winsz):
     return convolve(x, win, 'same')
 
 
-def _simulate_single_trial(params):
+def _simulate_single_trial(params, trial_idx):
     """Simulate one trial."""
     from .parallel import rank, pc, cvode
     from .network import Network
 
     from neuron import h
     h.load_file("stdrun.hoc")
+
+    if trial_idx != 1:
+        params['prng_*'] = trial_idx
     net = Network(params, n_jobs=1)
 
     # global variables, should be node-independent
@@ -114,7 +117,7 @@ def simulate_dipole(net, n_trials=1, n_jobs=1):
         The dipole object or list of dipole objects if n_trials > 1
     """
     parallel, myfunc = _parallel_func(_simulate_single_trial, n_jobs=n_jobs)
-    dpl = parallel(myfunc(net.params) for _ in range(n_trials))
+    dpl = parallel(myfunc(net.params, idx) for idx in range(n_trials))
 
     return dpl