diff --git a/src/pySODM/models/utils.py b/src/pySODM/models/utils.py
index 8333790..394031c 100644
--- a/src/pySODM/models/utils.py
+++ b/src/pySODM/models/utils.py
@@ -5,10 +5,9 @@ def int_to_date(actual_start_date, t):
     date = actual_start_date + timedelta(days=t)
     return date
 
-def list_to_dict(y, shape_dictionary, retain_floats=True):
+def list_to_dict(y, shape_dictionary, retain_floats=False):
     """
     A function to reconstruct model states, given a flat array of values and a dictionary with the variables name and desired shapes
-    Retains floats.
 
     Parameters
     ----------
@@ -19,9 +18,6 @@ def list_to_dict(y, shape_dictionary, retain_floats=True):
     shape_dictionary: dict
         A dictionary containing the desired names and shapes of the output dictionary
     
-    retain_floats: bool
-        If False, float/np.float are converted in a 0D np.ndarray. Default: True.
-
     Returns
     -------
 
@@ -33,10 +29,11 @@ def list_to_dict(y, shape_dictionary, retain_floats=True):
     offset=0
     for s in shape_dictionary.values():
         n = np.prod(s)
-        # Reshape changes type of floats to np.ndarray which is not desirable
-        if ((n == 1) & (retain_floats==True)):
+        if ((n == 1) & (retain_floats == True)):
             restoredArray.append(y[offset])
+        elif ((n == 1) & (retain_floats == False)):
+            restoredArray.append(y[offset].reshape(s))
         else:
             restoredArray.append(y[offset:(offset+n)].reshape(s))
         offset+=n
-    return dict(zip(shape_dictionary.keys(), restoredArray))
\ No newline at end of file
+    return dict(zip(shape_dictionary.keys(), restoredArray))
diff --git a/tutorials/IDD/spatial_SIR/models.py b/tutorials/IDD/spatial_SIR/models.py
index cb5d93d..180ce74 100644
--- a/tutorials/IDD/spatial_SIR/models.py
+++ b/tutorials/IDD/spatial_SIR/models.py
@@ -61,7 +61,7 @@ def integrate(t, S, I, R, beta, gamma, f_v, N, M):
         dI_v = beta * S_v * np.transpose(matmul_2D_3D_matrix(np.transpose(I_v/T_v), f_v*N))
 
         # distribute the number of new infections on visited patch to the home patch 
-        dI_v = S * np.transpose(M @ np.transpose(dI_v/S_v))
+        dI_v = S * np.transpose(np.atleast_2d(M) @ np.transpose(dI_v/S_v))
 
         # Calculate differentials
         dS = - (dI_h + dI_v)
diff --git a/tutorials/IDD/spatial_SIR/simulation.py b/tutorials/IDD/spatial_SIR/simulation.py
index e7d4615..f54fb8e 100644
--- a/tutorials/IDD/spatial_SIR/simulation.py
+++ b/tutorials/IDD/spatial_SIR/simulation.py
@@ -35,17 +35,18 @@
 
 # OR
 # Exactly the same model without age groups
+# Validated the reduction to one spatial unit as well (becomes a plain SIR)
 
 coordinates = {'age': ['0+'],  
                'location': ['A', 'B']
                 }
-init_states = {'S': np.array([[500, 5000]]),    
-               'I': np.array([[1, 0]])
+init_states = {'S': np.atleast_2d([500, 5000]),    
+               'I': np.atleast_2d([1, 0])
                }
-params = {'beta': 0.03,                                # infectivity (-)
+params = {'beta': 0.03,                                 # infectivity (-)
           'gamma': 5,                                   # duration of infection (d)
           'f_v': 0.1,                                   # fraction of total contacts on visited patch
-          'N': np.mean(np.sum(params['N'], axis=1)),    # contact matrix (BE, Van Hoang, 2020)
+          'N': np.mean(np.sum(params['N'], axis=1)),    # contact matrix
           'M': np.array([[0.8, 0.2], [0, 1]])           # origin-destination mobility matrix
           }