diff --git a/nipype/interfaces/spm/model.py b/nipype/interfaces/spm/model.py
index 89f2be33ea..a2e4a195ef 100644
--- a/nipype/interfaces/spm/model.py
+++ b/nipype/interfaces/spm/model.py
@@ -12,7 +12,7 @@
 
 # Local imports
 from ... import logging
-from ...utils.filemanip import ensure_list, simplify_list, split_filename
+from ...utils.filemanip import ensure_list, simplify_list, split_filename, load_spm_mat
 from ..base import (
     Bunch,
     traits,
@@ -312,12 +312,10 @@ def _parse_inputs(self):
         return einputs
 
     def _list_outputs(self):
-        import scipy.io as sio
-
         outputs = self._outputs().get()
         pth = os.path.dirname(self.inputs.spm_mat_file)
         outtype = "nii" if "12" in self.version.split(".")[0] else "img"
-        spm = sio.loadmat(self.inputs.spm_mat_file, struct_as_record=False)
+        spm = load_spm_mat(self.inputs.spm_mat_file, struct_as_record=False)
 
         betas = [vbeta.fname[0] for vbeta in spm["SPM"][0, 0].Vbeta[0]]
         if (
@@ -502,6 +500,10 @@ def _make_matlab_command(self, _):
 load(jobs{1}.stats{1}.con.spmmat{:});
 SPM.swd = '%s';
 save(jobs{1}.stats{1}.con.spmmat{:},'SPM');
+[msg,id] = lastwarn('');
+if strcmp(id,'MATLAB:save:sizeTooBigForMATFile')
+   save(jobs{1}.stats{1}.con.spmmat{:},'SPM','-v7.3');
+end
 names = SPM.xX.name;"""
             % (self.inputs.spm_mat_file, os.getcwd())
         ]
@@ -580,11 +582,9 @@ def _make_matlab_command(self, _):
         return "\n".join(script)
 
     def _list_outputs(self):
-        import scipy.io as sio
-
         outputs = self._outputs().get()
         pth, _ = os.path.split(self.inputs.spm_mat_file)
-        spm = sio.loadmat(self.inputs.spm_mat_file, struct_as_record=False)
+        spm = load_spm_mat(self.inputs.spm_mat_file, struct_as_record=False)
         con_images = []
         spmT_images = []
         for con in spm["SPM"][0, 0].xCon[0]:
diff --git a/nipype/utils/filemanip.py b/nipype/utils/filemanip.py
index ef52d00cc6..7e47d1519c 100644
--- a/nipype/utils/filemanip.py
+++ b/nipype/utils/filemanip.py
@@ -19,6 +19,7 @@
 from pathlib import Path
 import simplejson as json
 from time import sleep, time
+import scipy.io as sio
 
 from .. import logging, config, __version__ as version
 from .misc import is_container
@@ -932,3 +933,51 @@ def indirectory(path):
         yield
     finally:
         os.chdir(cwd)
+
+def load_spm_mat(spm_mat_file, **kwargs):
+    try:
+        mat = sio.loadmat(spm_mat_file, **kwargs)
+    except NotImplementedError:
+        import h5py
+        import numpy as np
+        mat = dict(SPM=np.array([[sio.matlab.mat_struct()]]))
+
+        # Get Vbeta, Vcon, and Vspm file names
+        with h5py.File(spm_mat_file, "r") as h5file:
+            fnames = dict()
+            try:
+                fnames["Vbeta"] = [u"".join(chr(c[0]) for c in h5file[obj_ref[0]]) for obj_ref in h5file["SPM"]["Vbeta"]["fname"]]
+            except Exception:
+                fnames["Vbeta"] = []
+            for contr_type in ["Vcon", "Vspm"]:
+                try:
+                    fnames[contr_type] = [u"".join(chr(c[0]) for c in h5file[obj_ref[0]]["fname"]) for obj_ref in h5file["SPM"]["xCon"][contr_type]]
+                except Exception:
+                    fnames[contr_type] = []
+
+        # Structure Vbeta as returned by scipy.io.loadmat
+        obj_list = []
+        for i in range(len(fnames["Vbeta"])):
+            obj = sio.matlab.mat_struct()
+            setattr(obj, "fname", np.array([fnames["Vbeta"][i]]))
+            obj_list.append(obj)
+        if len(obj_list) > 0:
+            setattr(mat["SPM"][0, 0], "Vbeta", np.array([obj_list]))
+        else:
+            setattr(mat["SPM"][0, 0], "Vbeta", np.empty((0, 0), dtype=object))
+
+        # Structure Vcon and Vspm as returned by scipy.io.loadmat
+        obj_list = []
+        for i in range(len(fnames["Vcon"])):
+            obj = sio.matlab.mat_struct()
+            for contr_type in ["Vcon", "Vspm"]:
+                temp = sio.matlab.mat_struct()
+                setattr(temp, "fname", np.array([fnames[contr_type][i]]))
+                setattr(obj, contr_type, np.array([[temp]]))
+            obj_list.append(obj)
+        if len(obj_list) > 0:
+            setattr(mat["SPM"][0, 0], "xCon", np.array([obj_list]))
+        else:
+            setattr(mat["SPM"][0, 0], "xCon", np.empty((0, 0), dtype=object))
+
+    return mat
\ No newline at end of file