Repository: myokit/myokit (language: python)

Issue:
Multiple protocols
There's some reasons it might be good to have multiple protocols:

- Clamping both membrane potential and a concentration to some time signal (e.g. Vm and Ca)

This would require quite major changes in Myokit.

~One way to implement this might be to remove the "binding" construct, and replace it with some special function (`myokit.Expression`), e.g.~

```
I_stim = input('pace') * stim_amplitude
I_diff = input('diffusion_current')
```
~where 'diffusion_current' would be a name set by the simulation (but customisable?), while the name 'pace' would be set by the user, e.g. by having a `name: pace` field in the protocol definition, or a constructor like `p = myokit.Protocol('pace')`.~

A different way to implement it would be to set it in the simulation, e.g.

```
s.set_protocol('pace', p)
s.set_protocol('other_label', q)
```
`Model_SetLiteralVariables` not used
function in `cmodel.h` doesn't seem to be used anywhere


Discussion / hints:
A "least resistance" way in could be to start with #386, then go
```
def set_protocol(protocol, label='pace')
```
but don't allow more than one protocol initially. Then it's up to simulations to decide if they want to support multiple, but there's a clear way to do it without model-syntax changes

A major difficulty with this is that it is no longer known at compile time whether a variable is "bound" or not. (Even if we would have e.g. `bind pace1`, `bind pace2` etc., the usual step is to _remove all unused bindings before compilation_). So we need to rewrite the simulations (and exports) so that they don't expect to know which variables will be used to input solver variables at compilation time

Roadmap:

- [ ] Get rid of the "binding" concept, replace it with labels
  - [ ] Remove `prepare_bindings` (and need for it) in simulations https://github.com/MichaelClerx/myokit/issues/601
  - [ ] Scan for is_bound checks listed here and figure out how to avoid them https://github.com/MichaelClerx/myokit/issues/386. Might also involve this one: https://github.com/MichaelClerx/myokit/issues/832
  - [ ] Update API docs https://github.com/MichaelClerx/myokit/issues/386
  - [ ] Update external examples https://github.com/MichaelClerx/myokit/issues/833
- [ ] Allow label to be specified, make sure sims expect this to be something that can change repeatedly post-compilation https://github.com/MichaelClerx/myokit/issues/732
- [ ] This issue https://github.com/MichaelClerx/myokit/issues/320 will probably disappear once we've done this?

Hi @MichaelClerx, I'd be happy to work on this if this would be helpful? Anywhere in particular I should start? should I follow the roadmap above, starting at #601 ?
That would be great if you wanted! The roadmap still makes sense to me. Might be good to start with removing prepare_bindings in some exports, like ansi c etc. It will be mostly older code that uses it.




Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 195eeed55..ea0c7c35a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,8 @@ This page lists the main changes made to Myokit in each release.
   - [#918](https://github.com/myokit/myokit/pull/918) Added `Variable.initial_value`, which returns a state variable's initial value as an expression or float.
   - [#918](https://github.com/myokit/myokit/pull/918) Added `Variable.set_initial_value`, which sets a state variable's initial value using an expression, string expression, or float.
   - [#918](https://github.com/myokit/myokit/pull/918) Added `IllegalReferenceInInitialValueError` for illegal references made in initial value expressions.
+  - [#947](https://github.com/myokit/myokit/pull/947) The `Simulation` class can now handle multiple protocols during a single simulation.
+  - [#947](https://github.com/myokit/myokit/pull/947) Fixed-form protocols are now implemented in a `TimeSeriesProtocol`.
 - Changed
   - [#898](https://github.com/myokit/myokit/pull/899) `Variable.promote` can now take an expression, a number, or a string (which will be parsed as an expression) for the argument `initial_value`.
   - [#909](https://github.com/MichaelClerx/myokit/pull/909) Myokit is now tested on Python 3.11.
@@ -25,14 +27,17 @@ This page lists the main changes made to Myokit in each release.
   - [#918](https://github.com/myokit/myokit/pull/918) `Variable.set_state_value` is deprecated in favour of `Variable.set_initial_value`.
   - [#918](https://github.com/myokit/myokit/pull/918) `Model.save_state` and `Model.load_state` are deprecated.
   - [#918](https://github.com/myokit/myokit/pull/918) The `Variable.promote` keyword argument `state_value` is deprecated in favour of `initial_value`.
+  - [#947](https://github.com/myokit/myokit/pull/947) `Simulation.set_fixed_form_protocol` is deprecated in favour of using `Simulation.set_protocol` with a `TimeSeriesProtocol`
 - Removed
   - [#898](https://github.com/myokit/myokit/pull/899) `NonLiteralValueError` is now no longer raised by Myokit and so has been removed.
   - [#921](https://github.com/myokit/myokit/pull/921) Removed deprecated `LegacySimulation` class.
   - [#921](https://github.com/myokit/myokit/pull/921) Removed deprecated `PSimulation` class.
   - [#921](https://github.com/myokit/myokit/pull/921) Removed deprecated `ICSimulation` class.
+  - [#957](https://github.com/myokit/myokit/pull/957) Removed the method `Model.prepare_bindings`.
 - Fixed
   - [#918](https://github.com/myokit/myokit/pull/918) Fixed error when `var()` was called on a `PartialDerivative` with a `Derivative` as first argument.
   - [#933](https://github.com/myokit/myokit/pull/933) Made `myokit.step` base its error classification (large/small/none) on numerics, not representation.
+  - [#947](https://github.com/myokit/myokit/pull/947) Fixed minor memory leak (reference counting issue) when using fixed form protocols.
 
 ## [1.33.9] - 2022-11-08
 - Fixed
diff --git a/docs/source/api_index/myokit.rst b/docs/source/api_index/myokit.rst
index c05c342e3..f64a22d17 100644
--- a/docs/source/api_index/myokit.rst
+++ b/docs/source/api_index/myokit.rst
@@ -155,6 +155,7 @@ myokit
 - :class:`myokit.Tan`
 - :meth:`myokit.time`
 - :class:`myokit.Timeout`
+- :class:`myokit.TimeSeriesProtocol`
 - :class:`myokit.Unit`
 - :class:`myokit.UnresolvedReferenceError`
 - :class:`myokit.UnusedVariableError`
diff --git a/docs/source/api_simulations/Protocol.rst b/docs/source/api_simulations/Protocol.rst
index e00f700e1..47752d950 100644
--- a/docs/source/api_simulations/Protocol.rst
+++ b/docs/source/api_simulations/Protocol.rst
@@ -12,6 +12,9 @@ Pacing protocols
 
 .. autoclass:: PacingSystem
 
+.. _api/myokit.TimeSeriesProtocol:
+.. autoclass:: TimeSeriesProtocol
+
 Protocol factory
 ================
 
diff --git a/myokit/__init__.py b/myokit/__init__.py
index c12f90b05..d0ec2f86f 100644
--- a/myokit/__init__.py
+++ b/myokit/__init__.py
@@ -462,6 +462,7 @@
     PacingSystem,
     Protocol,
     ProtocolEvent,
+    TimeSeriesProtocol,
 )
 from . import pacing  # noqa
 
@@ -505,6 +506,7 @@
     default_script,
     ModelComparison,
     numpy_writer,
+    _prepare_bindings,
     python_writer,
     run,
     step,
diff --git a/myokit/_aux.py b/myokit/_aux.py
index 6a0732d7b..4936c5945 100644
--- a/myokit/_aux.py
+++ b/myokit/_aux.py
@@ -434,6 +434,50 @@ def name(x):
     return _pywriter_
 
 
+def _prepare_bindings(model, labels):
+    """
+    Takes a mapping of binding labels to internal references and returns a map
+    from variable instances to the same internal references.
+
+    This method also modifies ``model``: All bindings that are not mapped to
+    an internal reference will be removed.
+
+    The argument ``mapping`` should take the form::
+
+        labels = {
+            binding_label_1 : internal_name_1,
+            binding_label_2 : internal_name_2,
+            ...
+            }
+
+    The returned dictionary will have the form::
+
+        variables = {
+            variable_x : internal_name_1,
+            variable_y : internal_name_2,
+            ...
+            }
+
+    Unsupported bindings (i.e. bindings not appearing in ``labels``) will
+    be ignored.
+    """
+    unused = []
+    variables = {}
+    for label, var in model.bindings():
+        try:
+            variables[var] = labels[label]
+        except KeyError:
+            unused.append(var)
+            continue
+        # TODO: Remove this line; https://github.com/myokit/myokit/issues/320
+        var.set_rhs(0)
+
+    for var in unused:
+        var.set_binding(None)
+
+    return variables
+
+
 def run(model, protocol, script, stdout=None, stderr=None, progress=None):
     """
     Runs a python ``script`` using the given ``model`` and ``protocol``.
diff --git a/myokit/_model_api.py b/myokit/_model_api.py
index 7cf9b6220..c44325e21 100644
--- a/myokit/_model_api.py
+++ b/myokit/_model_api.py
@@ -2634,45 +2634,6 @@ def name(self):
         except KeyError:
             return None
 
-    def prepare_bindings(self, labels):
-        """
-        Takes a mapping of binding labels to internal references as input and
-        returns a mapping of variables to internal references. All variables
-        appearing in the map will have their right hand side set to zero. All
-        bindings not mapped to any internal reference will be deleted.
-
-        The argument ``mapping`` should take the form::
-
-            labels = {
-                'binding_label_1' : internal_name_1,
-                'binding_label_2' : internal_name_2,
-                ...
-                }
-
-        The returned dictionary will have the form::
-
-            variables = {
-                variable_x : internal_name_1,
-                variable_y : internal_name_2,
-                ...
-                }
-
-        Unsupported bindings (i.e. bindings not appearing in ``labels``) will
-        be ignored.
-        """
-        unused = []
-        variables = {}
-        for label, var in self._bindings.items():
-            try:
-                variables[var] = labels[label]
-            except KeyError:
-                unused.append(var)
-                continue
-            var.set_rhs(0)
-        for var in unused:
-            var.set_binding(None)
-        return variables
-
     def __reduce__(self):
         """
         Pickles the model.
@@ -3461,7 +3422,7 @@ def user_functions(self):
 
     def validate(self, remove_unused_variables=False):
         """
-        Attempts to check model validity, raises errors if it isn't.
+        Validates this model and raises errors if any issues are found.
 
         Small issues (e.g. unused variables) will generate warnings, which
         can be retrieved using :meth:`Model.warnings()` or
diff --git a/myokit/_protocol.py b/myokit/_protocol.py
index 6c6267ca4..d3ed11428 100644
--- a/myokit/_protocol.py
+++ b/myokit/_protocol.py
@@ -1,5 +1,5 @@
 #
-# Defines the python classes that represent a pacing protocol.
+# Defines the python classes that represent pacing protocols.
 #
 # This file is part of Myokit.
 # See http://myokit.org for copyright, sharing, and licensing details.
@@ -11,6 +11,8 @@
 
 import numpy as np
 
+from bisect import bisect_right
+
 
 class Protocol(object):
     """
@@ -850,3 +852,94 @@ class NotAnUnbrokenSequenceError(myokit.MyokitError):
     """ Error raised exclusively by is_unbroken_sequence_exception(). """
     pass
 
+
+class TimeSeriesProtocol(object):
+    """
+    Represents a pacing protocol as a sequence of time value pairs and an
+    interpolation method (currently only linear interpolation is supported).
+
+    A 1D time-series should be given as input. During the simulation, the value
+    of the pacing variable will be determined by interpolating between the two
+    nearest points in the series. If the simulation time is outside the bounds
+    of the time-series, the first or last value in the series will be used.
+
+    Protocols can be compared with ``==``, which will check if the sequence of
+    time value pairs is the same, and the interpolation method is the same.
+    Protocols can be serialised with ``pickle``.
+
+    """
+
+    def __init__(self, times, values, method=None):
+        super(TimeSeriesProtocol, self).__init__()
+
+        if len(times) != len(values):
+            raise ValueError('Times and values array must have same size.')
+        times_tpl, values_tpl = zip(*[
+            (float(t), float(v)) for t, v in sorted(zip(times, values))
+        ])
+        self._times = list(times_tpl)
+        self._values = list(values_tpl)
+
+        if method is None:
+            self._method = 'linear'
+        else:
+            self._method = str(method).lower()
+            if self._method not in ('linear', ):
+                raise ValueError(
+                    'Unknown interpolation method: ' + self._method)
+
+    def __eq__(self, other):
+        if self is other:
+            return True
+        if not isinstance(other, TimeSeriesProtocol):
+            return False
+        return (
+            self._method == other._method
+            and self._values == other._values
+            and self._times == other._times
+        )
+
+    def __getstate__(self):
+        return {
+            'times': self._times,
+            'values': self._values,
+            'method': self._method,
+        }
+
+    def __setstate__(self, values):
+        self._times = values['times']
+        self._values = values['values']
+        self._method = values['method']
+
+    def clone(self):
+        """
+        Returns a clone of this protocol.
+        """
+        return TimeSeriesProtocol(self._times, self._values, self._method)
+
+    def pace(self, t):
+        """
+        Returns the value of the pacing variable at time ``t``.
+        """
+        if t < self._times[0]:
+            return self._values[0]
+        if t > self._times[-1]:
+            return self._values[-1]
+        i = bisect_right(self._times, t) - 1
+        if i == len(self._times) - 1:
+            return self._values[i]
+        return self._values[i] + (t - self._times[i]) * (
+            self._values[i + 1] - self._values[i]
+        ) / (self._times[i + 1] - self._times[i])
+
+    def times(self):
+        """
+        Returns a list of the times in this protocol.
+        """
+        return self._times
+
+    def values(self):
+        """
+        Returns a list of the values in this protocol.
+        """
+        return self._values
diff --git a/myokit/_sim/cable.c b/myokit/_sim/cable.c
index 88059769c..474957f1c 100644
--- a/myokit/_sim/cable.c
+++ b/myokit/_sim/cable.c
@@ -26,11 +26,11 @@ model.create_unique_names()
 w = ansic.AnsiCExpressionWriter()
 
 # Process bindings, remove unsupported bindings.
-bound_variables = model.prepare_bindings({
+bound_variables = myokit._prepare_bindings(model, {
     'time'         : 'engine_time',
     'pace'         : 'engine_pace',
     'diffusion_current' : 'diffusion_current',
-    })
+})
 
 # Define var/lhs function
 def v(var, pre='cell->'):
diff --git a/myokit/_sim/cmodel.h b/myokit/_sim/cmodel.h
index 5481f12c7..df92f691d 100644
--- a/myokit/_sim/cmodel.h
+++ b/myokit/_sim/cmodel.h
@@ -45,7 +45,8 @@ Model variables are divided into several (non-overlapping) groups:
 State derivatives:
     Calculated by the model.
 Bound variables:
-    External inputs to the model (e.g. time and pacing).
+    External inputs to the model (e.g. time and pacing). There are an arbitrary
+    number of these, and the model stores the current value of each
 Intermediary variables:
     The remaining variables that depend on state variables.
 Constants:
@@ -214,6 +215,8 @@ typedef int Model_Flag;
 /* Logging sensitivities */
 #define Model_NO_SENSITIVITIES_TO_LOG       -300
 #define Model_SENSITIVITY_LOG_APPEND_FAILED -303
+/* Pacing */
+#define Model_INVALID_PACING                -400
 
 /* Caching doesn't help much when running without jacobians etc., so disabled
    for now
@@ -259,7 +262,9 @@ Model_SetPyErr(Model_Flag flag)
     case Model_SENSITIVITY_LOG_APPEND_FAILED:
         PyErr_SetString(PyExc_Exception, "CModel error: Call to append() failed on sensitivity matrix logging list.");
         break;
-
+    case Model_INVALID_PACING:
+        PyErr_SetString(PyExc_Exception, "CModel error: Invalid pacing provided.");
+        break;
     /* Unknown */
     default:
         PyErr_Format(PyExc_Exception, "CModel error: Unlisted error %d", (int)flag);
@@ -277,9 +282,12 @@ struct Model_Memory {
     /* If this model has sensitivities this will be 1, otherwise 0. */
     int has_sensitivities;
 
+    /* Pacing */
+    realtype *pace_values;
+    int n_pace;
+
     /* Bound variables */
     realtype time;
-    realtype pace;
     realtype realtime;
     realtype evaluations;
 
@@ -444,6 +452,41 @@ Model_ClearCache(Model model)
     return Model_OK;
 }
 
+/*
+ * Sets up (i.e. allocates memory for) array of protocol-determined values
+ *
+ * Arguments
+ *  n_pace: the number of pacing values to use.
+ *
+ * Returns a model flag.
+ *
+ */
+Model_Flag
+Model_SetupPacing(Model model, int n_pace)
+{
+    if (model == NULL) return Model_INVALID_MODEL;
+    if (n_pace < 0) return Model_INVALID_PACING;
+
+    /* Free any existing pacing */
+    if (model->n_pace > 0) {
+        free(model->pace_values);
+    }
+
+    /* Allocate new pacing */
+    model->n_pace = n_pace;
+    model->pace_values = (realtype*)malloc(n_pace * sizeof(realtype));
+    if (model->pace_values == NULL) {
+        return Model_OUT_OF_MEMORY;
+    }
+
+    /* Clear values */
+    for (int i = 0; i < n_pace; i++) {
+        model->pace_values[i] = 0;
+    }
+
+    return Model_OK;
+}
+
 /*
  * (Re)calculates the values of all constants that are derived from other
  * constants.
@@ -503,14 +546,14 @@ for eq in parameter_derived.values():
  *  literals : An array of size model->n_literals
  *
  * Returns a model flag.
- */
+ *
 Model_Flag
 Model_SetLiteralVariables(Model model, const realtype* literals)
 {
     int i;
     if (model == NULL) return Model_INVALID_MODEL;
 
-    /* Scan for changes */
+    /* Scan for changes *
     i = 0;
     #ifdef Model_CACHING
     if (Model__ValidCache(model)) {
@@ -522,7 +565,7 @@ Model_SetLiteralVariables(Model model, const realtype* literals)
     }
     #endif
 
-    /* Update remaining */
+    /* Update remaining *
     if (i < model->n_literals) {
         for (; i<model->n_literals; i++) {
             model->literals[i] = literals[i];
@@ -535,7 +578,7 @@ Model_SetLiteralVariables(Model model, const realtype* literals)
     }
 
     return Model_OK;
-}
+}*/
 
 /*
  * Updates the parameter variables to the values given in `parameters`.
@@ -657,8 +700,10 @@ Model_SetParametersFromIndependents(Model model, const realtype* independents)
 Model_Flag
 Model_SetBoundVariables(
     Model model,
-    const realtype time, const realtype pace,
-    const realtype realtime, const realtype evaluations)
+    const realtype time,
+    const realtype *pace_values,
+    const realtype realtime,
+    const realtype evaluations)
 {
     int changed;
     if (model == NULL) return Model_INVALID_MODEL;
@@ -669,13 +714,13 @@ Model_SetBoundVariables(
         changed = 1;
     }
 
-<?
-if model.binding('pace') is not None:
-    print(tab + 'if (pace != model->pace) {')
-    print(tab + '    model->pace = pace;')
-    print(tab + '    changed = 1;')
-    print(tab + '}')
-?>
+    for (int i = 0; i < model->n_pace; i++) {
+        if (pace_values[i] != model->pace_values[i]) {
+            model->pace_values[i] = pace_values[i];
+            changed = 1;
+        }
+    }
+
     #ifdef Model_CACHING
     if (changed) {
         Model__InvalidateCache(model);
@@ -1132,6 +1177,10 @@ Model Model_Create(Model_Flag* flagp)
     model->parameters = (realtype*)malloc(model->n_parameters * sizeof(realtype));
     model->parameter_derived = (realtype*)malloc(model->n_parameter_derived * sizeof(realtype));
 
+    /* Pacing */
+    model->n_pace = 0;
+    model->pace_values = NULL;
+
     /* Literals */
     model->n_literals = <?= len(literals) ?>;
     model->n_literal_derived = <?= len(literal_derived) ?>;
@@ -1194,9 +1243,8 @@ for i, expr in enumerate(s_independents):
      * Default values
      */
 
-    /* Bound variables */
+    /* Bound variables (except pacing) */
     model->time = 0;
-    model->pace = 0;
     model->realtime = 0;
     model->evaluations = 0;
 
diff --git a/myokit/_sim/cmodel.py b/myokit/_sim/cmodel.py
index 30e39926b..cec9829b6 100644
--- a/myokit/_sim/cmodel.py
+++ b/myokit/_sim/cmodel.py
@@ -35,6 +35,9 @@ class CModel(object):
 
     ``model``
         A :class:`myokit.Model`.
+    ``pacing_labels``
+        A list of ``str`` variable labels, each corresponding to a paced
+        variable.
     ``sensitivities``
         Either ``None`` or a tuple ``(dependents, independents)``. See
         :class:`myokit.Simulation` for details.
@@ -76,7 +79,7 @@ class CModel(object):
         Constants that depend on literals, but not on parameters.
 
     """
-    def __init__(self, model, sensitivities):
+    def __init__(self, model, pacing_labels, sensitivities):
 
         # Parse sensitivity arguments
         has_sensitivities, dependents, independents = \
@@ -86,17 +89,15 @@ def __init__(self, model, sensitivities):
         # bother with keywords.
         model.create_unique_names()
 
-        # Remove any unused bindings, and get mapping from variables to C
-        # variable names as used in model.h
-        #TODO: Think about best way to do this for model re-use with different
-        # sets of bound variables... Presumably the model would simply support
-        # all bindings used by any of the simulations based on model.h ?
-        bound_variables = model.prepare_bindings({
+        # Get mapping from variables to C variable names as used in model.h
+        labels = {
             'time': 'time',
-            'pace': 'pace',
             'realtime': 'realtime',
             'evaluations': 'evaluations',
-        })
+        }
+        for i, label in enumerate(pacing_labels):
+            labels[label] = 'pace_values[' + str(i) + ']'
+        bound_variables = myokit._prepare_bindings(model, labels)
 
         # Get equations in solvable order (grouped by component)
         equations = model.solvable_order()
diff --git a/myokit/_sim/cvodessim.c b/myokit/_sim/cvodessim.c
index 0ded674a6..33683fadb 100644
--- a/myokit/_sim/cvodessim.c
+++ b/myokit/_sim/cvodessim.c
@@ -195,11 +195,19 @@ Model model;        /* A model object */
 /*
  * Pacing
  */
-ESys epacing;           /* Event-based pacing system */
-PyObject* eprotocol;    /* An event-based pacing protocol */
-FSys fpacing;           /* Fixed-form pacing system */
-PyObject* fprotocol;    /* A fixed-form pacing protocol */
-double pace = 0;        /* Pacing value */
+union PSys {
+    ESys event;
+    FSys fixed;
+};
+enum PSysType {
+    EVENT,
+    FIXED
+};
+union PSys *pacing_systems;   /* Array of pacing system (event or fixed) */
+enum PSysType *pacing_types;  /* Array of pacing system types */
+PyObject *protocols;          /* The protocols used to generate the pacing systems */
+double* pacing;               /* Pacing values, same size as pacing_systems and pacing_types */
+int n_pace;                   /* The number of pacing systems */
 
 /*
  * CVODE Memory
@@ -350,12 +358,14 @@ rhs(realtype t, N_Vector y, N_Vector ydot, void *user_data)
     UserData fdata;
     int i;
 
-    /* Fixed-form pacing? Then look-up correct value of pacing variable! */
-    if (fpacing != NULL) {
-        pace = FSys_GetLevel(fpacing, t, &flag_fpacing);
-        if (flag_fpacing != FSys_OK) { /* This should never happen */
-            FSys_SetPyErr(flag_fpacing);
-            return -1;  /* Negative value signals irrecoverable error to CVODE */
+    /* Fixed-form pacing? Then look-up correct value of pacing variable */
+    for (int i = 0; i < n_pace; i++) {
+        if (pacing_types[i] == FIXED) {
+            pacing[i] = FSys_GetLevel(pacing_systems[i].fixed, t, &flag_fpacing);
+            if (flag_fpacing != FSys_OK) { /* This should never happen */
+                FSys_SetPyErr(flag_fpacing);
+                return -1;  /* Negative value signals irrecoverable error to CVODE */
+            }
         }
     }
 
@@ -363,7 +373,7 @@ rhs(realtype t, N_Vector y, N_Vector ydot, void *user_data)
 
     /* Set time, pace, evaluations and realtime */
     evaluations++;
-    Model_SetBoundVariables(model, t, pace, realtime, evaluations);
+    Model_SetBoundVariables(model, t, pacing, realtime, evaluations);
 
     /* Set sensitivity parameters */
     if (model->has_sensitivities) {
@@ -463,8 +473,16 @@ sim_clean()
         }
 
         /* Pacing systems */
-        ESys_Destroy(epacing); epacing = NULL;
-        FSys_Destroy(fpacing); fpacing = NULL;
+        for (int i = 0; i < n_pace; i++) {
+            if (pacing_types[i] == FIXED) {
+                FSys_Destroy(pacing_systems[i].fixed);
+            } else if (pacing_types[i] == EVENT) {
+                ESys_Destroy(pacing_systems[i].event);
+            }
+        }
+        free(pacing_systems); pacing_systems = NULL;
+        free(pacing_types); pacing_types = NULL;
+        free(pacing); pacing = NULL;
 
         /* CModel */
         Model_Destroy(model); model = NULL;
@@ -557,8 +575,9 @@ sim_init(PyObject *self, PyObject *args)
     /* Set all global pointers to null */
     /* Model and pacing */
     model = NULL;
-    epacing = NULL;
-    fpacing = NULL;
+    pacing_types = NULL;
+    pacing_systems = NULL;
+    pacing = NULL;
     /* User data and parameter scaling */
     udata = NULL;
     pbar = NULL;
@@ -586,8 +605,8 @@ sim_init(PyObject *self, PyObject *args)
     sundials_context = NULL;
     #endif
 
-    /* Check input arguments     0123456789012345678 */
-    if (!PyArg_ParseTuple(args, "ddOOOOOOOOdOOidOOi",
+    /* Check input arguments     01234567890123456 */
+    if (!PyArg_ParseTuple(args, "ddOOOOOOOdOOidOOi",
             &tmin,              /*  0. Float: initial time */
             &tmax,              /*  1. Float: final time */
             &state_py,          /*  2. List: initial and final state */
@@ -595,17 +614,16 @@ sim_init(PyObject *self, PyObject *args)
             &bound_py,          /*  4. List: store final bound variables here */
             &literals,          /*  5. List: literal constant values */
             &parameters,        /*  6. List: parameter values */
-            &eprotocol,         /*  7. Event-based protocol */
-            &fprotocol,         /*  8. Fixed-form protocol (tuple) */
-            &log_dict,          /*  9. DataLog */
-            &log_interval,      /* 10. Float: log interval, or 0 */
-            &log_times,         /* 11. List of logging times, or None */
-            &sens_list,         /* 12. List to store sensitivities in */
-            &rf_indice,         /* 13. Int: root-finding state variable */
-            &rf_threshold,      /* 14. Float: root-finding threshold */
-            &rf_list,           /* 15. List to store roots in or None */
-            &benchmarker,       /* 16. myokit.tools.Benchmarker object */
-            &log_realtime       /* 17. Int: 1 if logging real time */
+            &protocols,        /*   7. Event-based or fixed protocols */
+            &log_dict,          /*  8. DataLog */
+            &log_interval,      /*  9. Float: log interval, or 0 */
+            &log_times,         /* 10. List of logging times, or None */
+            &sens_list,         /* 11. List to store sensitivities in */
+            &rf_indice,         /* 12. Int: root-finding state variable */
+            &rf_threshold,      /* 13. Float: root-finding threshold */
+            &rf_list,           /* 14. List to store roots in or None */
+            &benchmarker,       /* 15. myokit.tools.Benchmarker object */
+            &log_realtime       /* 16. Int: 1 if logging real time */
     )) {
         PyErr_SetString(PyExc_Exception, "Incorrect input arguments.");
         return 0;
@@ -929,48 +947,74 @@ sim_init(PyObject *self, PyObject *args)
     }
 
     /*
-     * Set up pacing system
+     * Set up pacing systems
      */
-
-    /* Set up event-based pacing */
-    if (eprotocol != Py_None) {
-        epacing = ESys_Create(&flag_epacing);
-        if (flag_epacing != ESys_OK) { ESys_SetPyErr(flag_epacing); return sim_clean(); }
-        flag_epacing = ESys_Populate(epacing, eprotocol);
-        if (flag_epacing != ESys_OK) { ESys_SetPyErr(flag_epacing); return sim_clean(); }
-        flag_epacing = ESys_AdvanceTime(epacing, tmin);
-        if (flag_epacing != ESys_OK) { ESys_SetPyErr(flag_epacing); return sim_clean(); }
-        tnext = ESys_GetNextTime(epacing, &flag_epacing);
-        pace = ESys_GetLevel(epacing, &flag_epacing);
-        tnext = (tnext < tmax) ? tnext : tmax;
-
-        #ifdef MYOKIT_DEBUG_PROFILING
-        benchmarker_print("CP Created event-based pacing system.");
-        #endif
-    } else {
-        tnext = tmax;
+    n_pace = 0;
+    if (protocols != Py_None) {
+        if (!PyList_Check(protocols)) {
+            return sim_cleanx(PyExc_TypeError, "'protocols' must be a list.");
+        }
+        n_pace = PyList_Size(protocols);
+    }
+    pacing_systems = (union PSys*)malloc(sizeof(union PSys) * n_pace);
+    if (pacing_systems == NULL) {
+        return sim_cleanx(PyExc_Exception, "Unable to allocate space to store pacing systems.");
+    }
+    pacing_types = (enum PSysType *)malloc(sizeof(enum PSysType) * n_pace);
+    if (pacing_types == NULL) {
+        return sim_cleanx(PyExc_Exception, "Unable to allocate space to store pacing types.");
     }
+    pacing = (realtype*)malloc(sizeof(realtype) * n_pace);
+    if (pacing == NULL) {
+        return sim_cleanx(PyExc_Exception, "Unable to allocate space to store pacing values.");
+    }
+    Model_SetupPacing(model, n_pace);
 
-    /* Set up fixed-form pacing */
-    if (eprotocol == Py_None && fprotocol != Py_None) {
-        /* Check 'protocol' is tuple (times, values) */
-        if (!PyTuple_Check(fprotocol)) {
-            return sim_cleanx(PyExc_TypeError, "Fixed-form pacing protocol should be tuple or None.");
-        }
-        if (PyTuple_Size(fprotocol) != 2) {
-            return sim_cleanx(PyExc_ValueError, "Fixed-form pacing protocol tuple should have size 2.");
-        }
-        /* Create fixed-form pacing object and populate */
-        fpacing = FSys_Create(&flag_fpacing);
-        if (flag_fpacing != FSys_OK) { FSys_SetPyErr(flag_fpacing); return sim_clean(); }
-        flag_fpacing = FSys_Populate(fpacing,
-            PyTuple_GetItem(fprotocol, 0),  /* Borrowed, no decref */
-            PyTuple_GetItem(fprotocol, 1));
-        if (flag_fpacing != FSys_OK) { FSys_SetPyErr(flag_fpacing); return sim_clean(); }
+    /*
+     *  Unless set by pacing, tnext is set to tmax
+     */
+    tnext = tmax;
 
-        #ifdef MYOKIT_DEBUG_PROFILING
-        benchmarker_print("CP Created fixed-form pacing system.");
-        #endif
+    /* 
+     * Set up event-based and/or fixed pacing.
+     */
+    if (protocols != Py_None) {
+        for (int i = 0; i < PyList_Size(protocols); i++) {
+            PyObject *protocol = PyList_GetItem(protocols, i);
+            const char* protocol_type_name = Py_TYPE(protocol)->tp_name;
+            if (strcmp(protocol_type_name, "Protocol") == 0) {
+                pacing_systems[i].event = ESys_Create(&flag_epacing);
+                pacing_types[i] = EVENT;
+                ESys epacing = pacing_systems[i].event;
+                if (flag_epacing != ESys_OK) { ESys_SetPyErr(flag_epacing); return sim_clean(); }
+                flag_epacing = ESys_Populate(epacing, protocol);
+                if (flag_epacing != ESys_OK) { ESys_SetPyErr(flag_epacing); return sim_clean(); }
+                flag_epacing = ESys_AdvanceTime(epacing, tmin);
+                if (flag_epacing != ESys_OK) { ESys_SetPyErr(flag_epacing); return sim_clean(); }
+                const double ptnext = ESys_GetNextTime(epacing, &flag_epacing);
+                pacing[i] = ESys_GetLevel(epacing, &flag_epacing);
+                tnext = fmin(ptnext, tnext);
+
+                #ifdef MYOKIT_DEBUG_PROFILING
+                benchmarker_print("CP Created event-based pacing system.");
+                #endif
+            } else if (strcmp(protocol_type_name, "TimeSeriesProtocol") == 0) {
+                pacing_systems[i].fixed = FSys_Create(&flag_fpacing);
+                pacing_types[i] = FIXED;
+                FSys fpacing = pacing_systems[i].fixed;
+                if (flag_fpacing != FSys_OK) { FSys_SetPyErr(flag_fpacing); return sim_clean(); }
+                flag_fpacing = FSys_Populate(fpacing, protocol);
+                if (flag_fpacing != FSys_OK) { FSys_SetPyErr(flag_fpacing); return sim_clean(); }
+
+
+                #ifdef MYOKIT_DEBUG_PROFILING
+                benchmarker_print("CP Created fixed-form pacing system.");
+                #endif
+            } else {
+                printf("protocol_type_name: %s", protocol_type_name);
+                return sim_cleanx(PyExc_TypeError, "Item %d in 'protocols' is not a myokit.Protocol or myokit.TimeSeriesProtocol object.", i);
+            }
+        }
     }
 
     /*
@@ -1242,8 +1286,8 @@ sim_step(PyObject *self, PyObject *args)
     /* Number of integration steps taken in this call */
     int steps_taken = 0;
 
-    /* Proposed next logging point */
-    double proposed_tlog;
+    /* Proposed next logging or pacing point */
+    double t_proposed;
 
     /* Multi-purpose Python objects */
     PyObject *val;
@@ -1288,9 +1332,11 @@ sim_step(PyObject *self, PyObject *args)
                     /* PyList_SetItem steals a reference: no need to decref the double! */
                 }
                 PyList_SetItem(bound_py, 0, PyFloat_FromDouble(tlast));
-                PyList_SetItem(bound_py, 1, PyFloat_FromDouble(pace));
-                PyList_SetItem(bound_py, 2, PyFloat_FromDouble(realtime));
-                PyList_SetItem(bound_py, 3, PyFloat_FromDouble(evaluations));
+                PyList_SetItem(bound_py, 1, PyFloat_FromDouble(realtime));
+                PyList_SetItem(bound_py, 2, PyFloat_FromDouble(evaluations));
+                for (int i = 0; i < n_pace; i++) {
+                    PyList_SetItem(bound_py, 3 + i, PyFloat_FromDouble(pacing[i]));
+                }
                 return sim_clean();
             }
 
@@ -1435,7 +1481,7 @@ sim_step(PyObject *self, PyObject *args)
                         if (ilog < PySequence_Size(log_times)) {
                             val = PySequence_GetItem(log_times, ilog); /* New reference */
                             if (PyFloat_Check(val)) {
-                                proposed_tlog = PyFloat_AsDouble(val);
+                                t_proposed = PyFloat_AsDouble(val);
                                 Py_DECREF(val);
                             } else if (PyNumber_Check(val)) {
                                 ret = PyNumber_Float(val);  /* New reference */
@@ -1443,17 +1489,17 @@ sim_step(PyObject *self, PyObject *args)
                                 if (ret == NULL) {
                                     return sim_cleanx(PyExc_ValueError, "Unable to cast entry in 'log_times' to float.");
                                 } else {
-                                    proposed_tlog = PyFloat_AsDouble(ret);
+                                    t_proposed = PyFloat_AsDouble(ret);
                                     Py_DECREF(ret);
                                 }
                             } else {
                                 Py_DECREF(val);
                                 return sim_cleanx(PyExc_ValueError, "Entries in 'log_times' must be floats.");
                             }
-                            if (proposed_tlog < tlog) {
+                            if (t_proposed < tlog) {
                                 return sim_cleanx(PyExc_ValueError, "Values in log_times must be non-decreasing.");
                             }
-                            tlog = proposed_tlog;
+                            tlog = t_proposed;
                             ilog++;
                             val = NULL;
                         } else {
@@ -1469,14 +1515,18 @@ sim_step(PyObject *self, PyObject *args)
              * At this point we have logged everything _before_ time t, so it
              * is safe to update the pacing mechanism to time t.
              */
-            if (epacing != NULL) {
-                flag_epacing = ESys_AdvanceTime(epacing, t);
-                if (flag_epacing != ESys_OK) {
-                    ESys_SetPyErr(flag_epacing); return sim_clean();
+            tnext = tmax;
+            for (int i = 0; i < n_pace; i++) {
+                if (pacing_types[i] == EVENT) {
+                    ESys epacing = pacing_systems[i].event;
+                    flag_epacing = ESys_AdvanceTime(epacing, t);
+                    if (flag_epacing != ESys_OK) {
+                        ESys_SetPyErr(flag_epacing); return sim_clean();
+                    }
+                    t_proposed = ESys_GetNextTime(epacing, NULL);
+                    tnext = fmin(tnext, t_proposed);
+                    pacing[i] = ESys_GetLevel(epacing, NULL);
                 }
-                tnext = ESys_GetNextTime(epacing, NULL);
-                tnext = (tnext < tmax) ? tnext : tmax;
-                pace = ESys_GetLevel(epacing, NULL);
             }
 
             /* Dynamic logging: Log every visited point */
@@ -1500,7 +1550,7 @@ sim_step(PyObject *self, PyObject *args)
                 } else if (model->logging_bound) {
                     /* Logging bounds but not derivs or inters: No need to run
                        full rhs, just update bound variables */
-                    Model_SetBoundVariables(model, t, pace, realtime, evaluations);
+                    Model_SetBoundVariables(model, t, pacing, realtime, evaluations);
                 }
 
                 /* Write to log */
@@ -1584,9 +1634,11 @@ sim_step(PyObject *self, PyObject *args)
 
     /* Set bound variable values */
     PyList_SetItem(bound_py, 0, PyFloat_FromDouble(t));
-    PyList_SetItem(bound_py, 1, PyFloat_FromDouble(pace));
-    PyList_SetItem(bound_py, 2, PyFloat_FromDouble(realtime));
-    PyList_SetItem(bound_py, 3, PyFloat_FromDouble(evaluations));
+    PyList_SetItem(bound_py, 1, PyFloat_FromDouble(realtime));
+    PyList_SetItem(bound_py, 2, PyFloat_FromDouble(evaluations));
+    for (int i = 0; i < n_pace; i++) {
+        PyList_SetItem(bound_py, 3 + i, PyFloat_FromDouble(pacing[i]));
+    }
 
     #ifdef MYOKIT_DEBUG_PROFILING
     benchmarker_print("CP Set final state and bound variable values.");
@@ -1606,7 +1658,8 @@ sim_eval_derivatives(PyObject *self, PyObject *args)
     int i;
     int success;
     double time_in;
-    double pace_in;
+    PyObject *pace_in;
+    double *pacing_in;
     Model model;
     Model_Flag flag_model;
     PyObject *state;
@@ -1620,9 +1673,9 @@ sim_eval_derivatives(PyObject *self, PyObject *args)
 
     /* Check input arguments */
     /* Check input arguments     0123456789ABCDEF*/
-    if (!PyArg_ParseTuple(args, "ddOOOO",
+    if (!PyArg_ParseTuple(args, "dOOOOO",
             &time_in,           /* 0. Float: time */
-            &pace_in,           /* 1. Float: pace */
+            &pace_in,           /* 1. List: pace */
             &state,             /* 2. List: state */
             &deriv,             /* 3. List: store derivatives here */
             &literals,          /* 4. List: literal constant values */
@@ -1634,6 +1687,10 @@ sim_eval_derivatives(PyObject *self, PyObject *args)
     }
 
     /* Check lists are sequences */
+    if (!PyList_Check(pace_in)) {
+        PyErr_SetString(PyExc_Exception, "Pace argument must be a list.");
+        return 0;
+    }
     if (!PyList_Check(state)) {
         PyErr_SetString(PyExc_Exception, "State argument must be a list.");
         return 0;
@@ -1665,8 +1722,25 @@ sim_eval_derivatives(PyObject *self, PyObject *args)
         goto error;
     }
 
+    flag_model = Model_SetupPacing(model, n_pace);
+    if (flag_model != Model_OK) {
+        Model_SetPyErr(flag_model);
+        goto error;
+    }
+
+    /* Set pacing values */
+    pacing_in = (double*)malloc(sizeof(double) * n_pace);
+    for (int i = 0; i < n_pace; i++) {
+        val = PyList_GetItem(pace_in, i); /* Don't decref */
+        if (!PyFloat_Check(val)) {
+            PyErr_Format(PyExc_Exception, "Item %d in pace vector is not a float.", i);
+            goto error;
+        }
+        pacing_in[i] = PyFloat_AsDouble(val);
+    }
+
     /* Set bound variables */
-    Model_SetBoundVariables(model, time_in, pace_in, 0, 0);
+    Model_SetBoundVariables(model, time_in, pacing_in, 0, 0);
 
     /* Set literal values */
     for (i=0; i<model->n_literals; i++) {
diff --git a/myokit/_sim/cvodessim.py b/myokit/_sim/cvodessim.py
index 9abbbf2bc..c79748d39 100644
--- a/myokit/_sim/cvodessim.py
+++ b/myokit/_sim/cvodessim.py
@@ -122,8 +122,11 @@ class Simulation(myokit.CModule):
     ``model``
         The model to simulate
     ``protocol``
-        An optional :class:`myokit.Protocol` to use as input for variables
-        bound to ``pace``.
+        A :class:`myokit.Protocol` or :class:`myokit.TimeSeriesProtocol` to use
+        for the variable with binding ``pace``. Atlernatively, a dictionary
+        mapping binding labels to :class:`myokit.Protocol` objects can be used
+        to run with multiple protocols. Finally, can be ``None`` to run without
+        a protocol.
     ``sensitivities``
         An optional tuple ``(dependents, independents)`` where ``dependents``
         is a list of variables or expressions to take derivatives of (``y`` in
@@ -158,13 +161,20 @@ def __init__(self, model, protocol=None, sensitivities=None, path=None):
         del model
 
         # Set protocol
-        self._protocol = None
-        self._fixed_form_protocol = None
-        self.set_protocol(protocol)
-        del protocol
+        self._protocols = []
+        self._pacing_labels = []
+        if isinstance(protocol, (myokit.Protocol, myokit.TimeSeriesProtocol)):
+            protocol = {'pace': protocol}
+        elif protocol is None:
+            # TODO: This can be an empty dict once #320 is resolved
+            protocol = {'pace': None}
+        for label, protocol in protocol.items():
+            self._pacing_labels.append(label)
+            self._protocols.append(myokit.Protocol())
+            self.set_protocol(protocol, label)
 
         # Generate C Model code, get sensitivity and constants info
-        cmodel = myokit.CModel(self._model, sensitivities)
+        cmodel = myokit.CModel(self._model, self._pacing_labels, sensitivities)
         if cmodel.has_sensitivities:
             self._sensitivities = (cmodel.dependents, cmodel.independents)
 
@@ -299,7 +309,8 @@ def _store_build(self, path, d_build, name):
             with open(fname, 'wb') as f:
                 pickle.dump(str(name), f)
                 pickle.dump(self._model, f)
-                pickle.dump(self._protocol, f)
+                pickle.dump(self._protocols, f)
+                pickle.dump(self._pacing_labels, f)
                 pickle.dump(sens_arg, f)
 
             # Zip it all in
@@ -345,7 +356,8 @@ def from_path(path):
             with open(fname, 'rb') as f:
                 name = pickle.load(f)
                 model = pickle.load(f)
-                protocol = pickle.load(f)
+                protocols = pickle.load(f)
+                pacing_labels = pickle.load(f)
                 sensitivities = pickle.load(f)
 
             # Load module
@@ -353,7 +365,12 @@ def from_path(path):
             module = load_module(name, d_build)
 
             # Create and return simulation
-            return Simulation(model, protocol, sensitivities, (path, module))
+            labeled_protocols = {
+                k: v for k, v in zip(pacing_labels, protocols)
+            }
+            return Simulation(
+                model, labeled_protocols, sensitivities, (path, module)
+            )
 
         finally:
             myokit.tools.rmtree(d_build, silent=True)
@@ -384,12 +401,16 @@ def last_state(self):
         """
         return list(self._error_state) if self._error_state else None
 
-    def eval_derivatives(self, y=None):
+    def eval_derivatives(self, y=None, pacing=None):
         """
         Evaluates and returns the state derivatives.
 
         The state to evaluate for can be given as ``y``. If no state is given
         the current simulation state is used.
+
+        An optional dict ``pacing`` can be passed in that maps labels to
+        numerical values. If a label is used but not provided in ``pacing``,
+        it's value will be set to 0.
         """
         # Get state
         if y is None:
@@ -397,6 +418,15 @@ def eval_derivatives(self, y=None):
         else:
             y = self._model.map_to_state(y)
 
+        pacing_values = [0.0] * len(self._pacing_labels)
+        if pacing is not None:
+            for k, v in pacing.items():
+                try:
+                    ki = self._pacing_labels.index(k)
+                    pacing_values[ki] = float(v)
+                except ValueError:
+                    pass
+
         # Create space to store derivatives
         dy = list(self._state)
 
@@ -405,7 +435,7 @@ def eval_derivatives(self, y=None):
             # 0. Time
             0,
             # 1. Pace
-            0,
+            pacing_values,
             # 2. State
             y,
             # 3. Space to store the state derivatives
@@ -496,16 +526,19 @@ def __reduce__(self):
             # etc.
             sens_arg = [[x.code() for x in y] for y in self._sensitivities]
 
+        protocols = {
+            k: p for k, p in zip(self._pacing_labels, self._protocols)
+        }
+
         return (
             self.__class__,
-            (self._model, self._protocol, sens_arg),
+            (self._model, protocols, sens_arg),
             (
                 self._time,
                 self._state,
                 self._default_state,
                 self._s_state,
                 self._s_default_state,
-                self._fixed_form_protocol,  # Can't be set in constructor
                 self._tolerance,
                 self._dtmin,
                 self._dtmax,
@@ -733,7 +766,7 @@ def _run(self, duration, log, log_interval, log_times, sensitivities,
                 s_state = [list(x) for x in self._s_state]
 
             # List to store final bound variables in (for debugging)
-            bound = [0, 0, 0, 0]
+            bound = [0, 0, 0] + [0] * len(self._pacing_labels)
 
             # Initialize
             if myokit.DEBUG_SP:
@@ -753,30 +786,28 @@ def _run(self, duration, log, log_interval, log_times, sensitivities,
                 list(self._literals.values()),
                 # 6. Parameter values
                 list(self._parameters.values()),
-                # 7. An event-based pacing protocol
-                self._protocol,
-                # 8. A fixed-form protocol
-                self._fixed_form_protocol,
-                # 9. A DataLog
+                # 7. Pacing protocols
+                self._protocols,
+                # 8. A DataLog
                 log,
-                # 10. The log interval, or 0
+                # 9. The log interval, or 0
                 log_interval,
-                # 11. A list of predetermind logging times, or None
+                # 10. A list of predetermind logging times, or None
                 log_times,
-                # 12. A list to store calculated sensitivities in
+                # 11. A list to store calculated sensitivities in
                 sensitivities,
-                # 13. The state variable indice for root finding (only used if
+                # 12. The state variable indice for root finding (only used if
                 #     root_list is a list)
                 root_indice,
-                # 14. The threshold for root crossing (can be 0 too, only used
+                # 13. The threshold for root crossing (can be 0 too, only used
                 #     if root_list is a list).
                 root_threshold,
-                # 15. A list to store calculated root crossing times and
+                # 14. A list to store calculated root crossing times and
                 #     directions in, or None
                 root_list,
-                # 16. A myokit.tools.Benchmarker or None (if not used)
+                # 15. A myokit.tools.Benchmarker or None (if not used)
                 b,
-                # 17. Boolean/int: 1 if we are logging realtime
+                # 16. Boolean/int: 1 if we are logging realtime
                 int(self._model.binding('realtime') is not None),
             )
             t = tmin
@@ -812,9 +843,12 @@ def _run(self, duration, log, log_interval, log_times, sensitivities,
                             in self._model.format_state(state).splitlines()])
                 txt.append('Inputs for binding:')
                 txt.append('  time        = ' + myokit.float.str(bound[0]))
-                txt.append('  pace        = ' + myokit.float.str(bound[1]))
-                txt.append('  realtime    = ' + myokit.float.str(bound[2]))
-                txt.append('  evaluations = ' + myokit.float.str(bound[3]))
+                txt.append('  realtime    = ' + myokit.float.str(bound[1]))
+                txt.append('  evaluations = ' + myokit.float.str(bound[2]))
+                for i, label in enumerate(self._pacing_labels):
+                    txt.append(
+                        '  ' + label + ' = ' + myokit.float.str(bound[3 + i])
+                    )
                 txt.append(str(e))
 
                 # Check if state derivatives can be evaluated in Python, if
@@ -951,66 +985,51 @@ def set_min_step_size(self, dtmin=None):
 
     def set_fixed_form_protocol(self, times=None, values=None):
         """
-        Configures this simulation to run with a predetermined protocol
-        instead of the usual event-based mechanism.
-
-        A 1D time-series should be given as input. During the simulation, the
-        value of the pacing variable will be determined by linearly
-        interpolating between the two nearest points in the series. If the
-        simulation time is outside the bounds of the time-series, the first or
-        last value in the series will be used.
-
-        Setting a predetermined protocol clears any previously set (event-based
-        or pre-determined) protocol. To clear all protocols, call this method
-        with `times=None`. When a simulation is run without any protocol, the
-        value of any variables bound to `pace` will be set to 0.
-
-        Arguments:
-
-        ``times``
-            A non-decreasing array of times. If any times appear more than
-            once, only the value at the highest index will be used.
-        ``values``
-            An array of values for the pacing variable. Must have the same size
-            as ``times``.
+        Sets a :class:`TimeSeriesProtocol` specified by ``times`` and
+        ``values`` for the label ``pace``.
 
+        This method is provided for backwards compatibility with older
+        versions, please use :meth:`set_protocol` and the
+        :class:`TimeSeriesProtocol` class instead.
         """
-        # Check input
-        if times is None:
-            if values is not None:
-                raise ValueError('Values array given, but no times array.')
-        else:
-            if values is None:
-                raise ValueError('Times array given, but no values array.')
-            if len(times) != len(values):
-                raise ValueError('Times and values array must have same size.')
-
-        # Clear event-based protocol, if set
-        self._protocol = None
+        # Deprecated on 2023-06-02
+        import warnings
+        warnings.warn(
+            'The method `myokit.Simulation.set_fixed_form_protocol` is '
+            'deprecated. It will be removed in future versions of Myokit.'
+        )
 
-        # Set new protocol
+        if times is None and values is None:
+            self.set_protocol(None)
+            return
         if times is None:
-            # Clear predetermined protocol
-            self._fixed_form_protocol = None
-        else:
-            # Copy data and set
-            self._fixed_form_protocol = (list(times), list(values))
+            raise ValueError('No times given.')
+        if values is None:
+            raise ValueError('No values given.')
+
+        self.set_protocol(myokit.TimeSeriesProtocol(times, values))
 
-    def set_protocol(self, protocol=None):
+    def set_protocol(self, protocol, label='pace'):
         """
-        Sets the pacing :class:`Protocol` used by this simulation.
+        Set an event-based pacing :class:`Protocol` or a :class:`FixedProtocol`
+        for the given ``label``.
 
-        To run without pacing call this method with ``protocol = None``. In
-        this case, the value of any variables bound to `pace` will be set to 0.
+        To remove a previously set binding call this method with ``protocol =
+        None``. In this case, the value of any variables bound to ``label``
+        will be set to 0.
+
+        The label must be one of the pacing labels set in the constructor.
         """
-        # Clear predetermined protocol, if set
-        self._fixed_form_protocol = None
+        try:
+            index = self._pacing_labels.index(label)
+        except ValueError:
+            raise ValueError('Unknown pacing label: ' + str(label))
 
         # Set new protocol
         if protocol is None:
-            self._protocol = None
+            self._protocols[index] = myokit.Protocol()
         else:
-            self._protocol = protocol.clone()
+            self._protocols[index] = protocol.clone()
 
     def __setstate__(self, state):
         """
@@ -1024,13 +1043,12 @@ def __setstate__(self, state):
         self._default_state = state[2]
         self._s_state = state[3]
         self._s_default_state = state[4]
-        self._fixed_form_protocol = state[5]
 
         # The following properties need to be set on the internal simulation
         # object
-        self.set_tolerance(*state[6])
-        self.set_min_step_size(state[7])
-        self.set_max_step_size(state[8])
+        self.set_tolerance(*state[5])
+        self.set_min_step_size(state[6])
+        self.set_max_step_size(state[7])
 
     def set_state(self, state):
         """
diff --git a/myokit/_sim/fiber_tissue.py b/myokit/_sim/fiber_tissue.py
index 42a1adff7..002ce55d2 100644
--- a/myokit/_sim/fiber_tissue.py
+++ b/myokit/_sim/fiber_tissue.py
@@ -310,12 +310,12 @@ def __init__(
 
         # Process bindings, remove unsupported bindings, get map of bound
         # variables to internal names.
-        self._bound_variablesf = self._modelf.prepare_bindings({
+        self._bound_variablesf = myokit._prepare_bindings(self._modelf, {
             'time': 'time',
             'pace': 'pace',
             'diffusion_current': 'idiff',
         })
-        self._bound_variablest = self._modelt.prepare_bindings({
+        self._bound_variablest = myokit._prepare_bindings(self._modelt, {
             'time': 'time',
             'pace': 'pace',
             'diffusion_current': 'idiff',
diff --git a/myokit/_sim/openclsim.py b/myokit/_sim/openclsim.py
index 66f609187..e8e489314 100644
--- a/myokit/_sim/openclsim.py
+++ b/myokit/_sim/openclsim.py
@@ -283,7 +283,7 @@ def __init__(
         inputs = {'time': 'time', 'pace': 'pace'}
         if self._diffusion_enabled:
             inputs['diffusion_current'] = 'idiff'
-        self._bound_variables = self._model.prepare_bindings(inputs)
+        self._bound_variables = myokit._prepare_bindings(self._model, inputs)
 
         # Create unique names
         self._model.create_unique_names()
diff --git a/myokit/_sim/pacing.h b/myokit/_sim/pacing.h
index 98b9088bd..b87778c7e 100644
--- a/myokit/_sim/pacing.h
+++ b/myokit/_sim/pacing.h
@@ -620,6 +620,7 @@ typedef int FSys_Flag;
 #define FSys_POPULATE_INVALID_TIMES_DATA    -24
 #define FSys_POPULATE_INVALID_VALUES_DATA   -25
 #define FSys_POPULATE_DECREASING_TIMES_DATA -26
+#define FSys_POPULATE_INVALID_PROTOCOL      -27
 
 /*
  * Sets a python exception based on a fixed-form pacing error flag.
@@ -647,6 +648,9 @@ FSys_SetPyErr(FSys_Flag flag)
         PyErr_SetString(PyExc_Exception, "F-Pacing error: Pacing system not populated.");
         break;
     // Populate
+    case FSys_POPULATE_INVALID_PROTOCOL:
+        PyErr_SetString(PyExc_Exception, "F-Pacing error: Invalid protocol python object passed.");
+        break;
     case FSys_POPULATE_INVALID_TIMES:
         PyErr_SetString(PyExc_Exception, "F-Pacing error: Invalid times array passed.");
         break;
@@ -750,7 +754,7 @@ FSys_Destroy(FSys sys)
  * Returns a fixed-form pacing error flag.
  */
 FSys_Flag
-FSys_Populate(FSys sys, PyObject* times_list, PyObject* values_list)
+FSys_Populate(FSys sys, PyObject* protocol)
 {
     int i;
     Py_ssize_t n;
@@ -758,20 +762,27 @@ FSys_Populate(FSys sys, PyObject* times_list, PyObject* values_list)
     // Check ESys
     if(sys == 0) return FSys_INVALID_SYSTEM;
     if (sys->n_points != -1) return FSys_POPULATED_SYSTEM;
+    if (protocol == Py_None) return FSys_POPULATE_INVALID_PROTOCOL;
+
+    // Get PyList from protocol (will need to decref!)
+    // Cast to (char*) happens because CallMethod accepts a mutable char*
+    // This should have been const char* and has been fixed in python 3
+    PyObject* times_list = PyObject_CallMethod(protocol, (char*)"times", NULL); // Returns a new reference
+    if(times_list == NULL) return FSys_POPULATE_INVALID_PROTOCOL;
+    if(!PyList_Check(times_list)) {
+        Py_DECREF(times_list);
+        return FSys_POPULATE_INVALID_TIMES;
+    }
 
-    // Check input lists
-    if(!PyList_Check(times_list)) return FSys_POPULATE_INVALID_TIMES;
-    if(!PyList_Check(values_list)) return FSys_POPULATE_INVALID_VALUES;
+    // Check and convert times list
     n = PyList_Size(times_list);
-    if (n != PyList_Size(values_list)) return FSys_POPULATE_SIZE_MISMATCH;
-    if (n < 2) return FSys_POPULATE_NOT_ENOUGH_DATA;
-
-    // Convert and check times list
     sys->times = (double*)malloc((size_t)n*sizeof(double));
     for(i=0; i<n; i++) {
         // GetItem and convert --> Borrowed reference so ok not to decref!
         sys->times[i] = PyFloat_AsDouble(PyList_GetItem(times_list, i));
     }
+    Py_DECREF(times_list);  // Finished with the times_list
+
     if (PyErr_Occurred()) {
         free(sys->times); sys->times = NULL;
         return FSys_POPULATE_INVALID_TIMES_DATA;
@@ -783,12 +794,24 @@ FSys_Populate(FSys sys, PyObject* times_list, PyObject* values_list)
         }
     }
 
-    // Convert values list
+    // Check and convert values list
+    PyObject* values_list = PyObject_CallMethod(protocol, (char*)"values", NULL); // Returns a new reference
+    if(values_list == NULL) {
+        free(sys->times); sys->times = NULL;
+        return FSys_POPULATE_INVALID_PROTOCOL;
+    }
+    if(!PyList_Check(values_list) || PyList_Size(values_list) != n) {
+        free(sys->times); sys->times = NULL;
+        Py_DECREF(values_list);
+        return FSys_POPULATE_INVALID_VALUES;
+    }
     sys->values = (double*)malloc((size_t)n*sizeof(double));
     for(i=0; i<n; i++) {
         // GetItem and convert --> Borrowed reference so ok not to decref!
         sys->values[i] = PyFloat_AsDouble(PyList_GetItem(values_list, i));
     }
+    Py_DECREF(values_list); // Finished with the values list
+
     if (PyErr_Occurred()) {
         free(sys->times); sys->times = NULL;
         free(sys->values); sys->values = NULL;
diff --git a/myokit/_sim/rhs.c b/myokit/_sim/rhs.c
index 5e8211107..5a4554a6b 100644
--- a/myokit/_sim/rhs.c
+++ b/myokit/_sim/rhs.c
@@ -49,9 +49,9 @@ tab = '    '
 
 # Process bindings, remove unsupported bindings, get map of bound variables to
 # internal names
-bound_variables = model.prepare_bindings({
+bound_variables = myokit._prepare_bindings(model, {
     'time' : 'engine_time',
-    })
+})
 
 # Get equations
 equations = model.solvable_order()
diff --git a/myokit/formats/ansic/template/cable.c b/myokit/formats/ansic/template/cable.c
index 856b1b553..9fda8168b 100644
--- a/myokit/formats/ansic/template/cable.c
+++ b/myokit/formats/ansic/template/cable.c
@@ -23,11 +23,11 @@ model.create_unique_names()
 w = ansic.AnsiCExpressionWriter()
 
 # Process bindings, remove unsupported bindings
-bound_variables = model.prepare_bindings({
+bound_variables = myokit._prepare_bindings(model, {
     'time'         : 't',
     'pace'         : 'pace',
     'diffusion_current' : 'current'
-    }).keys()
+}).keys()
 
 if model.binding('diffusion_current') is None:
     raise Exception('This exporter requires a variable to be bound to'
diff --git a/myokit/formats/ansic/template/euler.c b/myokit/formats/ansic/template/euler.c
index 689b8b0ad..e7d96c2e6 100644
--- a/myokit/formats/ansic/template/euler.c
+++ b/myokit/formats/ansic/template/euler.c
@@ -29,10 +29,10 @@ model.create_unique_names()
 
 # Process bindings, remove unsupported bindings, get map of bound variables to
 # internal names
-bound_variables = model.prepare_bindings({
+bound_variables = myokit._prepare_bindings(model, {
     'time' : 'time',
     'pace' : 'pace',
-    })
+})
 
 # Get equations
 equations = model.solvable_order()
diff --git a/myokit/formats/ansic/template/sim.c b/myokit/formats/ansic/template/sim.c
index 9b45a8434..08811cb3e 100644
--- a/myokit/formats/ansic/template/sim.c
+++ b/myokit/formats/ansic/template/sim.c
@@ -38,10 +38,10 @@ w.set_lhs_function(v)
 
 # Process bindings, remove unsupported bindings, get map of bound variables to
 # internal names
-bound_variables = model.prepare_bindings({
+bound_variables = myokit._prepare_bindings(model, {
     'time' : 't',
     'pace' : 'pace',
-    })
+})
 
 # Tab
 tab = '    '
diff --git a/myokit/formats/cuda/template/kernel.cu b/myokit/formats/cuda/template/kernel.cu
index 60c73132b..d443cd289 100644
--- a/myokit/formats/cuda/template/kernel.cu
+++ b/myokit/formats/cuda/template/kernel.cu
@@ -71,11 +71,11 @@ initial_values = model.initial_values(True)
 
 # Process bindings, remove unsupported bindings, get map of bound variables to
 # internal names
-bound_variables = model.prepare_bindings({
+bound_variables = myokit._prepare_bindings(model, {
     #'time' : 'time',
     #'pace' : 'pace',
     'diffusion_current' : 'I_diff',
-    })
+})
 
 # Get equations
 equations = model.solvable_order()
diff --git a/myokit/formats/matlab/_exporter.py b/myokit/formats/matlab/_exporter.py
index 01e74b8bf..bf51ede99 100644
--- a/myokit/formats/matlab/_exporter.py
+++ b/myokit/formats/matlab/_exporter.py
@@ -81,7 +81,7 @@ def v(var):
         ew.set_condition_function('ifthenelse')
 
         # Process bound variables
-        bound_variables = model.prepare_bindings({
+        bound_variables = myokit._prepare_bindings(model, {
             'time': 't',
             'pace': 'pace',
         })
diff --git a/myokit/formats/opencl/_exporter.py b/myokit/formats/opencl/_exporter.py
index e5cd10590..a5608e5a8 100644
--- a/myokit/formats/opencl/_exporter.py
+++ b/myokit/formats/opencl/_exporter.py
@@ -80,7 +80,7 @@ def _vars(self, model, protocol):
 
         # Process bindings, remove unsupported bindings, get map of bound
         # variables to internal names.
-        bound_variables = model.prepare_bindings({
+        bound_variables = myokit._prepare_bindings(model, {
             'time': 'time',
             'pace': 'pace',
             'diffusion_current': 'idiff',
diff --git a/myokit/formats/python/template/sim.py b/myokit/formats/python/template/sim.py
index 4989aed5e..d8c2bc714 100644
--- a/myokit/formats/python/template/sim.py
+++ b/myokit/formats/python/template/sim.py
@@ -45,7 +45,7 @@ def c(comp):
 e = w.eq
 
 # Process bindings
-bound_variables = model.prepare_bindings({
+bound_variables = myokit._prepare_bindings(model, {
     'time' : 'engine.time',
     'pace' : 'engine.pace',
     })
diff --git a/myokit/formats/stan/_exporter.py b/myokit/formats/stan/_exporter.py
index 3b689e909..966160f53 100644
--- a/myokit/formats/stan/_exporter.py
+++ b/myokit/formats/stan/_exporter.py
@@ -109,7 +109,7 @@ def v(var):
         ew.set_lhs_function(v)
 
         # Process bound variables
-        bound_variables = model.prepare_bindings({
+        bound_variables = myokit._prepare_bindings(model, {
             'time': 'time',
             'pace': 'pace',
         })


Test patch:
diff --git a/myokit/tests/ansic_fixed_form_pacing.c b/myokit/tests/ansic_fixed_form_pacing.c
index 605e4bc08..4a1227068 100644
--- a/myokit/tests/ansic_fixed_form_pacing.c
+++ b/myokit/tests/ansic_fixed_form_pacing.c
@@ -56,8 +56,7 @@ static PyObject*
 fpacing_init(PyObject *self, PyObject *args)
 {
     // Input arguments
-    PyObject* times;
-    PyObject* values;
+    PyObject* protocol;
 
     FSys_Flag flag;
 
@@ -71,7 +70,7 @@ fpacing_init(PyObject *self, PyObject *args)
     pacing = NULL;
 
     // Check input arguments
-    if (!PyArg_ParseTuple(args, "OO", &times, &values)) {
+    if (!PyArg_ParseTuple(args, "O", &protocol)) {
         PyErr_SetString(PyExc_Exception, "Incorrect input arguments.");
         // Nothing allocated yet, no pyobjects _created_, return directly
         return 0;
@@ -85,7 +84,7 @@ fpacing_init(PyObject *self, PyObject *args)
     // Set up fixed-form pacing
     pacing = FSys_Create(&flag);
     if (flag != FSys_OK) { FSys_SetPyErr(flag); return fpacing_clean(); }
-    flag = FSys_Populate(pacing, times, values);
+    flag = FSys_Populate(pacing, protocol);
     if (flag != FSys_OK) { FSys_SetPyErr(flag); return fpacing_clean(); }
 
     // Done!
diff --git a/myokit/tests/ansic_fixed_form_pacing.py b/myokit/tests/ansic_fixed_form_pacing.py
index e46ee5a28..cfd6a69ad 100644
--- a/myokit/tests/ansic_fixed_form_pacing.py
+++ b/myokit/tests/ansic_fixed_form_pacing.py
@@ -24,7 +24,7 @@ class AnsicFixedFormPacing(myokit.CModule):
     """
     _index = 0
 
-    def __init__(self, times=None, values=None):
+    def __init__(self, protocol):
         super(AnsicFixedFormPacing, self).__init__()
 
         # Unique id
@@ -43,7 +43,7 @@ def __init__(self, times=None, values=None):
         self._sys = self._compile(module_name, fname, args, libs, libd, incd)
 
         # Initialize
-        self._sys.init(times, values)
+        self._sys.init(protocol)
 
     def pace(self, time):
         return self._sys.pace(time)
diff --git a/myokit/tests/test_cmodel.py b/myokit/tests/test_cmodel.py
index d16c08023..df55f27bb 100755
--- a/myokit/tests/test_cmodel.py
+++ b/myokit/tests/test_cmodel.py
@@ -33,7 +33,10 @@ def setUpClass(cls):
         m, p, x = myokit.load(os.path.join(DIR_DATA, 'lr-1991.mmt'))
         cls.model = m
         cls.sensitivities = (['ik1.gK1', 'ikp.IKp'], ['cell.K_o', 'ikp.gKp'])
-        cls.cmodel = myokit.CModel(cls.model, cls.sensitivities)
+        cls.pacing_labels = {}
+        cls.cmodel = myokit.CModel(
+            cls.model, cls.pacing_labels, cls.sensitivities
+        )
 
     def test_sensitivities(self):
         # Test instantiation of cmodel with sensitivities
@@ -41,14 +44,14 @@ def test_sensitivities(self):
         # Bad type
         sens = 'Bad type'
         with self.assertRaisesRegex(ValueError, 'The argument `sensitivities'):
-            myokit.CModel(self.model, sens)
+            myokit.CModel(self.model, self.pacing_labels, sens)
 
         # Empty deps or indeps
         sens = ([], ['some parameter'])
-        m = myokit.CModel(self.model, sens)
+        m = myokit.CModel(self.model, self.pacing_labels, sens)
         self.assertFalse(m.has_sensitivities)
         sens = (['some state'], [])
-        m = myokit.CModel(self.model, sens)
+        m = myokit.CModel(self.model, self.pacing_labels, sens)
         self.assertFalse(m.has_sensitivities)
 
         # Provide sensitivies as Variables
@@ -57,14 +60,14 @@ def test_sensitivities(self):
         p1 = self.model.get('cell.K_o')
         p2 = self.model.get('ikp.gKp')
         sens = ([s1, s2], [p1, p2])
-        m = myokit.CModel(self.model, sens)
+        m = myokit.CModel(self.model, self.pacing_labels, sens)
         self.assertTrue(m.has_sensitivities)
 
         # Provide sensitivities as Names
         sens = (
             [myokit.Name(s1), myokit.Name(s2)],
             [myokit.Name(p1), myokit.Name(p2)])
-        m = myokit.CModel(self.model, sens)
+        m = myokit.CModel(self.model, self.pacing_labels, sens)
         self.assertTrue(m.has_sensitivities)
 
         # Sensitivity of derivative
@@ -72,13 +75,13 @@ def test_sensitivities(self):
         sens = (
             [myokit.Derivative(myokit.Name(s3)), myokit.Name(s2)],
             [myokit.Name(p1), myokit.Name(p2)])
-        m = myokit.CModel(self.model, sens)
+        m = myokit.CModel(self.model, self.pacing_labels, sens)
         self.assertTrue(m.has_sensitivities)
         s3 = 'dot(ik.x)'
         sens = (
             [s3, myokit.Name(s2)],
             [myokit.Name(p1), myokit.Name(p2)])
-        m = myokit.CModel(self.model, sens)
+        m = myokit.CModel(self.model, self.pacing_labels, sens)
         self.assertTrue(m.has_sensitivities)
 
         # Sensitivity of derivative of non-state
@@ -86,21 +89,21 @@ def test_sensitivities(self):
             [myokit.Derivative(myokit.Name(s1)), myokit.Name(s2)],
             [myokit.Name(p1), myokit.Name(p2)])
         with self.assertRaisesRegex(ValueError, 'Sensitivity of '):
-            myokit.CModel(self.model, sens)
+            myokit.CModel(self.model, self.pacing_labels, sens)
 
         # Sensitivity of bound variable
         sens = (
             ['engine.time', myokit.Name(s2)],
             [myokit.Name(p1), myokit.Name(p2)])
         with self.assertRaisesRegex(ValueError, 'Sensitivities cannot'):
-            myokit.CModel(self.model, sens)
+            myokit.CModel(self.model, self.pacing_labels, sens)
 
         # Sensitivity w.r.t. Initial value
         s3 = self.model.get('ik.x')
         sens = (
             [s3, myokit.Name(s2)],
             [myokit.Name(p1), myokit.InitialValue(myokit.Name(s3))])
-        m = myokit.CModel(self.model, sens)
+        m = myokit.CModel(self.model, self.pacing_labels, sens)
         self.assertTrue(m.has_sensitivities)
 
         # Sensitivity w.r.t. initial value of non-state
@@ -108,14 +111,14 @@ def test_sensitivities(self):
             [myokit.Name(s1), myokit.Name(s2)],
             [myokit.Name(p1), myokit.InitialValue(myokit.Name(p2))])
         with self.assertRaisesRegex(ValueError, 'Sensitivity with respect to'):
-            myokit.CModel(self.model, sens)
+            myokit.CModel(self.model, self.pacing_labels, sens)
 
         # Sensitivity w.r.t. non-literal
         sens = (
             [myokit.Name(s1), myokit.Name(s2)],
             [myokit.Name(p1), 'ik.E'])
         with self.assertRaisesRegex(ValueError, 'Sensitivity with respect to'):
-            myokit.CModel(self.model, sens)
+            myokit.CModel(self.model, self.pacing_labels, sens)
 
 
 if __name__ == '__main__':
diff --git a/myokit/tests/test_pacing_system_c.py b/myokit/tests/test_pacing_system_c.py
index 7e3c1e2d5..ac90423a9 100755
--- a/myokit/tests/test_pacing_system_c.py
+++ b/myokit/tests/test_pacing_system_c.py
@@ -144,24 +144,11 @@ def test_fixed_form_pacing_ansic(self):
             import sys
             sys.exit(1)
 
-        # Test input checking
-        times = 1
-        values = [1, 2]
-        self.assertRaises(Exception, AnsicFixedFormPacing)
-        self.assertRaises(Exception, AnsicFixedFormPacing, 1)
-        self.assertRaises(Exception, AnsicFixedFormPacing, 1, 2)
-        self.assertRaises(Exception, AnsicFixedFormPacing, [1], [2])
-        self.assertRaises(
-            Exception, AnsicFixedFormPacing, [1, 2], [2])
-        self.assertRaises(
-            Exception, AnsicFixedFormPacing, [2, 1], [2, 2])
-        AnsicFixedFormPacing([1, 2], [1, 2])
-
         # Test with small lists
         values = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
         times = [0, 0, 1, 1, 1, 2, 2, 2, 3, 4, 5, 7]
         values = list(range(len(times)))
-        pacing = AnsicFixedFormPacing(times, values)
+        pacing = AnsicFixedFormPacing(myokit.TimeSeriesProtocol(times, values))
 
         def test(value, index):
             self.assertEqual(pacing.pace(value), index)
diff --git a/myokit/tests/test_protocol_time_series.py b/myokit/tests/test_protocol_time_series.py
new file mode 100755
index 000000000..9a63e6cdf
--- /dev/null
+++ b/myokit/tests/test_protocol_time_series.py
@@ -0,0 +1,92 @@
+#!/usr/bin/env python3
+#
+# Tests the time series protocol class
+#
+# This file is part of Myokit.
+# See http://myokit.org for copyright, sharing, and licensing details.
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+
+import unittest
+import pickle
+
+import myokit
+
+
+# Unit testing in Python 2 and 3
+try:
+    unittest.TestCase.assertRaisesRegex
+except AttributeError:
+    unittest.TestCase.assertRaisesRegex = unittest.TestCase.assertRaisesRegexp
+
+
+class TimeSeriesProtocolTest(unittest.TestCase):
+    """
+    Tests the TimeSeriesProtocol class.
+    """
+
+    def test_eq(self):
+        # Equality testing
+        p1 = myokit.TimeSeriesProtocol([1, 2], [2, 4])
+        p2 = myokit.TimeSeriesProtocol([1, 2], [2, 4])
+        self.assertEqual(p1, p2)
+        self.assertNotEqual(p1, 1)
+        self.assertEqual(p1, p1)
+
+    def test_pickle(self):
+        # Pickling and unpickling
+        p1 = myokit.TimeSeriesProtocol([1, 2], [2, 4])
+        p2 = pickle.loads(pickle.dumps(p1))
+        self.assertEqual(p1, p2)
+
+    def test_constructor(self):
+        p = myokit.TimeSeriesProtocol([1], [2])
+        self.assertEqual(p.times(), [1])
+        self.assertEqual(p.values(), [2])
+
+        p = myokit.TimeSeriesProtocol([1, 2], [1, 2])
+        self.assertEqual(p.times(), [1, 2])
+        self.assertEqual(p.values(), [1, 2])
+
+        p = myokit.TimeSeriesProtocol([2, 1], [2, 1])
+        self.assertEqual(p.times(), [1, 2])
+        self.assertEqual(p.values(), [1, 2])
+
+        self.assertRaises(  # No specific error from class!!
+            TypeError, myokit.TimeSeriesProtocol, 1, 2)
+        self.assertRaisesRegex(
+            ValueError, 'same size', myokit.TimeSeriesProtocol, [1, 2], [2])
+        self.assertRaisesRegex(
+            ValueError, 'nknown interpolation', myokit.TimeSeriesProtocol,
+            [1, 2], [2, 4], method='cubic'
+        )
+
+    def test_values(self):
+        values = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
+        times = [0, 0, 1, 1, 1, 2, 2, 2, 3, 4, 5, 7]
+        values = list(range(len(times)))
+        pacing = myokit.TimeSeriesProtocol(times, values)
+
+        def test(value, index):
+            self.assertEqual(pacing.pace(value), index)
+
+        test(-1, 0)
+        test(0, 1)
+        test(1, 4)
+        test(2, 7)
+        test(3, 8)
+        test(4, 9)
+        test(5, 10)
+        test(7, 11)
+        test(8, 11)
+        test(1.5, 4.5)
+        test(1.75, 4.75)
+        test(6, 10.5)
+        test(5.5, 10.25)
+
+
+if __name__ == '__main__':
+    import warnings
+    warnings.simplefilter('always')
+    unittest.main()
diff --git a/myokit/tests/test_simulation_cvodes.py b/myokit/tests/test_simulation_cvodes.py
index eed42eb86..7a327d1ea 100755
--- a/myokit/tests/test_simulation_cvodes.py
+++ b/myokit/tests/test_simulation_cvodes.py
@@ -80,6 +80,49 @@ def test_simple(self):
         d2 = self.sim.run(5, log_interval=-5)
         self.assertEqual(d1.time(), d2.time())
 
+    def test_multiple_protocols(self):
+        model = myokit.Model()
+        c = model.add_component('c')
+
+        a = c.add_variable('a')
+        a.set_binding('a')
+        a.set_rhs(0)
+        b = c.add_variable('b')
+        b.set_binding('b')
+        b.set_rhs(0)
+
+        y = c.add_variable('y')
+        t = c.add_variable('t')
+        t.set_binding('time')
+        t.set_rhs(0)
+        y.promote(1)
+        y.set_rhs('- a * y - b * y')
+
+        pa = myokit.Protocol()
+        pa.schedule(1, 0, 0.5)
+
+        pb = myokit.Protocol()
+        pb.schedule(2, 1.0, 0.5)
+
+        sim = myokit.Simulation(model, {'a': pa, 'b': pb})
+        sol = sim.run(2)
+        times = np.array(sol['c.t'])
+        a = np.array(sol['c.a'])
+        a_expect = np.where(times < 0.5, 1., 0.)
+        b = np.array(sol['c.b'])
+        b_expect = np.where((times >= 1.0) & (times < 1.5), 2., 0.)
+        y = np.array(sol['c.y'])
+        y_expect = np.where(times < 0.5, np.exp(-times), 0)
+        y_expect += np.where((times >= 0.5) & (times < 1.0), np.exp(-0.5), 0)
+        y_expect += np.where(
+            (times >= 1.0) & (times < 1.5),
+            np.exp(-0.5) * np.exp(-2 * (times - 1.0)), 0
+        )
+        y_expect += np.where((times >= 1.5), np.exp(-0.5) * np.exp(-1.0), 0)
+        np.testing.assert_array_equal(a, a_expect)
+        np.testing.assert_array_equal(b, b_expect)
+        np.testing.assert_array_almost_equal(y, y_expect, decimal=3)
+
     def test_no_protocol(self):
         # Test running without a protocol.
 
@@ -91,50 +134,73 @@ def test_no_protocol(self):
         # Check if pace was set to zero (see prop 651 / technical docs).
         self.assertTrue(np.all(d['engine.pace'] == 0.0))
 
-    def test_fixed_form_protocol(self):
-        # Test running with a fixed form protocol.
+    def test_wrong_label_set_pacing(self):
+        # Test set_pacing with incorrect label
+        self.sim.reset()
+        self.sim.pre(50)
+        with self.assertRaisesRegex(ValueError, 'Unknown pacing label'):
+            self.sim.set_protocol(None, label='does not exist')
+
+    def test_time_series_protocol(self):
+        # Test running with a time series protocol
 
         n = 10
-        time = list(range(n))
-        pace = [0] * n
-        pace[2:4] = [0.5, 0.5]
+        times = list(range(n))
+        values = [0] * n
+        values[2:4] = [0.5, 0.5]
+        p = myokit.TimeSeriesProtocol(times, values)
 
-        self.sim.set_fixed_form_protocol(time, pace)
+        self.sim.set_protocol(p)
         self.sim.reset()
         d = self.sim.run(n, log_interval=1)
-        self.assertEqual(list(d.time()), time)
-        self.assertEqual(list(d['engine.pace']), pace)
+        self.assertEqual(list(d.time()), times)
+        self.assertEqual(list(d['engine.pace']), values)
 
         # Unset
-        self.sim.set_fixed_form_protocol(None)
+        self.sim.set_protocol(None)
         self.sim.reset()
         d = self.sim.run(n, log_interval=1)
         self.assertEqual(list(d['engine.pace']), [0] * n)
 
         # Reset
-        self.sim.set_fixed_form_protocol(time, pace)
+        self.sim.set_protocol(p)
         self.sim.reset()
         d = self.sim.run(n, log_interval=1)
-        self.assertEqual(list(d.time()), time)
-        self.assertEqual(list(d['engine.pace']), pace)
+        self.assertEqual(list(d.time()), times)
+        self.assertEqual(list(d['engine.pace']), values)
 
         # Unset, replace with original protocol
         self.sim.set_protocol(self.protocol)
         self.sim.reset()
         d = self.sim.run(n, log_interval=1)
-        self.assertNotEqual(list(d['engine.pace']), pace)
+        self.assertNotEqual(list(d['engine.pace']), values)
         self.assertNotEqual(list(d['engine.pace']), [0] * n)
 
-        # Invalid protocols
-        self.assertRaisesRegex(
-            ValueError, 'no times', self.sim.set_fixed_form_protocol,
-            values=pace)
+        # Deprecated version
+        with WarningCollector() as w:
+            self.sim.set_fixed_form_protocol(times, values)
+        self.assertIn('eprecated', w.text())
+        self.sim.reset()
+        d = self.sim.run(n, log_interval=1)
+        self.assertEqual(list(d.time()), times)
+        self.assertEqual(list(d['engine.pace']), values)
         self.assertRaisesRegex(
-            ValueError, 'no values', self.sim.set_fixed_form_protocol,
-            times=time)
+            ValueError, 'No times',
+            self.sim.set_fixed_form_protocol, values=values)
         self.assertRaisesRegex(
-            ValueError, 'same size', self.sim.set_fixed_form_protocol,
-            time, pace[:-1])
+            ValueError, 'No values',
+            self.sim.set_fixed_form_protocol, times=times)
+        self.sim.set_fixed_form_protocol(None)
+        self.sim.reset()
+        d = self.sim.run(n, log_interval=1)
+        self.assertEqual(list(d['engine.pace']), [0] * n)
+
+        # Reset original protocol
+        self.sim.set_protocol(self.protocol)
+        self.sim.reset()
+        d = self.sim.run(n, log_interval=1)
+        self.assertNotEqual(list(d['engine.pace']), values)
+        self.assertNotEqual(list(d['engine.pace']), [0] * n)
 
     def test_in_parts(self):
         # Test running the simulation in parts.
@@ -359,6 +425,47 @@ def test_eval_derivatives(self):
         self.sim.set_state(s1)
         self.assertEqual(d1, self.sim.eval_derivatives())
 
+    def test_eval_derivatives_with_pacing(self):
+        # Test :meth:`Simulation.eval_derivatives()`.
+
+        model = myokit.Model()
+        c = model.add_component('c')
+
+        a = c.add_variable('a')
+        a.set_binding('a')
+        a.set_rhs(0)
+        b = c.add_variable('b')
+        b.set_binding('b')
+        b.set_rhs(0)
+
+        y = c.add_variable('y')
+        t = c.add_variable('t')
+        t.set_binding('time')
+        t.set_rhs(0)
+        y.promote(1)
+        y.set_rhs('- a * y - b * y')
+
+        pa = myokit.Protocol()
+        pa.schedule(1, 0, 0.5)
+
+        pb = myokit.Protocol()
+        pb.schedule(2, 1.0, 0.5)
+
+        sim = myokit.Simulation(model, {'a': pa, 'b': pb})
+        sim.run(1)
+        d1 = sim.eval_derivatives(pacing={'a': 0.5, 'b': 0.5})
+        d2 = sim.eval_derivatives(pacing={'a': 1.5, 'b': 0.5})
+        self.assertNotEqual(d1, d2)
+        d1 = sim.eval_derivatives(pacing={'b': 0.5})
+        d2 = sim.eval_derivatives(pacing={'a': 0.0, 'b': 0.5})
+        self.assertEqual(d1, d2)
+        d1 = sim.eval_derivatives()
+        d2 = sim.eval_derivatives(pacing={'a': 0.0, 'b': 0.0})
+        self.assertEqual(d1, d2)
+        d1 = sim.eval_derivatives(pacing={'a': 0.0, 'b': 0.5})
+        d2 = sim.eval_derivatives(pacing={'aaaaa': 1.0, 'b': 0.5})
+        self.assertEqual(d1, d2)
+
     def test_sensitivities_initial(self):
         # Test setting initial sensitivity values.
 
@@ -830,6 +937,7 @@ def test_cvode_simulation_with_zero_states(self):
         z.promote(0)
 
         # Test without protocol and dynamic logging
+        #myokit.DEBUG_WG = True
         s1 = myokit.Simulation(m1)
         d1 = s1.run(5)
         self.assertEqual(len(d1.time()), 2)