diff --git a/niworkflows/conftest.py b/niworkflows/conftest.py
index 0408bc8155e..75a19fe26db 100644
--- a/niworkflows/conftest.py
+++ b/niworkflows/conftest.py
@@ -48,6 +48,14 @@ def find_resource_or_skip(resource):
     return pathlike
 
 
+@pytest.fixture(scope="session", autouse=True)
+def legacy_printoptions():
+    from packaging.version import Version
+
+    if Version(np.__version__) >= Version("1.22"):
+        np.set_printoptions(legacy="1.21")
+
+
 @pytest.fixture(autouse=True)
 def add_np(doctest_namespace):
     from .utils.bids import collect_data
diff --git a/niworkflows/interfaces/images.py b/niworkflows/interfaces/images.py
index 2669904196a..916c200c7ab 100644
--- a/niworkflows/interfaces/images.py
+++ b/niworkflows/interfaces/images.py
@@ -749,6 +749,6 @@ def _run_interface(self, runtime):
 
         output = np.vstack((self.inputs.class_labels, series.astype(str)))
         self._results["out_file"] = os.path.join(runtime.cwd, self.inputs.out_file)
-        np.savetxt(self._results["out_file"], output, fmt=b"%s", delimiter="\t")
+        np.savetxt(self._results["out_file"], output, fmt="%s", delimiter="\t")
 
         return runtime
diff --git a/niworkflows/tests/test_confounds.py b/niworkflows/tests/test_confounds.py
index 0cd871a417f..e09aac30756 100644
--- a/niworkflows/tests/test_confounds.py
+++ b/niworkflows/tests/test_confounds.py
@@ -98,15 +98,15 @@ def test_expansion_derivatives_and_powers():
         {
             "a": [-1, -2, -3, -4, -5],
             "a_power2": [1, 4, 9, 16, 25],
-            "a_derivative1": [np.NaN, -1, -1, -1, -1],
-            "a_derivative1_power2": [np.NaN, 1, 1, 1, 1],
-            "b_derivative1": [np.NaN, 0, 0, 0, 0],
-            "b_derivative1_power2": [np.NaN, 0, 0, 0, 0],
-            "c_power2_derivative1": [np.NaN, 1, -1, 1, -1],
-            "c_power2_derivative2": [np.NaN, np.NaN, -2, 2, -2],
+            "a_derivative1": [np.nan, -1, -1, -1, -1],
+            "a_derivative1_power2": [np.nan, 1, 1, 1, 1],
+            "b_derivative1": [np.nan, 0, 0, 0, 0],
+            "b_derivative1_power2": [np.nan, 0, 0, 0, 0],
+            "c_power2_derivative1": [np.nan, 1, -1, 1, -1],
+            "c_power2_derivative2": [np.nan, np.nan, -2, 2, -2],
             "d": [9, 7, 5, 3, 1],
             "e": [0, 0, 0, 0, 0],
-            "f": [np.NaN, 6, 4, 2, 0],
+            "f": [np.nan, 6, 4, 2, 0],
         }
     )
     exp_data = _expand_test(model_formula)
@@ -122,10 +122,10 @@ def test_expansion_na_robustness():
     model_formula = "(dd1(f))^^2"
     expected_data = pd.DataFrame(
         {
-            "f": [np.NaN, 6, 4, 2, 0],
-            "f_power2": [np.NaN, 36, 16, 4, 0],
-            "f_derivative1": [np.NaN, np.NaN, -2, -2, -2],
-            "f_derivative1_power2": [np.NaN, np.NaN, 4, 4, 4],
+            "f": [np.nan, 6, 4, 2, 0],
+            "f_power2": [np.nan, 36, 16, 4, 0],
+            "f_derivative1": [np.nan, np.nan, -2, -2, -2],
+            "f_derivative1_power2": [np.nan, np.nan, 4, 4, 4],
         }
     )
     exp_data = _expand_test(model_formula)