diff --git a/niworkflows/anat/ants.py b/niworkflows/anat/ants.py
index 3f2c04940e6..6b460396047 100644
--- a/niworkflows/anat/ants.py
+++ b/niworkflows/anat/ants.py
@@ -507,6 +507,8 @@ def init_atropos_wf(
     padding=10,
     in_segmentation_model=tuple(ATROPOS_MODELS["T1w"].values()),
     bspline_fitting_distance=200,
+    adaptive_bspline_grid=False,
+    n4_iter=5,
     wm_prior=False,
 ):
     """
@@ -556,6 +558,12 @@ def init_atropos_wf(
         ``(4,4,2,3)`` uses K=4, CSF=4, GM=2, WM=3.
     bspline_fitting_distance : float
         The size of the b-spline mesh grid elements, in mm (default: 200)
+    adaptive_bspline_grid : :obj:`bool`
+        If true, defines the number of B-Spline mesh grid elements in each dimension rather
+        than using the isotropic distance given in ``bspline_fitting_distance``.
+    n4_iter : :obj:`int`
+        The number of B-Spline fitting iterations (default: 5). Fewer (e.g. 4) are recommended
+        for rodents and other non-human/non-adult cases.
     wm_prior : :obj:`bool`
         Whether the WM posterior obtained with ATROPOS should be regularized with a prior
         map (typically, mapped from the template). When ``wm_prior`` is ``True`` the input
@@ -755,7 +763,7 @@ def init_atropos_wf(
             dimension=3,
             save_bias=True,
             copy_header=True,
-            n_iterations=[50] * 5,
+            n_iterations=[50]*n4_iter,
             convergence_threshold=1e-7,
             shrink_factor=4,
             bspline_fitting_distance=bspline_fitting_distance,
@@ -875,6 +883,19 @@ def _argmax(in_dice):
             (apply_wm_prior, inu_n4_final, [("out", "weight_image")]),
         ])
         # fmt: on
+
+    if adaptive_bspline_grid:
+        # set INU bspline grid based on image shape
+        from ..utils.images import _bspline_grid
+        bspline_grid = pe.Node(niu.Function(function=_bspline_grid), name="bspline_grid")
+
+        # fmt:off
+        wf.connect([
+            (inputnode, bspline_grid, [(("in_files", _pop), "in_file")]),
+            (bspline_grid, inu_n4_final, [("out", "args")])
+        ])
+        # fmt:on
+
     return wf
 
 
@@ -1045,7 +1066,6 @@ def init_n4_only_wf(
             ]),
         ])
         # fmt: on
-
     return wf
 
 
diff --git a/niworkflows/utils/images.py b/niworkflows/utils/images.py
index c75a238b1b8..7a5ebf21bd6 100644
--- a/niworkflows/utils/images.py
+++ b/niworkflows/utils/images.py
@@ -324,3 +324,25 @@ def nii_ones_like(in_file, value, dtype, newpath=None):
     nii.to_filename(out_file)
 
     return out_file
+
+
+def _bspline_grid(in_file):
+    """
+    Estimate B-Spline fitting distance grid using the number of slices of ``in_file``.
+
+    Using slice number to determine grid sparsity is inspired by the conversation found at
+    https://itk.org/pipermail/community/2014-February/005036.html
+    
+    """
+    import nibabel as nb
+    import numpy as np
+    import math
+
+    img = nb.load(in_file)
+    zooms = img.header.get_zooms()[:3]
+    # find extent of each dimension wrt voxel sizes
+    extent = np.array(img.shape[:3]) * zooms
+    # convert ratio to integers
+    retval = [f"{math.ceil(i / extent[np.argmin(extent)])}" for i in extent]
+    # return string for command line call
+    return f"-b [{'x'.join(retval)}]"
diff --git a/niworkflows/workflows/epi/refmap.py b/niworkflows/workflows/epi/refmap.py
index ae2991467c7..915c88b3901 100644
--- a/niworkflows/workflows/epi/refmap.py
+++ b/niworkflows/workflows/epi/refmap.py
@@ -34,6 +34,8 @@
 def init_epi_reference_wf(
     omp_nthreads,
     auto_bold_nss=False,
+    adaptive_bspline_grid=False,
+    n4_iter=5,
     name="epi_reference_wf",
 ):
     """
@@ -84,6 +86,12 @@ def init_epi_reference_wf(
         If ``True``, determines nonsteady states in the beginning of the timeseries
         and selects them for the averaging of each run.
         IMPORTANT: this option applies only to BOLD EPIs.
+    adaptive_bspline_grid : :obj:`bool`
+        If ``True``, determines the number of B-Spline grid elements from data shape
+        and feeds them into N4BiasFieldCorrection, rather than setting an isotropic distance.
+    n4_iter : :obj:`int`
+        The number of B-Spline fitting iterations (default: 5). Fewer (e.g. 4) are recommended
+        for rodents and other non-human/non-adult cases.
 
     Inputs
     ------
@@ -158,7 +166,7 @@ def init_epi_reference_wf(
         N4BiasFieldCorrection(
             dimension=3,
             copy_header=True,
-            n_iterations=[50] * 5,
+            n_iterations=[50]*n4_iter,
             convergence_threshold=1e-7,
             shrink_factor=4,
         ),
@@ -166,6 +174,7 @@ def init_epi_reference_wf(
         name="n4_avgs",
         iterfield=["input_image"],
     )
+
     clip_bg_noise = pe.MapNode(
         IntensityClip(p_min=2.0, p_max=100.0),
         name="clip_bg_noise",
@@ -225,6 +234,20 @@ def _set_threads(in_list, maximum):
     else:
         wf.connect(inputnode, "t_masks", per_run_avgs, "t_mask")
 
+    # rodent-specific N4 settings
+    if adaptive_bspline_grid:
+        from ...utils.images import _bspline_grid
+        from ...utils.connections import pop_file as _pop
+        # set INU bspline grid based on voxel size
+        bspline_grid = pe.Node(niu.Function(function=_bspline_grid), name="bspline_grid")
+
+        # fmt:off
+        wf.connect([
+            (clip_avgs, bspline_grid, [(("out_file", _pop), "in_file")]),
+            (bspline_grid, n4_avgs, [("out", "args")])
+        ])
+        # fmt:on
+
     return wf
 
 
diff --git a/setup.cfg b/setup.cfg
index 07a6c149ee5..26a6e3e1dd4 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -30,7 +30,7 @@ install_requires =
     matplotlib >= 2.2.0
     nibabel >= 3.0.1
     nilearn >= 0.2.6, != 0.5.0, != 0.5.1
-    nipype >= 1.5.1
+    nipype @ git+https://github.com/nipy/nipype.git@master
     nitransforms >= 20.0.0rc3,<20.2
     numpy
     packaging