Support stimulus_template as optional column in IntracellularRecordingsTable

CHANGELOG.md

@@ -6,6 +6,7 @@
 - Added support for NWB schema 2.7.0.
   - ...
   - Modified `OptogeneticSeries` to allow 2D data, primarily in extensions of `OptogeneticSeries`. @rly [#1812](https://github.com/NeurodataWithoutBorders/pynwb/pull/1812)
+  - Support `stimulus_template` as optional predefined column in `IntracellularStimuliTable`. @stephprince [#1815](https://github.com/NeurodataWithoutBorders/pynwb/pull/1815)
   - ...
   - ...
 - For `NWBHDF5IO()`, change the default of arg `load_namespaces` from `False` to `True`. @bendichter [#1748](https://github.com/NeurodataWithoutBorders/pynwb/pull/1748)

docs/gallery/domain/plot_icephys.py

@@ -97,6 +97,7 @@
 
 # Import additional core datatypes used in the example
 from pynwb.core import DynamicTable, VectorData
+from pynwb.base import TimeSeriesReference, TimeSeriesReferenceVectorData
 
 # Import icephys TimeSeries types used
 from pynwb.icephys import VoltageClampSeries, VoltageClampStimulusSeries
@@ -457,6 +458,59 @@
     category="electrodes",
 )
 
+#####################################################################
+# Adding stimulus templates
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+#
+# One predefined subcategory column is the ``stimulus_template`` column in the stimuli table. This column is
+# used to store template waveforms of stimuli in addition to the actual recorded stimulus that is stored in the
+# ``stimulus`` column. The ``stimulus_template`` column contains an idealized version of the template waveform used as
+# the stimulus. This can be useful as a noiseless version of the stimulus for data analysis or to validate that the
+# recorded stimulus matches the expected waveform of the template. Similar to the ``stimulus`` and ``response``
+# columns, we can specify a relevant time range.
+
+stimulus_template = VoltageClampStimulusSeries(
+    name="ccst",
+    data=[0, 1, 2, 3, 4],
+    starting_time=0.0,
+    rate=10e3,
+    electrode=electrode,
+    gain=0.02,
+)
+nwbfile.add_stimulus_template(stimulus_template)
+
+nwbfile.intracellular_recordings.add_column(
+    name="stimulus_template",
+    data=[TimeSeriesReference(0, 5, stimulus_template),  # (start_index, index_count, stimulus_template)
+          TimeSeriesReference(1, 3, stimulus_template),
+          TimeSeriesReference.empty(stimulus_template)],  # if there was no data for that recording, use empty reference
+    description="Column storing the reference to the stimulus template for the recording (rows).",
+    category="stimuli",
+    col_cls=TimeSeriesReferenceVectorData
+)
+
+# we can also add stimulus template data as follows
+rowindex = nwbfile.add_intracellular_recording(
+    electrode=electrode,
+    stimulus=stimulus,
+    stimulus_template=stimulus_template,  # the full time range of the stimulus template will be used unless specified
+    recording_tag='A4',
+    recording_lab_data={'location': 'Isengard'},
+    electrode_metadata={'voltage_threshold': 0.14},
+    id=13,
+)
+
+#####################################################################
+# .. note:: If a stimulus template column exists but there is no stimulus template data for that recording, then
+#           :py:meth:`~pynwb.file.NWBFile.add_intracellular_recording` will internally set the stimulus template to the
+#           provided stimulus or response TimeSeries and the start_index and index_count for the missing parameter are
+#           set to -1. The missing values will be represented via masked numpy arrays.
+
+#####################################################################
+# .. note:: Since stimulus templates are often reused across many recordings, the timestamps in the templates are not
+#           usually aligned with the recording nor with the reference time of the file. The timestamps often start
+#           at 0 and are relative to the time of the application of the stimulus.
+
 #####################################################################
 # Add a simultaneous recording
 # ---------------------------------

src/pynwb/base.py

@@ -497,6 +497,26 @@ def data(self):
         # load the data from the timeseries
         return self.timeseries.data[self.idx_start: (self.idx_start + self.count)]
 
+    @classmethod
+    @docval({'name': 'timeseries', 'type': TimeSeries, 'doc': 'the timeseries object to reference.'})
+    def empty(cls, timeseries):
+        """
+        Creates an empty TimeSeriesReference object to represent missing data.
+
+        When missing data needs to be represented, NWB defines ``None`` for the complex data type ``(idx_start,
+        count, TimeSeries)`` as (-1, -1, TimeSeries) for storage. The exact timeseries object will technically not
+        matter since the empty reference is a way of indicating a NaN value in a
+        :py:class:`~pynwb.base.TimeSeriesReferenceVectorData` column.
+
+        An example where this functionality is used is :py:class:`~pynwb.icephys.IntracellularRecordingsTable`
+        where only one of stimulus or response data was recorded. In such cases, the timeseries object for the
+        empty stimulus :py:class:`~pynwb.base.TimeSeriesReference` could be set to the response series, or vice versa.
+
+        :returns: Returns :py:class:`~pynwb.base.TimeSeriesReference`
+        """
+
+        return cls(-1, -1, timeseries)
+
 
 @register_class('TimeSeriesReferenceVectorData', CORE_NAMESPACE)
 class TimeSeriesReferenceVectorData(VectorData):

src/pynwb/icephys.py

@@ -415,6 +415,12 @@ class IntracellularStimuliTable(DynamicTable):
          'index': False,
          'table': False,
          'class': TimeSeriesReferenceVectorData},
+        {'name': 'stimulus_template',
+         'description': 'Column storing the reference to the stimulus template for the recording (rows)',
+         'required': False,
+         'index': False,
+         'table': False,
+         'class': TimeSeriesReferenceVectorData},
     )
 
     @docval(*get_docval(DynamicTable.__init__, 'id', 'columns', 'colnames'))
@@ -518,6 +524,13 @@ def __init__(self, **kwargs):
         {'name': 'stimulus', 'type': TimeSeries,
          'doc': 'The TimeSeries (usually a PatchClampSeries) with the stimulus',
          'default': None},
+        {'name': 'stimulus_template_start_index', 'type': int, 'doc': 'Start index of the stimulus template',
+         'default': None},
+        {'name': 'stimulus_template_index_count', 'type': int, 'doc': 'Stop index of the stimulus template',
+         'default': None},
+        {'name': 'stimulus_template', 'type': TimeSeries,
+         'doc': 'The TimeSeries (usually a PatchClampSeries) with the stimulus template waveforms',
+         'default': None},
         {'name': 'response_start_index', 'type': int, 'doc': 'Start index of the response', 'default': None},
         {'name': 'response_index_count', 'type': int, 'doc': 'Stop index of the response', 'default': None},
         {'name': 'response', 'type': TimeSeries,
@@ -553,6 +566,11 @@ def add_recording(self, **kwargs):
                                                                        'response',
                                                                        kwargs)
         electrode = popargs('electrode', kwargs)
+        stimulus_template_start_index, stimulus_template_index_count, stimulus_template = popargs(
+            'stimulus_template_start_index',
+            'stimulus_template_index_count',
+            'stimulus_template',
+            kwargs)
 
         # if electrode is not provided, take from stimulus or response object
         if electrode is None:
@@ -572,6 +590,15 @@ def add_recording(self, **kwargs):
         response_start_index, response_index_count = self.__compute_index(response_start_index,
                                                                           response_index_count,
                                                                           response, 'response')
+        stimulus_template_start_index, stimulus_template_index_count = self.__compute_index(
+            stimulus_template_start_index,
+            stimulus_template_index_count,
+            stimulus_template, 'stimulus_template')
+
+        # if stimulus template is already a column in the stimuli table, but stimulus_template was None
+        if 'stimulus_template' in self.category_tables['stimuli'].colnames and stimulus_template is None:
+            stimulus_template = stimulus if stimulus is not None else response  # set to stimulus if it was provided
+
         # If either stimulus or response are None, then set them to the same TimeSeries to keep the I/O happy
         response = response if response is not None else stimulus
         stimulus_provided_is_not_none = stimulus is not None  # Store if stimulus is None for error checks later
@@ -612,6 +639,9 @@ def add_recording(self, **kwargs):
             stimuli = {}
         stimuli['stimulus'] = TimeSeriesReferenceVectorData.TIME_SERIES_REFERENCE_TUPLE(
             stimulus_start_index, stimulus_index_count, stimulus)
+        if stimulus_template is not None:
+            stimuli['stimulus_template'] = TimeSeriesReferenceVectorData.TIME_SERIES_REFERENCE_TUPLE(
+                stimulus_template_start_index, stimulus_template_index_count, stimulus_template)
 
         # Compile the responses table data
         responses = copy(popargs('response_metadata', kwargs))

tests/unit/test_base.py

@@ -877,3 +877,9 @@ def test_data_property_bad_reference(self):
             IndexError, "'idx_start + count' out of range for timeseries 'test'"
         ):
             tsr.data
+
+    def test_empty_reference_creation(self):
+        tsr = TimeSeriesReference.empty(self._create_time_series_with_rate())
+        self.assertFalse(tsr.isvalid())
+        self.assertIsNone(tsr.data)
+        self.assertIsNone(tsr.timestamps)

tests/unit/test_icephys_metadata_tables.py

@@ -419,6 +419,15 @@ def test_add_row_index_out_of_range(self):
                 response=self.response,
                 id=np.int64(10)
             )
+        with self.assertRaises(IndexError):
+            ir = IntracellularRecordingsTable()
+            ir.add_recording(
+                electrode=self.electrode,
+                stimulus_template=self.stimulus,
+                stimulus_template_start_index=10,
+                response=self.response,
+                id=np.int64(10)
+            )
         # Stimulus/Response index count too large
         with self.assertRaises(IndexError):
             ir = IntracellularRecordingsTable()
@@ -438,6 +447,15 @@ def test_add_row_index_out_of_range(self):
                 response=self.response,
                 id=np.int64(10)
             )
+        with self.assertRaises(IndexError):
+            ir = IntracellularRecordingsTable()
+            ir.add_recording(
+                electrode=self.electrode,
+                stimulus_template=self.stimulus,
+                stimulus_template_index_count=10,
+                response=self.response,
+                id=np.int64(10)
+            )
         # Stimulus/Response start+count combination too large
         with self.assertRaises(IndexError):
             ir = IntracellularRecordingsTable()
@@ -459,6 +477,16 @@ def test_add_row_index_out_of_range(self):
                 response=self.response,
                 id=np.int64(10)
             )
+        with self.assertRaises(IndexError):
+            ir = IntracellularRecordingsTable()
+            ir.add_recording(
+                electrode=self.electrode,
+                stimulus_template=self.stimulus,
+                stimulus_template_start_index=3,
+                stimulus_template_index_count=4,
+                response=self.response,
+                id=np.int64(10)
+            )
 
     def test_add_row_no_stimulus_and_response(self):
         with self.assertRaises(ValueError):
@@ -469,6 +497,40 @@ def test_add_row_no_stimulus_and_response(self):
                 response=None
             )
 
+    def test_add_row_with_stimulus_template(self):
+        ir = IntracellularRecordingsTable()
+        ir.add_recording(
+            electrode=self.electrode,
+            stimulus=self.stimulus,
+            stimulus_template=self.stimulus,
+            response=self.response,
+            id=np.int64(10)
+        )
+
+    def test_add_stimulus_template_column(self):
+        ir = IntracellularRecordingsTable()
+        ir.add_column(name='stimulus_template',
+                      description='test column',
+                      category='stimuli',
+                      col_cls=TimeSeriesReferenceVectorData)
+
+    def test_add_row_with_no_stimulus_template_when_stimulus_template_column_exists(self):
+        ir = IntracellularRecordingsTable()
+        ir.add_recording(electrode=self.electrode,
+                         stimulus=self.stimulus,
+                         response=self.response,
+                         stimulus_template=self.stimulus,
+                         id=np.int64(10))
+
+        # add row with only stimulus when stimulus template column already exists
+        ir.add_recording(electrode=self.electrode,
+                         stimulus=self.stimulus,
+                         id=np.int64(20))
+        # add row with only response when stimulus template column already exists
+        ir.add_recording(electrode=self.electrode,
+                         response=self.stimulus,
+                         id=np.int64(30))
+
     def test_add_column(self):
         ir = IntracellularRecordingsTable()
         ir.add_recording(