diff --git a/src/pypulseq/Sequence/block.py b/src/pypulseq/Sequence/block.py
index 9b9d0da2..9077d0d5 100644
--- a/src/pypulseq/Sequence/block.py
+++ b/src/pypulseq/Sequence/block.py
@@ -226,7 +226,7 @@ def set_block(self, block_index: int, *args: SimpleNamespace) -> None:
                     if prev_type == 't':
                         last = 0
                     elif prev_type == 'g':
-                        last = prev_lib['data'][5]
+                        last = prev_lib['data'][2]  # v150
 
             # Check whether the difference between the last gradient value and
             # the first value of the new gradient is achievable with the
@@ -245,7 +245,7 @@ def set_block(self, block_index: int, *args: SimpleNamespace) -> None:
                     if next_type == 't':
                         first = 0
                     elif next_type == 'g':
-                        first = next_lib['data'][4]
+                        first = next_lib['data'][1]  # v150
                 else:
                     first = 0
 
@@ -260,7 +260,7 @@ def set_block(self, block_index: int, *args: SimpleNamespace) -> None:
             raise RuntimeError('First gradient in the the first block has to start at 0.')
 
         if (
-            grad_to_check.stop[1] > self.system.max_slew * self.system.grad_raster_time
+            abs(grad_to_check.stop[1]) > self.system.max_slew * self.system.grad_raster_time
             and abs(grad_to_check.stop[0] - duration) > 1e-7
         ):
             raise RuntimeError("A gradient that doesn't end at zero needs to be aligned to the block boundary.")
@@ -331,9 +331,9 @@ def get_block(self, block_index: int) -> SimpleNamespace:
             grad.channel = grad_channels[i][1]
             if grad.type == 'grad':
                 amplitude = lib_data[0]
-                shape_id = lib_data[1]
-                time_id = lib_data[2]
-                delay = lib_data[3]
+                shape_id = lib_data[3]  # change in v150
+                time_id = lib_data[4]  # change in v150
+                delay = lib_data[5]  # change in v150
                 shape_data = self.shape_library.data[shape_id]
                 compressed.num_samples = shape_data[0]
                 compressed.data = shape_data[1:]
@@ -343,22 +343,36 @@ def get_block(self, block_index: int) -> SimpleNamespace:
                 if time_id == 0:
                     grad.tt = (np.arange(1, len(g) + 1) - 0.5) * self.grad_raster_time
                     t_end = len(g) * self.grad_raster_time
+                    grad.area = sum(grad.waveform) * self.grad_raster_time
+                elif time_id == -1:
+                    # Gradient with oversampling by a factor of 2
+                    grad.tt = 0.5 * (np.arange(1, len(g) + 1)) * self.grad_raster_time
+                    if len(grad.tt) != len(grad.waveform):
+                        raise ValueError(
+                            f'Mismatch between time shape length ({len(grad.tt)}) and gradient shape length ({len(grad.waveform)}).'
+                        )
+                    if len(grad.waveform) % 2 != 1:
+                        raise ValueError('Oversampled gradient waveforms must have odd number of samples')
+                    t_end = (len(g) + 1) * self.grad_raster_time
+                    grad.area = sum(grad.waveform[::2]) * self.grad_raster_time  # remove oversampling
                 else:
                     t_shape_data = self.shape_library.data[time_id]
                     compressed.num_samples = t_shape_data[0]
                     compressed.data = t_shape_data[1:]
                     grad.tt = decompress_shape(compressed) * self.grad_raster_time
-
-                    assert len(grad.waveform) == len(grad.tt)
+                    if len(grad.tt) != len(grad.waveform):
+                        raise ValueError(
+                            f'Mismatch between time shape length ({len(grad.tt)}) and gradient shape length ({len(grad.waveform)}).'
+                        )
                     t_end = grad.tt[-1]
+                    grad.area = 0.5 * sum((grad.tt[1:] - grad.tt[:-1]) * (grad.waveform[1:] + grad.waveform[:-1]))
 
                 grad.shape_id = shape_id
                 grad.time_id = time_id
                 grad.delay = delay
                 grad.shape_dur = t_end
-                if len(lib_data) > 5:
-                    grad.first = lib_data[4]
-                    grad.last = lib_data[5]
+                grad.first = lib_data[1]  # change in v150 - we always have first/last now
+                grad.last = lib_data[2]  # change in v150 - we always have first/last now
             else:
                 grad.amplitude = lib_data[0]
                 grad.rise_time = lib_data[1]
@@ -368,6 +382,8 @@ def get_block(self, block_index: int) -> SimpleNamespace:
                 grad.area = grad.amplitude * (grad.flat_time + grad.rise_time / 2 + grad.fall_time / 2)
                 grad.flat_area = grad.amplitude * grad.flat_time
 
+            # TODO: add optional grad ID from raw_block
+
             setattr(block, grad_channels[i], grad)
 
     # ADC
@@ -572,37 +588,45 @@ def register_grad_event(self, event: SimpleNamespace) -> Union[int, Tuple[int, L
     """
     may_exist = True
     any_changed = False
+
     if event.type == 'grad':
-        amplitude = np.abs(event.waveform).max()
+        amplitude = np.max(np.abs(event.waveform))
         if amplitude > 0:
             fnz = event.waveform[np.nonzero(event.waveform)[0][0]]
-            amplitude *= np.sign(fnz) if fnz != 0 else 1  # Workaround for np.sign(0) = 0
+            amplitude *= np.sign(fnz) if fnz != 0 else 1
 
+        # Shape ID initialization
         if hasattr(event, 'shape_IDs'):
             shape_IDs = event.shape_IDs
         else:
             shape_IDs = [0, 0]
-            if amplitude != 0:
-                g = event.waveform / amplitude
-            else:
-                g = event.waveform
+
+            # Shape for waveform
+            g = event.waveform / amplitude if amplitude != 0 else event.waveform
             c_shape = compress_shape(g)
             s_data = np.concatenate(([c_shape.num_samples], c_shape.data))
             shape_IDs[0], found = self.shape_library.find_or_insert(s_data)
-            may_exist = may_exist & found
+            may_exist = may_exist and found
             any_changed = any_changed or found
 
-            # Check whether tt == np.arange(len(event.tt)) * self.grad_raster_time + 0.5
-            tt_regular = (np.floor(event.tt / self.grad_raster_time) == np.arange(len(event.tt))).all()
+            # Shape for timing
+            c_time = compress_shape(event.tt / self.grad_raster_time)
+            t_data = np.concatenate(([c_time.num_samples], c_time.data))
 
-            if not tt_regular:
-                c_time = compress_shape(event.tt / self.grad_raster_time)
-                t_data = np.concatenate(([c_time.num_samples], c_time.data))
+            if len(c_time.data) == 4 and np.allclose(c_time.data, [0.5, 1, 1, c_time.num_samples - 3]):
+                # Standard raster → leave shape_IDs[1] as 0
+                pass
+            elif len(c_time.data) == 3 and np.allclose(c_time.data, [0.5, 0.5, c_time.num_samples - 2]):
+                # Half-raster → set to -1 as special flag
+                shape_IDs[1] = -1
+            else:
                 shape_IDs[1], found = self.shape_library.find_or_insert(t_data)
-                may_exist = may_exist & found
+                may_exist = may_exist and found
                 any_changed = any_changed or found
 
-        data = (amplitude, *shape_IDs, event.delay, event.first, event.last)
+        # Updated data layout to match MATLAB v1.5.0 ordering
+        data = (amplitude, event.first, event.last, *shape_IDs, event.delay)
+
     elif event.type == 'trap':
         data = (
             event.amplitude,
@@ -625,6 +649,9 @@ def register_grad_event(self, event: SimpleNamespace) -> Union[int, Tuple[int, L
     if self.use_block_cache and any_changed:
         self.block_cache.clear()
 
+    if hasattr(event, 'name'):
+        self.grad_id_to_name_map[grad_id] = event.name
+
     if event.type == 'grad':
         return grad_id, shape_IDs
     elif event.type == 'trap':
diff --git a/src/pypulseq/Sequence/read_seq.py b/src/pypulseq/Sequence/read_seq.py
index 6a160289..ec96a9ad 100644
--- a/src/pypulseq/Sequence/read_seq.py
+++ b/src/pypulseq/Sequence/read_seq.py
@@ -152,7 +152,9 @@ def read(self, path: str, detect_rf_use: bool = False, remove_duplicates: bool =
                 else:  # 1.3.x and below
                     self.rf_library = __read_events(input_file, (1, 1, 1, 1e-6, 1, 1), event_library=self.rf_library)
         elif section == '[GRADIENTS]':
-            if version_combined >= 1004000:  # 1.4.x format
+            if version_combined >= 1005000:  # 1.5.x format
+                self.grad_library = __read_events(input_file, (1, 1, 1, 1, 1, 1e-6), 'g', self.grad_library)
+            elif version_combined >= 1004000:  # 1.4.x format
                 self.grad_library = __read_events(input_file, (1, 1, 1, 1e-6), 'g', self.grad_library)
             else:  # 1.3.x and below
                 self.grad_library = __read_events(input_file, (1, 1, 1e-6), 'g', self.grad_library)
diff --git a/src/pypulseq/Sequence/sequence.py b/src/pypulseq/Sequence/sequence.py
index 259b5cd6..614088d3 100644
--- a/src/pypulseq/Sequence/sequence.py
+++ b/src/pypulseq/Sequence/sequence.py
@@ -94,6 +94,7 @@ def __init__(self, system: Union[Opts, None] = None, use_block_cache: bool = Tru
         self.signature_value = ''
         self.rf_id_to_name_map = {}
         self.adc_id_to_name_map = {}
+        self.grad_id_to_name_map = {}
 
         self.block_durations = {}
         self.extension_numeric_idx = []
@@ -1197,7 +1198,7 @@ def remove_duplicates(self, in_place: bool = False) -> Self:
         for grad_id in seq_copy.grad_library.data:
             if seq_copy.grad_library.type[grad_id] == 'g':
                 data = seq_copy.grad_library.data[grad_id]
-                new_data = (data[0],) + (mapping[data[1]], mapping[data[2]]) + data[3:]
+                new_data = data[0:3] + (mapping[data[3]], mapping[data[4]]) + (data[5],)
                 if data != new_data:
                     seq_copy.grad_library.update(grad_id, None, new_data)
 
diff --git a/src/pypulseq/Sequence/write_seq.py b/src/pypulseq/Sequence/write_seq.py
index 48f661c2..ba3584a4 100644
--- a/src/pypulseq/Sequence/write_seq.py
+++ b/src/pypulseq/Sequence/write_seq.py
@@ -127,16 +127,16 @@ def write(self, file_name: Union[str, Path], create_signature, remove_duplicates
             output_file.write(
                 '#   time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster)\n'
             )
-            output_file.write('# id amplitude amp_shape_id time_shape_id delay\n')
-            output_file.write('# ..      Hz/m       ..         ..          us\n')
+            output_file.write('# id amplitude first last amp_shape_id time_shape_id delay\n')
+            output_file.write('# ..      Hz/m  Hz/m Hz/m        ..         ..          us\n')
             output_file.write('[GRADIENTS]\n')
-            id_format_str = '{:.0f} {:12g} {:.0f} {:.0f} {:.0f}\n'  # Refer lines 20-21
+            id_format_str = '{:.0f} {:12g} {:12g} {:12g} {:.0f} {:.0f} {:.0f}\n'  # Refer lines 20-21
             keys = np.array(list(self.grad_library.data.keys()))
             for k in keys[arb_grad_mask]:
                 s = id_format_str.format(
                     k,
-                    *self.grad_library.data[k][:3],
-                    round(self.grad_library.data[k][3] * 1e6),
+                    *self.grad_library.data[k][:5],
+                    round(self.grad_library.data[k][5] * 1e6),
                 )
                 output_file.write(s)
             output_file.write('\n')
diff --git a/src/pypulseq/add_gradients.py b/src/pypulseq/add_gradients.py
index c4e2aacd..13bfef07 100644
--- a/src/pypulseq/add_gradients.py
+++ b/src/pypulseq/add_gradients.py
@@ -86,24 +86,29 @@ def add_gradients(
         return grad
 
     # Find out the general delay of all gradients and other statistics
-    delays, firsts, lasts, durs, is_trap, is_arb = [], [], [], [], [], []
+    delays, firsts, lasts, durs, is_trap, is_arb, is_osa = [], [], [], [], [], [], []
     for ii in range(len(grads)):
         if grads[ii].channel != channel:
             raise ValueError('Cannot add gradients on different channels.')
 
         delays.append(grads[ii].delay)
-        firsts.append(grads[ii].first)
-        lasts.append(grads[ii].last)
         durs.append(calc_duration(grads[ii]))
         is_trap.append(grads[ii].type == 'trap')
         if is_trap[-1]:
             is_arb.append(False)
+            is_osa.append(False)
+            firsts.append(0.0)
+            lasts.append(0.0)
         else:
-            tt_rast = grads[ii].tt / system.grad_raster_time - 0.5
-            is_arb.append(np.all(np.abs(tt_rast - np.arange(len(tt_rast)))) < eps)
+            tt_rast = grads[ii].tt / system.grad_raster_time
+            is_arb.append(np.all(np.abs(tt_rast + 0.5 - np.arange(1, len(tt_rast) + 1))) < eps)
+            is_osa.append(np.all(np.abs(tt_rast - 0.5 * np.arange(1, len(tt_rast) + 1)) < eps))
+            firsts.append(grads[ii].first)
+            lasts.append(grads[ii].last)
 
     # Check if we only have arbitrary grads on irregular time samplings, optionally mixed with trapezoids
-    if np.all(np.logical_or(is_trap, np.logical_not(is_arb))):
+    is_etrap = np.logical_and.reduce((np.logical_not(is_trap), np.logical_not(is_arb), np.logical_not(is_osa)))
+    if np.all(np.logical_or(is_trap, is_etrap)):
         # Keep shapes still rather simple
         times = []
         for ii in range(len(grads)):
@@ -161,21 +166,35 @@ def add_gradients(
     # Convert to numpy.ndarray for fancy-indexing later on
     firsts, lasts = np.array(firsts), np.array(lasts)
     common_delay = np.min(delays)
+    total_duration = np.max(durs)
     durs = np.array(durs)
 
     # Convert everything to a regularly-sampled waveform
     waveforms = {}
     max_length = 0
+
+    if np.any(is_osa):
+        target_raster = 0.5 * system.grad_raster_time
+    else:
+        target_raster = system.grad_raster_time
+
     for ii in range(len(grads)):
         g = grads[ii]
         if g.type == 'grad':
-            if is_arb[ii]:
-                waveforms[ii] = g.waveform
+            if is_arb[ii] or is_osa[ii]:
+                if np.any(is_osa) and is_arb[ii]:  # Porting MATLAB here, maybe a bit ugly
+                    # Interpolate missing samples
+                    idx = np.arange(0, len(g.waveform) - 0.5 + eps, 0.5)
+                    wf = g.waveform
+                    interp_waveform = 0.5 * (wf[np.floor(idx).astype(int)] + wf[np.ceil(idx).astype(int)])
+                    waveforms[ii] = interp_waveform
+                else:
+                    waveforms[ii] = g.waveform
             else:
                 waveforms[ii] = points_to_waveform(
                     amplitudes=g.waveform,
                     times=g.tt,
-                    grad_raster_time=system.grad_raster_time,
+                    grad_raster_time=target_raster,
                 )
         elif g.type == 'trap':
             if g.flat_time > 0:  # Triangle or trapezoid
@@ -200,7 +219,7 @@ def add_gradients(
             waveforms[ii] = points_to_waveform(
                 amplitudes=amplitudes,
                 times=times,
-                grad_raster_time=system.grad_raster_time,
+                grad_raster_time=target_raster,
             )
         else:
             raise ValueError('Unknown gradient type')
@@ -228,6 +247,9 @@ def add_gradients(
         max_slew=max_slew,
         max_grad=max_grad,
         delay=common_delay,
+        oversampling=np.any(is_osa),
+        first=np.sum(firsts[delays == common_delay]),
+        last=np.sum(lasts[durs == total_duration]),
     )
     # Fix the first and the last values
     # First is defined by the sum of firsts with the minimal delay (common_delay)
diff --git a/src/pypulseq/add_ramps.py b/src/pypulseq/add_ramps.py
index c1c3c8dc..0bffa6b4 100644
--- a/src/pypulseq/add_ramps.py
+++ b/src/pypulseq/add_ramps.py
@@ -14,6 +14,7 @@ def add_ramps(
     max_slew: int = 0,
     rf: Union[SimpleNamespace, None] = None,
     system: Union[Opts, None] = None,
+    oversampling: bool = False,
 ) -> List[np.ndarray]:
     """
     Add segments to the trajectory to ramp to and from the given trajectory.
@@ -32,6 +33,8 @@ def add_ramps(
         Maximum gradient amplitude.
     max_slew : int, default=0
         Maximum slew rate.
+    oversampling : bool, default=False
+        Boolean flag to indicate if gradient is oversampled by a factor of 2.
 
     Returns
     -------
@@ -64,8 +67,8 @@ def add_ramps(
     num_channels = k.shape[0]
     k = np.vstack((k, np.zeros((3 - num_channels, k.shape[1]))))  # Pad with zeros if needed
 
-    k_up, ok1 = calc_ramp(k0=np.zeros((3, 2)), k_end=k[:, :2], system=system)
-    k_down, ok2 = calc_ramp(k0=k[:, -2:], k_end=np.zeros((3, 2)), system=system)
+    k_up, ok1 = calc_ramp(k0=np.zeros((3, 2)), k_end=k[:, :2], system=system, oversampling=oversampling)
+    k_down, ok2 = calc_ramp(k0=k[:, -2:], k_end=np.zeros((3, 2)), system=system, oversampling=oversampling)
     if not (ok1 and ok2):
         raise RuntimeError('Failed to calculate gradient ramps')
 
diff --git a/src/pypulseq/calc_ramp.py b/src/pypulseq/calc_ramp.py
index fd3342e1..0b3bbd3b 100644
--- a/src/pypulseq/calc_ramp.py
+++ b/src/pypulseq/calc_ramp.py
@@ -12,6 +12,7 @@ def calc_ramp(
     max_points: int = 500,
     max_slew: Union[np.ndarray, None] = None,
     system: Union[Opts, None] = None,
+    oversampling: bool = False,
 ) -> Tuple[np.ndarray, bool]:
     """
     Join the points `k0` and `k_end` in three-dimensional  k-space in minimal time, observing the gradient and slew
@@ -34,6 +35,8 @@ def calc_ramp(
         Maximum total slew rate. Either a single value or one value for each coordinate, of shape `[3, 1]`.
     system : Opts, default=Opts()
         System limits.
+    oversampling : bool, default=False
+        Boolean flag to indicate if gradient is oversampled by a factor of 2.
 
     Returns
     -------
@@ -321,7 +324,10 @@ def __joinright1(k0, k_end, use_points, G0, G_end):
     if np.all(np.where(max_slew <= 0)):
         max_slew = [system.max_slew]
 
-    grad_raster = system.grad_raster_time
+    if oversampling:
+        grad_raster = 0.5 * system.grad_raster_time
+    else:
+        grad_raster = system.grad_raster_time
 
     if len(max_grad) == 1 and len(max_slew) == 1:
         mode = 0
diff --git a/src/pypulseq/make_arbitrary_grad.py b/src/pypulseq/make_arbitrary_grad.py
index 720e94c6..41793143 100644
--- a/src/pypulseq/make_arbitrary_grad.py
+++ b/src/pypulseq/make_arbitrary_grad.py
@@ -1,3 +1,4 @@
+# import warnings
 from types import SimpleNamespace
 from typing import Union
 
@@ -17,6 +18,7 @@ def make_arbitrary_grad(
     max_grad: Union[float, None] = None,
     max_slew: Union[float, None] = None,
     system: Union[Opts, None] = None,
+    oversampling: bool = False,
 ) -> SimpleNamespace:
     """
     Creates a gradient event from an arbitrary waveform.
@@ -51,6 +53,8 @@ def make_arbitrary_grad(
         Will default to `system.max_slew` if not provided.
     delay : float, default=0
         Delay in seconds (s).
+    oversampling : bool, default=False
+        Boolean flag to indicate if gradient is oversampled by a factor of 2.
 
     Returns
     -------
@@ -76,28 +80,57 @@ def make_arbitrary_grad(
     if channel not in ['x', 'y', 'z']:
         raise ValueError(f'Invalid channel. Must be one of x, y or z. Passed: {channel}')
 
-    slew_rate = np.diff(waveform) / system.grad_raster_time
+    if first is None:
+        # warnings.warn(
+        #     'it will be compulsory to provide the first point of the gradient shape in the future releases; finding the first by extrapolation for now...',
+        #     FutureWarning,
+        # )
+        if oversampling:
+            first = 2 * waveform[0] - waveform[1]  # extrapolate by 1 gradient raster
+        else:
+            first = (3 * waveform[0] - waveform[1]) * 0.5  # extrapolate by 1/2 gradient of the raster
+
+    if last is None:
+        # warnings.warn(
+        #     'it will be compulsory to provide the last point of the gradient shape in the future releases; finding the last by extrapolation for now...',
+        #     FutureWarning,
+        # )
+        if oversampling:
+            last = 2 * waveform[-1] - waveform[-2]  # extrapolate by 1 gradient raster
+        else:
+            last = (3 * waveform[-1] - waveform[-2]) * 0.5  # extrapolate by 1/2 gradient of the raster
+
+    # Slew rate calculation
+    if oversampling:
+        slew_rate = np.concatenate([[first - waveform[0]], np.diff(waveform), [last - waveform[-1]]]) / (
+            system.grad_raster_time * 2
+        )
+    else:
+        slew_rate = (
+            np.concatenate([[2 * (first - waveform[0])], np.diff(waveform), [2 * (waveform[-1] - last)]])
+            / system.grad_raster_time
+        )
+
     if max(abs(slew_rate)) > max_slew * (1 + eps):
         raise ValueError(f'Slew rate violation {max(abs(slew_rate)) / max_slew * 100}')
     if max(abs(waveform)) > max_grad + eps:
         raise ValueError(f'Gradient amplitude violation {max(abs(waveform)) / max_grad * 100}')
 
-    if first is None:
-        first = (3 * waveform[0] - waveform[1]) * 0.5  # linear extrapolation
-
-    if last is None:
-        last = (3 * waveform[-1] - waveform[-2]) * 0.5  # linear extrapolation
-
     grad = SimpleNamespace()
     grad.type = 'grad'
     grad.channel = channel
     grad.waveform = waveform
     grad.delay = delay
-    grad.tt = (np.arange(len(waveform)) + 0.5) * system.grad_raster_time
-    grad.shape_dur = len(waveform) * system.grad_raster_time
-    grad.first = first
-    grad.last = last
-    grad.area = (waveform * system.grad_raster_time).sum()
+    if oversampling:
+        if len(waveform) % 2 == 0:
+            raise ValueError('When oversampling is active, waveform must have an odd number of samples')
+        grad.area = (waveform[::2] * system.grad_raster_time).sum()
+        grad.tt = np.arange(1, len(waveform) + 1) * 0.5 * system.grad_raster_time
+        grad.shape_dur = (len(waveform) + 1) * 0.5 * system.grad_raster_time
+    else:
+        grad.area = (waveform * system.grad_raster_time).sum()
+        grad.tt = (np.arange(len(waveform)) + 0.5) * system.grad_raster_time
+        grad.shape_dur = len(waveform) * system.grad_raster_time
 
     if trace_enabled():
         grad.trace = trace()
diff --git a/src/pypulseq/scale_grad.py b/src/pypulseq/scale_grad.py
index 37af6e22..32458798 100644
--- a/src/pypulseq/scale_grad.py
+++ b/src/pypulseq/scale_grad.py
@@ -1,8 +1,14 @@
 from copy import copy
 from types import SimpleNamespace
+from typing import Union
 
+import numpy as np
 
-def scale_grad(grad: SimpleNamespace, scale: float) -> SimpleNamespace:
+from pypulseq import eps
+from pypulseq.opts import Opts
+
+
+def scale_grad(grad: SimpleNamespace, scale: float, system: Union[Opts, None] = None) -> SimpleNamespace:
     """
     Scales the gradient with the scalar.
 
@@ -12,6 +18,8 @@ def scale_grad(grad: SimpleNamespace, scale: float) -> SimpleNamespace:
         Gradient event to be scaled.
     scale : float
         Scaling factor.
+    system : Opts, default=Opts()
+        System limits.
 
     Returns
     -------
@@ -22,12 +30,36 @@ def scale_grad(grad: SimpleNamespace, scale: float) -> SimpleNamespace:
     scaled_grad = copy(grad)
     if scaled_grad.type == 'trap':
         scaled_grad.amplitude = scaled_grad.amplitude * scale
-        scaled_grad.area = scaled_grad.area * scale
         scaled_grad.flat_area = scaled_grad.flat_area * scale
+        if system is not None:
+            if abs(scaled_grad.amplitude) > system.max_grad:
+                raise ValueError(
+                    f'scale_grad: maximum amplitude exceeded {100 * abs(scaled_grad.amplitude) / system.max_grad} %'
+                )
+            if (
+                abs(grad.amplitude) > eps
+                and abs(scaled_grad.amplitude) / min(scaled_grad.rise_time, scaled_grad.fall_time) > system.max_slew
+            ):
+                raise ValueError(
+                    'mr.scale_grad: maximum slew rate exceeded {100 * abs(scaled_grad.amplitude) / min(scaled_grad.rise_time, scaled_grad.fall_time) / system.max_slew} %'
+                )
+
     else:
         scaled_grad.waveform = scaled_grad.waveform * scale
         scaled_grad.first = scaled_grad.first * scale
         scaled_grad.last = scaled_grad.last * scale
+        if system is not None:
+            if max(abs(scaled_grad.waveform)) > system.max_grad:
+                raise ValueError(
+                    f'scale_grad: maximum amplitude exceeded {100 * max(abs(scaled_grad.waveform)) / system.max_grad} %'
+                )
+            if max(abs(scaled_grad.waveform)) > eps:
+                scaled_grad_max_abs_slew = max(abs(np.diff(scaled_grad.waveform) / np.diff(grad.tt)))
+                if scaled_grad_max_abs_slew > system.max_slew:
+                    raise ValueError(
+                        f'scale_grad: maximum slew rate exceeded {100 * scaled_grad_max_abs_slew / system.max_slew} %'
+                    )
+    scaled_grad.area = scaled_grad.area * scale
 
     if hasattr(scaled_grad, 'id'):
         delattr(scaled_grad, 'id')
diff --git a/tests/expected_output/write_epi_se_rs.seq b/tests/expected_output/write_epi_se_rs.seq
index a82513dd..ebcdd783 100644
--- a/tests/expected_output/write_epi_se_rs.seq
+++ b/tests/expected_output/write_epi_se_rs.seq
@@ -202,23 +202,23 @@ TotalDuration 0.21735
 # ..   Hz      ..       ..            ..     us    us     ppm  rad/MHz   Hz   rad  ..
 # Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined
 [RF]
-1      129.712 1 2 0 4000 100 0 0 -424.504 0 u
-2      493.727 3 4 0 1000 130 0 0 -2000 0 u
+1      129.712 1 2 0 4000 100 0 0 -424.504 0 s
+2      493.727 3 4 0 1000 130 0 0 -2000 0 e
 3      987.454 3 4 0 1000 1740 0 0 -2000 1.5708 r
-4      493.727 3 4 0 1000 130 0 0 0 0 u
+4      493.727 3 4 0 1000 130 0 0 0 0 e
 5      987.454 3 4 0 1000 1740 0 0 0 1.5708 r
-6      493.727 3 4 0 1000 130 0 0 2000 0 u
+6      493.727 3 4 0 1000 130 0 0 2000 0 e
 7      987.454 3 4 0 1000 1740 0 0 2000 1.5708 r
 
 # Format of arbitrary gradients:
 #   time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster)
-# id amplitude amp_shape_id time_shape_id delay
-# ..      Hz/m       ..         ..          us
+# id amplitude first last amp_shape_id time_shape_id delay
+# ..      Hz/m  Hz/m Hz/m        ..         ..          us
 [GRADIENTS]
-3  1.34624e+06 5 6 0
-7      -133333 7 8 450
-9      -133333 9 10 0
-10      -133333 11 8 0
+3  1.34624e+06            0            0 5 6 0
+7      -133333            0      -133333 7 8 450
+9      -133333      -133333      -133333 9 10 0
+10      -133333      -133333            0 11 8 0
 
 # Format of trapezoid gradients:
 # id amplitude rise flat fall delay
@@ -232,10 +232,10 @@ TotalDuration 0.21735
  8 -1.00036e+06 220   40 220   0
 
 # Format of ADC events:
-# id num dwell delay freq phase
-# ..  ..    ns    us   Hz   rad
+# id num dwell delay freqPPM phasePPM freq phase phase_id
+# ..  ..    ns    us     ppm  rad/MHz   Hz   rad       ..
 [ADC]
-1 104 3900 37 0 0
+1 104 3900 37 0 0 0 0 0
 
 # Format of extension lists:
 # id type ref next_id
@@ -10336,4 +10336,4 @@ num_samples 2
 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE]
 # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification)
 Type md5
-Hash ae9d2b98e7fc357a1f15425a7fe6d952
+Hash 7754f3ee61085517dc7839965fb2b3fd
diff --git a/tests/expected_output/write_haste.seq b/tests/expected_output/write_haste.seq
index e52550b8..27cded80 100644
--- a/tests/expected_output/write_haste.seq
+++ b/tests/expected_output/write_haste.seq
@@ -189,23 +189,23 @@ TotalDuration 7
 # ..   Hz      ..       ..            ..     us    us     ppm  rad/MHz   Hz   rad  ..
 # Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined
 [RF]
-1      394.982 3 4 0 1250 100 0 0 0 1.5708 u
+1      394.982 3 4 0 1250 100 0 0 0 1.5708 e
 2      987.454 9 10 0 1000 100 0 0 0 0 r
 
 # Format of arbitrary gradients:
 #   time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster)
-# id amplitude amp_shape_id time_shape_id delay
-# ..      Hz/m       ..         ..          us
+# id amplitude first last amp_shape_id time_shape_id delay
+# ..      Hz/m  Hz/m Hz/m        ..         ..          us
 [GRADIENTS]
-1       320000 1 2 0
-2       320000 5 6 0
-3       320000 7 8 0
-5       320000 5 11 0
-6       180372 12 13 0
-8       491667 14 13 0
-9      38940.8 5 15 0
-10       180372 16 13 0
-12       491667 17 13 0
+1       320000            0       320000 1 2 0
+2       320000       320000       320000 5 6 0
+3       320000       320000       320000 7 8 0
+5       320000       320000       320000 5 11 0
+6       180372            0      38940.8 12 13 0
+8       491667       320000            0 14 13 0
+9      38940.8      38940.8      38940.8 5 15 0
+10       180372      38940.8            0 16 13 0
+12       491667            0       320000 17 13 0
 
 # Format of trapezoid gradients:
 # id amplitude rise flat fall delay
@@ -277,10 +277,10 @@ TotalDuration 7
 73     -84092.9 250 1190 250   0
 
 # Format of ADC events:
-# id num dwell delay freq phase
-# ..  ..    ns    us   Hz   rad
+# id num dwell delay freqPPM phasePPM freq phase phase_id
+# ..  ..    ns    us     ppm  rad/MHz   Hz   rad       ..
 [ADC]
-1 64 100000 10 0 0
+1 64 100000 10 0 0 0 0 0
 
 # Sequence Shapes
 [SHAPES]
@@ -4906,4 +4906,4 @@ num_samples 4
 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE]
 # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification)
 Type md5
-Hash eed03ff8e324d1c0f9f0e8e919fdc3eb
+Hash af032808145852407caf689951a9cd32
diff --git a/tests/expected_output/write_mprage.seq b/tests/expected_output/write_mprage.seq
index 282c053c..866ff97b 100644
--- a/tests/expected_output/write_mprage.seq
+++ b/tests/expected_output/write_mprage.seq
@@ -5963,38 +5963,38 @@ TotalDuration 150
 # Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined
 [RF]
 1      563.681 1 2 3 5125 100 0 0 0 0 i
-2      194.444 4 5 6 50 610 0 0 0 0 u
-3      194.444 4 5 6 50 610 0 0 0 2.04204 u
-4      194.444 4 5 6 50 610 0 0 0 6.12611 u
-5      194.444 4 5 6 50 610 0 0 0 5.96903 u
-6      194.444 4 5 6 50 610 0 0 0 1.5708 u
-7      194.444 4 5 6 50 610 0 0 0 5.49779 u
-8      194.444 4 5 6 50 610 0 0 0 5.18363 u
-9      194.444 4 5 6 50 610 0 0 0 0.628319 u
-10      194.444 4 5 6 50 610 0 0 0 4.39823 u
-11      194.444 4 5 6 50 610 0 0 0 3.92699 u
-12      194.444 4 5 6 50 610 0 0 0 2.82743 u
-13      194.444 4 5 6 50 610 0 0 0 2.19911 u
-14      194.444 4 5 6 50 610 0 0 0 3.61283 u
-15      194.444 4 5 6 50 610 0 0 0 0.785398 u
-16      194.444 4 5 6 50 610 0 0 0 1.25664 u
-17      194.444 4 5 6 50 610 0 0 0 4.55531 u
-18      194.444 4 5 6 50 610 0 0 0 4.71239 u
-19      194.444 4 5 6 50 610 0 0 0 0.471239 u
-20      194.444 4 5 6 50 610 0 0 0 1.41372 u
-21      194.444 4 5 6 50 610 0 0 0 3.14159 u
-22      194.444 4 5 6 50 610 0 0 0 5.34071 u
-23      194.444 4 5 6 50 610 0 0 0 2.35619 u
-24      194.444 4 5 6 50 610 0 0 0 3.76991 u
-25      194.444 4 5 6 50 610 0 0 0 2.98451 u
+2      194.444 4 5 6 50 610 0 0 0 0 e
+3      194.444 4 5 6 50 610 0 0 0 2.04204 e
+4      194.444 4 5 6 50 610 0 0 0 6.12611 e
+5      194.444 4 5 6 50 610 0 0 0 5.96903 e
+6      194.444 4 5 6 50 610 0 0 0 1.5708 e
+7      194.444 4 5 6 50 610 0 0 0 5.49779 e
+8      194.444 4 5 6 50 610 0 0 0 5.18363 e
+9      194.444 4 5 6 50 610 0 0 0 0.628319 e
+10      194.444 4 5 6 50 610 0 0 0 4.39823 e
+11      194.444 4 5 6 50 610 0 0 0 3.92699 e
+12      194.444 4 5 6 50 610 0 0 0 2.82743 e
+13      194.444 4 5 6 50 610 0 0 0 2.19911 e
+14      194.444 4 5 6 50 610 0 0 0 3.61283 e
+15      194.444 4 5 6 50 610 0 0 0 0.785398 e
+16      194.444 4 5 6 50 610 0 0 0 1.25664 e
+17      194.444 4 5 6 50 610 0 0 0 4.55531 e
+18      194.444 4 5 6 50 610 0 0 0 4.71239 e
+19      194.444 4 5 6 50 610 0 0 0 0.471239 e
+20      194.444 4 5 6 50 610 0 0 0 1.41372 e
+21      194.444 4 5 6 50 610 0 0 0 3.14159 e
+22      194.444 4 5 6 50 610 0 0 0 5.34071 e
+23      194.444 4 5 6 50 610 0 0 0 2.35619 e
+24      194.444 4 5 6 50 610 0 0 0 3.76991 e
+25      194.444 4 5 6 50 610 0 0 0 2.98451 e
 
 # Format of arbitrary gradients:
 #   time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster)
-# id amplitude amp_shape_id time_shape_id delay
-# ..      Hz/m       ..         ..          us
+# id amplitude first last amp_shape_id time_shape_id delay
+# ..      Hz/m  Hz/m Hz/m        ..         ..          us
 [GRADIENTS]
-4      -708063 7 8 0
-5       984720 9 10 0
+4      -708063            0      49824.6 7 8 0
+5       984720      49824.6            0 9 10 0
 
 # Format of trapezoid gradients:
 # id amplitude rise flat fall delay
@@ -6100,33 +6100,33 @@ TotalDuration 150
 100     -23148.1 180    0 180   0
 
 # Format of ADC events:
-# id num dwell delay freq phase
-# ..  ..    ns    us   Hz   rad
+# id num dwell delay freqPPM phasePPM freq phase phase_id
+# ..  ..    ns    us     ppm  rad/MHz   Hz   rad       ..
 [ADC]
-1 64 78400 380 0 0
-2 64 78400 380 0 2.04204
-3 64 78400 380 0 6.12611
-4 64 78400 380 0 5.96903
-5 64 78400 380 0 1.5708
-6 64 78400 380 0 5.49779
-7 64 78400 380 0 5.18363
-8 64 78400 380 0 0.628319
-9 64 78400 380 0 4.39823
-10 64 78400 380 0 3.92699
-11 64 78400 380 0 2.82743
-12 64 78400 380 0 2.19911
-13 64 78400 380 0 3.61283
-14 64 78400 380 0 0.785398
-15 64 78400 380 0 1.25664
-16 64 78400 380 0 4.55531
-17 64 78400 380 0 4.71239
-18 64 78400 380 0 0.471239
-19 64 78400 380 0 1.41372
-20 64 78400 380 0 3.14159
-21 64 78400 380 0 5.34071
-22 64 78400 380 0 2.35619
-23 64 78400 380 0 3.76991
-24 64 78400 380 0 2.98451
+1 64 78400 380 0 0 0 0 0
+2 64 78400 380 0 0 0 2.04204 0
+3 64 78400 380 0 0 0 6.12611 0
+4 64 78400 380 0 0 0 5.96903 0
+5 64 78400 380 0 0 0 1.5708 0
+6 64 78400 380 0 0 0 5.49779 0
+7 64 78400 380 0 0 0 5.18363 0
+8 64 78400 380 0 0 0 0.628319 0
+9 64 78400 380 0 0 0 4.39823 0
+10 64 78400 380 0 0 0 3.92699 0
+11 64 78400 380 0 0 0 2.82743 0
+12 64 78400 380 0 0 0 2.19911 0
+13 64 78400 380 0 0 0 3.61283 0
+14 64 78400 380 0 0 0 0.785398 0
+15 64 78400 380 0 0 0 1.25664 0
+16 64 78400 380 0 0 0 4.55531 0
+17 64 78400 380 0 0 0 4.71239 0
+18 64 78400 380 0 0 0 0.471239 0
+19 64 78400 380 0 0 0 1.41372 0
+20 64 78400 380 0 0 0 3.14159 0
+21 64 78400 380 0 0 0 5.34071 0
+22 64 78400 380 0 0 0 2.35619 0
+23 64 78400 380 0 0 0 3.76991 0
+24 64 78400 380 0 0 0 2.98451 0
 
 # Format of extension lists:
 # id type ref next_id
@@ -8263,4 +8263,4 @@ num_samples 4
 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE]
 # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification)
 Type md5
-Hash 92db8dfb899437b8caab1977b9c478f0
+Hash aa9547958f21808d0bc90c8b0ba90fda
diff --git a/tests/expected_output/write_tse.seq b/tests/expected_output/write_tse.seq
index cb546032..3704c86b 100644
--- a/tests/expected_output/write_tse.seq
+++ b/tests/expected_output/write_tse.seq
@@ -372,23 +372,23 @@ TotalDuration 10
 # ..   Hz      ..       ..            ..     us    us     ppm  rad/MHz   Hz   rad  ..
 # Field "use" is the initial of: excitation refocusing inversion saturation preparation other undefined
 [RF]
-1      394.982 3 4 0 1250 100 0 0 0 1.5708 u
+1      394.982 3 4 0 1250 100 0 0 0 1.5708 e
 2      987.454 9 10 0 1000 100 0 0 0 0 r
 
 # Format of arbitrary gradients:
 #   time_shape_id of 0 means default timing (stepping with grad_raster starting at 1/2 of grad_raster)
-# id amplitude amp_shape_id time_shape_id delay
-# ..      Hz/m       ..         ..          us
+# id amplitude first last amp_shape_id time_shape_id delay
+# ..      Hz/m  Hz/m Hz/m        ..         ..          us
 [GRADIENTS]
-1       320000 1 2 0
-2       320000 5 6 0
-3       320000 7 8 0
-5       320000 5 11 0
-6       180372 12 13 0
-8       491667 14 13 0
-9      38940.8 5 15 0
-10       180372 16 13 0
-11       491667 17 13 0
+1       320000            0       320000 1 2 0
+2       320000       320000       320000 5 6 0
+3       320000       320000       320000 7 8 0
+5       320000       320000       320000 5 11 0
+6       180372            0      38940.8 12 13 0
+8       491667       320000            0 14 13 0
+9      38940.8      38940.8      38940.8 5 15 0
+10       180372      38940.8            0 16 13 0
+11       491667            0       320000 17 13 0
 
 # Format of trapezoid gradients:
 # id amplitude rise flat fall delay
@@ -462,10 +462,10 @@ TotalDuration 10
 75      86805.6 250 1190 250   0
 
 # Format of ADC events:
-# id num dwell delay freq phase
-# ..  ..    ns    us   Hz   rad
+# id num dwell delay freqPPM phasePPM freq phase phase_id
+# ..  ..    ns    us     ppm  rad/MHz   Hz   rad       ..
 [ADC]
-1 64 100000 10 0 0
+1 64 100000 10 0 0 0 0 0
 
 # Sequence Shapes
 [SHAPES]
@@ -5091,4 +5091,4 @@ num_samples 4
 # It can be reproduced/verified with md5sum if the file trimmed to the position right above [SIGNATURE]
 # The new line character preceding [SIGNATURE] BELONGS to the signature (and needs to be stripped away for recalculating/verification)
 Type md5
-Hash f8003a45f3ee6e0f4999aab132c804b2
+Hash 143ed0af8842a255835442c1ba7df1e7
diff --git a/tests/test_scale_grad.py b/tests/test_scale_grad.py
new file mode 100644
index 00000000..7d98170f
--- /dev/null
+++ b/tests/test_scale_grad.py
@@ -0,0 +1,132 @@
+import numpy as np
+import pypulseq as pp
+import pytest
+from pypulseq import scale_grad
+from pypulseq.opts import Opts
+
+# Updated gradients with realistic hardware limits
+grad_list = [
+    pp.make_trapezoid(channel='x', amplitude=10, duration=13, max_grad=30, max_slew=200),
+    pp.make_trapezoid(channel='y', amplitude=10, duration=13, max_grad=30, max_slew=200),
+    pp.make_trapezoid(channel='z', amplitude=10, duration=13, max_grad=30, max_slew=200),
+    pp.make_trapezoid(channel='x', amplitude=20, duration=5, max_grad=25, max_slew=150),
+    pp.make_trapezoid(channel='y', amplitude=20, duration=5, max_grad=25, max_slew=150),
+    pp.make_trapezoid(channel='z', amplitude=20, duration=5, max_grad=25, max_slew=150),
+    pp.make_extended_trapezoid(
+        'x', [0, 15, 5, 10], convert_to_arbitrary=True, times=[1, 3, 4, 7], max_grad=40, max_slew=300
+    ),
+    pp.make_extended_trapezoid(
+        'y', [0, 15, 5, 10], convert_to_arbitrary=True, times=[1, 3, 4, 7], max_grad=40, max_slew=300
+    ),
+    pp.make_extended_trapezoid(
+        'z', [0, 15, 5, 10], convert_to_arbitrary=True, times=[1, 3, 4, 7], max_grad=40, max_slew=300
+    ),
+    pp.make_extended_trapezoid(
+        'x', [0, 20, 10, 15], convert_to_arbitrary=False, times=[1, 3, 4, 7], max_grad=25, max_slew=150
+    ),
+    pp.make_extended_trapezoid(
+        'y', [0, 20, 10, 15], convert_to_arbitrary=False, times=[1, 3, 4, 7], max_grad=25, max_slew=150
+    ),
+    pp.make_extended_trapezoid(
+        'z', [0, 20, 10, 15], convert_to_arbitrary=False, times=[1, 3, 4, 7], max_grad=25, max_slew=150
+    ),
+    pp.make_extended_trapezoid(
+        'x', [0, 10, 5, 10], convert_to_arbitrary=False, times=[1, 2, 3, 4], max_grad=15, max_slew=80
+    ),
+    pp.make_extended_trapezoid(
+        'y', [0, 10, 5, 10], convert_to_arbitrary=False, times=[1, 2, 3, 4], max_grad=15, max_slew=80
+    ),
+    pp.make_extended_trapezoid(
+        'z', [0, 10, 5, 10], convert_to_arbitrary=False, times=[1, 2, 3, 4], max_grad=15, max_slew=80
+    ),
+]
+
+# ----------- TEST SCALING IS CORRECT -----------
+
+
+@pytest.mark.parametrize('grad', grad_list)
+def test_scale_grad_correct_scaling(grad):
+    scale = 0.5  # Safe scale
+    system = Opts(max_grad=40, max_slew=300)
+    scaled = scale_grad(grad, scale, system)
+
+    if grad.type == 'trap':
+        assert np.isclose(scaled.amplitude, grad.amplitude * scale)
+        assert np.isclose(scaled.flat_area, grad.flat_area * scale)
+    else:
+        assert np.allclose(scaled.waveform, grad.waveform * scale)
+        assert np.isclose(scaled.first, grad.first * scale)
+        assert np.isclose(scaled.last, grad.last * scale)
+
+    assert np.isclose(scaled.area, grad.area * scale)
+
+    # ID must be removed
+    assert not hasattr(scaled, 'id')
+
+
+# ----------- TEST AMPLITUDE VIOLATION -----------
+
+
+def test_scale_grad_violates_amplitude():
+    scale = 100.0
+    system = Opts(max_grad=40, max_slew=999999999)  # make sure we have infinite slew rate
+
+    expected_failures = 0
+    actual_failures = 0
+
+    for grad in grad_list:
+        if grad.type == 'trap':
+            should_fail = abs(grad.amplitude) * scale > system.max_grad
+        else:
+            should_fail = max(abs(grad.waveform)) * scale > system.max_grad
+
+        if should_fail:
+            expected_failures += 1
+            with pytest.raises(ValueError, match='maximum amplitude exceeded'):
+                scale_grad(grad, scale, system)
+            actual_failures += 1
+        else:
+            scale_grad(grad, scale, system)
+
+    assert expected_failures == actual_failures, (
+        f'Expected {expected_failures} amplitude violations, but got {actual_failures}.'
+    )
+
+
+# ----------- TEST SLEW RATE VIOLATION -----------
+
+
+def test_scale_grad_violates_slew():
+    scale = 100.0
+    system = Opts(max_grad=999999999, max_slew=300)  # make sure we have infinite gradient amp
+
+    expected_failures = 0
+    actual_failures = 0
+
+    for grad in grad_list:
+        if grad.type == 'trap':
+            if abs(grad.amplitude) > 1e-6:
+                approx_slew = abs(grad.amplitude * scale) / min(grad.rise_time, grad.fall_time)
+                should_fail = approx_slew > system.max_slew
+            else:
+                should_fail = False
+        else:
+            if max(abs(grad.waveform)) > 1e-6:
+                waveform = np.array(grad.waveform) * scale
+                tt = np.array(grad.tt)
+                diffs = np.abs(np.diff(waveform) / np.diff(tt))
+                should_fail = max(diffs) > system.max_slew
+            else:
+                should_fail = False
+
+        if should_fail:
+            expected_failures += 1
+            with pytest.raises(ValueError, match='maximum slew rate exceeded'):
+                scale_grad(grad, scale, system)
+            actual_failures += 1
+        else:
+            scale_grad(grad, scale, system)
+
+    assert expected_failures == actual_failures, (
+        f'Expected {expected_failures} slew violations, but got {actual_failures}.'
+    )
diff --git a/tests/test_sequence.py b/tests/test_sequence.py
index 4f07cfc4..b8047a6a 100644
--- a/tests/test_sequence.py
+++ b/tests/test_sequence.py
@@ -215,14 +215,14 @@ def seq4():
 seq_examples = [
     'write_gre',
     'write_gre_label',
-    # 'write_haste', # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last')
+    'write_haste',  # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last')
     'write_radial_gre',
-    # 'write_tse', # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last')
+    'write_tse',  # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last')
     'write_epi',
     'write_epi_label',
     'write_epi_se',
-    # 'write_epi_se_rs', # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last')
-    # 'write_mprage', # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last')
+    'write_epi_se_rs',  # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last')
+    'write_mprage',  # TODO: re-enable when bumping grad storage to v1.5.x (i.e., 'first', 'last')
     'write_ute',
 ]