diff --git a/examples/preprocessing/muscle_ica.py b/examples/preprocessing/muscle_ica.py
index 81956d7acda..6a758f8a08f 100644
--- a/examples/preprocessing/muscle_ica.py
+++ b/examples/preprocessing/muscle_ica.py
@@ -92,10 +92,9 @@
 
 # %%
 # Compare to EEGBCI dataset
-# Note that the ICA for the third subject got
 
 muscle_idxs = {1: [9], 2: [5, 10, 12, 13, 14]}
-for sub in (1, 2, 3):
+for sub, muscle_idx in muscle_idxs.items():
     raw = mne.io.read_raw_edf(
         mne.datasets.eegbci.load_data(subject=sub, runs=(1,))[0], preload=True)
     mne.datasets.eegbci.standardize(raw)  # set channel names
diff --git a/mne/preprocessing/ica.py b/mne/preprocessing/ica.py
index 78b46b5a5a1..b0a659e9c71 100644
--- a/mne/preprocessing/ica.py
+++ b/mne/preprocessing/ica.py
@@ -1648,16 +1648,16 @@ def find_bads_muscle(self, inst, threshold=0.5, start=None,
         # compute metric #2: distance from the vertex of focus
         components_norm = abs(components) / np.max(abs(components), axis=0)
         pos = _find_topomap_coords(inst.info, picks='data', sphere=sphere)
-        assert pos.shape[0] == components.shape[0]  # pos for each component
+        assert pos.shape[0] == components.shape[0]  # pos for each sensor
         pos -= pos.mean(axis=0)  # center
         dists = np.linalg.norm(pos, axis=1)
         dists /= dists.max()
         focus_dists = np.dot(dists, components_norm)
 
         # compute metric #3: smoothness
+        smoothnesses = np.zeros((components.shape[1],))
         dists = squareform(pdist(pos))
         dists = 1 - (dists / dists.max())  # invert
-        smoothnesses = np.zeros((components.shape[1],))
         for idx, comp in enumerate(components.T):
             comp_dists = squareform(pdist(comp[:, np.newaxis]))
             comp_dists /= comp_dists.max()
diff --git a/mne/preprocessing/tests/test_ica.py b/mne/preprocessing/tests/test_ica.py
index d0df00b701f..a56543d9e74 100644
--- a/mne/preprocessing/tests/test_ica.py
+++ b/mne/preprocessing/tests/test_ica.py
@@ -1324,8 +1324,8 @@ def test_ica_labels():
 
     scores = ica.find_bads_muscle(raw)[1]
     assert 'muscle' in ica.labels_
-    assert ica.labels_['muscle'] == [0, 1, 2, 3]
-    assert_allclose(scores, [0.2, 0.25, 0.35, 0.55], atol=0.03)
+    assert ica.labels_['muscle'] == [0]
+    assert_allclose(scores, [0.56, 0.01, 0.03, 0.00], atol=0.03)
 
 
 @requires_sklearn
diff --git a/tutorials/preprocessing/40_artifact_correction_ica.py b/tutorials/preprocessing/40_artifact_correction_ica.py
index 729c01510e2..157307c7628 100644
--- a/tutorials/preprocessing/40_artifact_correction_ica.py
+++ b/tutorials/preprocessing/40_artifact_correction_ica.py
@@ -459,6 +459,11 @@
 # of automated approaches like `~mne.preprocessing.ICA.find_bads_ecg`
 # before accepting them.
 
+# %%
+# For EEG, activation of muscles for postural control of the head and neck
+# contaminate the signal as well. This is usually not detected by MEG. For
+# an example showing how to remove these components, see :ref:`ex-muscle-ica`.
+
 # clean up memory before moving on
 del raw, ica, new_ica