diff --git a/.flake8 b/.flake8
index 23d029b..3b3b4cb 100644
--- a/.flake8
+++ b/.flake8
@@ -1,14 +1,4 @@
 [flake8]
-max_line_length = 120
-ignore =
-    # allow empty line at end of file
-    W391,
-    # break before binary operator - allow either style
-    W503,
-    # break after binary operator - allow either style
-    W504,
-    # missing whitespace around arithmetic operator
-    E226,
-exclude=
-    .git,
-    venv,
+max-line-length = 120
+ignore = W391,W503,W504,E226,E501
+exclude = .git,.*, .venv,docs,build,doctrees,test_output,**/__pycache__,*.egg-info
diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml
index 5ad1ec5..5e5e670 100644
--- a/.github/workflows/pytest.yml
+++ b/.github/workflows/pytest.yml
@@ -23,14 +23,16 @@ jobs:
         python-version: ${{ matrix.python-version }}
     - name: Install dependencies
       run: |
+        sudo apt-get install libsundials-dev -y
         python -m pip install --upgrade pip
         python -m pip install .[test]
     - name: Lint with flake8
       run: |
-        python -m flake8
+        python -m flake8 --exclude=".venv"
     - name: Import sorting with isort
       run: |
-        python -m isort --verbose --check-only --diff markov_builder tests setup.py
+        python -m isort . --show-config
+        python -m isort --verbose --check-only --diff . --skip .venv
     - name: Test with pytest
       run: |
         python -m pytest --cov --cov-config=.coveragerc
diff --git a/.gitignore b/.gitignore
index 02fd1e3..e6619e3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -25,6 +25,8 @@ wheels/
 *.egg
 MANIFEST
 
+markov_builder/_version.py
+
 # PyInstaller
 #  Usually these files are written by a python script from a template
 #  before PyInstaller builds the exe, so as to inject date/other infos into it.
@@ -106,3 +108,8 @@ venv.bak/
 # MacOS extensions
 **.DS_Store
 
+# Test output
+test_output/*
+
+# output
+output/*
diff --git a/.isort.cfg b/.isort.cfg
deleted file mode 100644
index a7b5d22..0000000
--- a/.isort.cfg
+++ /dev/null
@@ -1,14 +0,0 @@
-[isort]
-line_length = 120
-multi_line_output = 3
-# ^ Vertical Hanging Indent
-force_grid_wrap = 4
-# ^ Don't allow more than 3 imports on a single line
-lines_after_imports = 2
-order_by_type = True
-include_trailing_comma = True
-force_single_line = False
-
-default_section = THIRDPARTY
-known_first_party = markov_builder
-skip =
diff --git a/README.md b/README.md
index cbf7c63..30b7650 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,11 @@
 
 
 ## Description
-`markov_builder` is a Python package for constructing Markov models of ion-channel kinetics. These models may be thought of as chemical reaction networks where the transition rate matrix depends on the cell's transmembrane potential (and sometimes the concentration of a drug). This package uses 'networkx' to allow the user to specify the graph associated with the model (model topology). Many models of ion channel kinetics can be expressed in this way. For example, Markov models equivalent to Hodgkin-Huxley style conductance models are straightforward to generate. The transition rates may then be parameterised with any choice of equation, and depend on any number of 'shared variables' such as transmembrane potential, drug concentration, temperature, etc... 
+`markov_builder` is a Python package for constructing Markov models of ion-channel kinetics.
+These models may be thought of as chemical reaction networks where the transition rate matrix depends on the cell's transmembrane potential (and sometimes the concentration of a drug). This package uses 'networkx' to allow the user to specify the graph associated with the model (model topology).
+Many models of ion channel kinetics can be expressed in this way.
+For example, Markov models equivalent to Hodgkin-Huxley style conductance models are straightforward to generate.
+The transition rates may then be parameterised with any choice of equation, and depend on any number of 'shared variables' such as transmembrane potential, drug concentration, temperature, etc...
 
 Given the graph for a model and a parameterisation, you may then:
 1. Output latex describing the model
@@ -10,7 +14,7 @@ Given the graph for a model and a parameterisation, you may then:
 3. Output the system of equations as `sympy` expressions for use in other code
 4. Model drug trapping by composing models together
 5. Simulate a small number of channels using the Gillespie algorithm
-6. Check the model satisfies properties such as microscopic reversibility and connectedness 
+6. Check the model satisfies properties such as microscopic reversibility and connectedness
 
 
 ## Getting Started
@@ -20,7 +24,8 @@ This package requires Python with version >= 3.6 with `numpy`. `pip` is required
 
 ### Installation
 
-It is recommended to install `markov_builder` in a virtual environment to avoid dependency conflicts. To do this, navigate to this repository and run:
+It is recommended to install `markov_builder` in a virtual environment to avoid dependency conflicts.
+To do this, navigate to this repository and run:
 
 ```
 python3 -m pip install --upgrade pip
diff --git a/markov_builder/MarkovChain.py b/markov_builder/MarkovChain.py
index 9c1f687..653e4f7 100644
--- a/markov_builder/MarkovChain.py
+++ b/markov_builder/MarkovChain.py
@@ -1,9 +1,9 @@
-import itertools
 import logging
 from dataclasses import asdict
-from typing import List, Tuple
+from typing import Any, List, Tuple
 
 import myokit
+import myokit.formats.sympy
 import networkx as nx
 import numpy as np
 import pandas as pd
@@ -14,61 +14,114 @@
 from .MarkovStateAttributes import MarkovStateAttributes
 
 
-class MarkovChain():
-    """A class used to construct continuous time Markov Chains/compartmental models using networkx.
+class MarkovChain:
+    """Class describing a CTMC or Markov model of an ion channel.
 
-    Various helper functions are included to generate equations, code for the
-    Markov chains and to test whether the model has certain properties.
+    Describes the topology, parameters and observation function (auxiliary
+    expression) of a ion-channel Markov model.
+
+    - The topology of the model is given by MarkovChain.graph, which has states
+    labels (str) as nodes, and transition rates (str). Optionally, edges may
+    also include a label (str) for diagram purposes.
+
+    - By default states have a flags as defined by the MarkovStateAttributes
+    class, but other properties may be included by specifying a different
+    state_attributes_class Additional properties of states may be included
+
+    - Model parameters relating to a particular transition rate (labels and
+    default values) are stored in rate_expressions, "global variables" which
+    can be referenced by multiple rates.
+
+    - Model parameters relating to the auxiliary expression are stored in
+    MarkovModels.shared_variables, and auxiliary_parameters map variables used
+    in the auxiliary expression to default values.
+
+    - Additional external variables, such as the transmembrane potential or a
+    drug concentration are found in auxiliary_variables with corresponding
+    default values.
+
+
+    :param states: List of states in the model.
+    :param state_attributes_class: A dataclass detailing what data is stored for each state
+    :param seed: Optional random seed to use for random simulations.
+    :param transition_rates: Details transitions and rates to include in the model, each is a tuple, (from_state, to_state, label).
+    :param rate_dictionary: Provides mathematical expressions for transition rates. Each dictionary value is a tuple including the expression and optionally dummy variables and corresponding values. See MarkovChain.parameterise_rates
+    :param auxiliary_expression: Symbol to use to assign to the models auxiliary expression (sometimes known as an observation function)
+    :param shared_variables: Parameters treated as global variables within the model and their corresponding default values
+    :param auxiliary_parameters: Parameters that appear in the auxiliary expression and their default values
+    :param auxiliary_variables: External variables that appear in the auxiliary expression but are not model parameters e.g. transmembrane voltage and drug concentrations
     """
 
-    def __init__(self, states: list = [], state_attributes_class:
-                 MarkovStateAttributes = None, seed: int = None,
-                 name: str = None):
+    def __init__(self, states: list = [], state_attributes_class: Any = MarkovStateAttributes,
+                 dataclass=None, seed: int = None, name: str = None, transition_rates: List = [],
+                 rate_expressions: dict = {}, auxiliary_expression: str = "",
+                 auxiliary_symbol: str = None, shared_variables: dict
+                 = {}, auxiliary_parameters: dict = {}, auxiliary_variables: dict = {}):
 
         # Initialise the graph representing the states. Each directed edge has
         # a `rate` attribute which is a string representing the transition rate
         # between the two nodes
 
         self.graph = nx.DiGraph()
-        if states:
-            self.graph.add_nodes_from(states)
         self.rates = set()
 
+        self.reserved_names = []
+
         # Initialise a random number generator for simulation. Optionally, a
         # seed can be specified.
         self.rng = default_rng(seed)
+
         self.name = name
-        self.shared_rate_variables = set()
 
+        self.shared_variables = shared_variables
         self.rate_expressions = {}
         self.default_values = {}
-
-        self.auxiliary_expression = None
+        self.auxiliary_parameters = {}
+        self.auxiliary_variables = auxiliary_variables
 
         if state_attributes_class is None:
             state_attributes_class = MarkovStateAttributes
 
         self.state_attributes_class = state_attributes_class
-        self.reserved_names = []
-        self.auxiliary_variable = 'auxiliary_expression'
+        self.auxiliary_variable = auxiliary_symbol
 
-        if not issubclass(self.state_attributes_class, MarkovStateAttributes):
-            raise Exception("state_attirbutes_class must be a dataclass")
+        for state in states:
+            self.add_state(state)
 
-    def mirror_model(self, prefix: str, new_rates: bool = False) -> None:
-        """ Duplicate all states and rates in the model such that there are two identical components.
+        # Format auxiliary_expression in case it's using a placeholder for the open state
+        auxiliary_function_defined = auxiliary_expression and auxiliary_symbol and auxiliary_parameters
 
-        The new nodes will be disconnected from the nodes in the original graph.
-        New nodes will be the same as the original nodes but with the prefix
-        prepended. This function may be used to construct drug trapping models.
+        if auxiliary_function_defined:
+            self.define_auxiliary_expression(sp.sympify(auxiliary_expression),
+                                             auxiliary_symbol,
+                                             auxiliary_parameters)
 
-        :param prefix: The prefix to prepend to the new (trapped) nodes and rates (if new_rates is True)
-        :param new_rates: Whether or not to add a prefix to the new transition rates
+        if not auxiliary_expression:
+            auxiliary_expression = sp.sympify("0")
 
-       """
+        if transition_rates:
+            for r in transition_rates:
+                self.add_both_transitions(*r)
 
+        self.parameterise_rates(rate_expressions, shared_variables)
+
+    def mirror_model(self, prefix: str, new_rates: bool = False) -> None:
+        """Duplicate all states and rates in the model such that there are two identical components.
+
+        The new nodes will be disconnected from the nodes in the
+        original graph. New nodes will be the same as the original nodes
+        but with the prefix prepended. This function may be used to
+        construct drug trapping models.
+
+        :param prefix: The prefix to prepend to the new (trapped) nodes
+            and rates (if new_rates is True)
+        :param new_rates: Whether or not to add a prefix to the new
+            transition rates
+        """
         trapped_graph = nx.relabel_nodes(self.graph, dict([(n, "{}{}".format(prefix, n)) for n in self.graph.nodes]))
+
         nx.set_node_attributes(trapped_graph, False, 'open_state')
+        nx.set_node_attributes(trapped_graph, True, 'drug_bound')
 
         if new_rates:
             for frm, to, attr in trapped_graph.edges(data=True):
@@ -86,12 +139,20 @@ def mirror_model(self, prefix: str, new_rates: bool = False) -> None:
 
         self.graph = new_graph
 
+    def set_state_attribute(self, state, **kwargs):
+
+        if state not in self.graph.nodes():
+            raise ValueError(f"State {state} not found in model")
+
+        nx.set_node_attributes(self.graph, {state: kwargs})
+
     def add_open_trapping(self, prefix: str = "d_", new_rates: bool = False) -> None:
         """Construct an open trapping model by mirroring the current model and connecting the open states.
 
-        :param prefix: The prefix to prepend onto the new (trapped) nodes and rates (if new_rates is true.)
-        :param new_rates: Whether or not to create new transition rates for the new mirrored edges
-
+        :param prefix: The prefix to prepend onto the new (trapped)
+            nodes and rates (if new_rates is true.)
+        :param new_rates: Whether or not to create new transition rates
+            for the new mirrored edges
         """
         self.mirror_model(prefix, new_rates)
         self.add_rates(("drug_on", "drug_off"))
@@ -114,9 +175,7 @@ def add_state(self, label, **kwargs) -> None:
         attribute `open_state`.
 
         Example: ``add_state('O', open_state=True)``
-
         """
-
         # Check that the label isn't already in use
         labels_in_use = list(self.graph.nodes())
 
@@ -144,15 +203,10 @@ def add_state(self, label, **kwargs) -> None:
         self.graph.add_node(str(label), **attributes)
 
     def add_rate(self, rate: str) -> None:
-        """
-
-        Add a new transition rate to the model. These are stored in self.rates.
+        """Add a new transition rate to the model. These are stored in self._rates.
 
         :param rate: A string defining the rate to be added
-
-
         """
-
         # Check that the new rate isn't some complicated sympy expression
         # TODO test this and add nice exception
 
@@ -161,36 +215,39 @@ def add_rate(self, rate: str) -> None:
 
         symrate = sp.sympify(rate)
         if len(symrate.atoms()) != 1:
-            raise Exception()
+            raise Exception(f"rate {symrate} consists of multiple symbols")
 
         if rate in self.rates:
-            # TODO
-            raise Exception()
+            raise Exception(f"rate {rate} already exists in model")
         else:
             self.rates.add(rate)
 
     def add_rates(self, rates: list) -> None:
-        """
-        Add a list of rates to the model
+        """Add a list of rates to the model.
 
         :param rates: A list of strings to be added to self.rates
-
         """
         for rate in rates:
             self.add_rate(rate)
 
     def add_transition(self, from_node: str, to_node: str, transition_rate: str,
                        label: str = None, update=False) -> None:
-        """Adds an edge describing the transition rate between `from_node` and `to_node`.
-
-        :param from_node: The state that the transition rate is incident from
-        :param to_node: The state that the transition rate is incident to
-        :param transition rate: A string identifying this transition with a rate from self.rates.
-        :param update: If false and exception will be thrown if an edge between from_node and to_node already exists
+        """Add an edge describing the transition rate between `from_node` and `to_node`.
+
+        :param from_node: The state that the transition rate is incident
+            from
+        :param to_node: The state that the transition rate is incident
+            to
+        :param transition rate: A string identifying this transition
+            with a rate from self._rates.
+        :param update: If True, an edge will be created even if an edge
+            between from_node and to_node already exists
+        :raises RuntimeError: if update is True and a duplicate already
+            edge exists
         """
-
         if from_node not in self.graph.nodes or to_node not in self.graph.nodes:
-            raise Exception("A node wasn't present in the graph ({} or {})".format(from_node, to_node))
+            not_present = [node for node in [from_node, to_node] if node not in self.graph.nodes]
+            raise ValueError(f"A node wasn't present in the graph ({not_present})")
 
         if not isinstance(transition_rate, str):
             transition_rate = str(transition_rate)
@@ -210,17 +267,16 @@ def add_transition(self, from_node: str, to_node: str, transition_rate: str,
             if update:
                 self.graph.add_edge(from_node, to_node, rate=transition_rate, label=label)
             else:
-                raise Exception(f"An edge already exists between {from_node} and {to_node}. \
+                raise RuntimeError(f"An edge already exists between {from_node} and {to_node}. \
                 Edges are {self.graph.edges()}")
         else:
             self.graph.add_edge(from_node, to_node, rate=transition_rate, label=label)
 
     def add_both_transitions(self, frm: str, to: str, fwd_rate: str = None,
                              bwd_rate: str = None, update=True) -> None:
-        """A helper function to add forwards and backwards rates between two
-        states.
+        """Add forwards and backwards rates between two states.
 
-        This is a convenient way to connect new states to the model.
+        This provides a convenient way to connect new states to the model.
 
         :param frm: Either of the two states to be connected.
         :param to: Either of the two states to be connected.
@@ -228,7 +284,6 @@ def add_both_transitions(self, frm: str, to: str, fwd_rate: str = None,
         :param bwd_rate: The transition rate from `to` to `frm`.
         :param update: If false and exception will be thrown if an edge between from_node and to_node already exists
         """
-
         self.add_transition(frm, to, fwd_rate, update=update)
         self.add_transition(to, frm, bwd_rate, update=update)
 
@@ -236,9 +291,12 @@ def get_transition_matrix(self, use_parameters: bool = False,
                               label_order: list = None) -> Tuple[List[str], sp.Matrix]:
         """Compute the Q Matrix of the Markov chain. Q[i,j] is the transition rate between states i and j.
 
-        :param use_parameters: If true substitute in parameters of the transition rates
-        :return: a 2-tuple, labels, and the transition matrix such that the labels column correspond.
-
+        :param use_parameters: If true substitute in parameters of the
+            transition rates
+        :param label_order: If not None, return a transition-rate matrix
+            with columns/rows in the provided order
+        :return: a 2-tuple, labels, and the transition matrix such that
+            the labels column correspond.
         """
         matrix = []
         for current_state in self.graph.nodes:
@@ -250,9 +308,10 @@ def get_transition_matrix(self, use_parameters: bool = False,
                 else:
                     edge = self.graph.get_edge_data(current_state, incident_state)
                     if edge is not None:
-                        row.append(edge["rate"])
+                        rate = edge['rate']
+                        row.append(rate)
                     else:
-                        row.append(0)
+                        row.append(sp.sympify('0'))
             matrix.append(row)
 
         matrix = sp.Matrix(matrix)
@@ -263,47 +322,54 @@ def get_transition_matrix(self, use_parameters: bool = False,
 
         if use_parameters:
             if len(self.rate_expressions) == 0:
-                raise Exception()
+                raise RuntimeError("No rate expressions provided")
             matrix = matrix.subs(self.rate_expressions)
         if label_order is None:
             return list(self.graph.nodes), matrix
         else:
-            if len(label_order) != self.graph.nodes():
+            if extra_labels := [lab for lab in label_order if lab not in self.graph.nodes()]:
+                raise Exception(f"labels: {extra_labels} provided but no found in model. \n"
+                                f"States in model are: {list(self.graph.nodes())}")
+
+            if len(label_order) != len(self.graph.nodes()):
                 raise Exception("Not all states accounted for in label order")
 
-            permutation = [label_order.index(n) for n in self.graph.nodes]
-            matrix_reordered = matrix[permutation, permutation]
-            return label_order, matrix_reordered
+            permutation = [list(self.graph.nodes()).index(s)
+                           for s in label_order]
+            matrix = matrix[permutation, permutation]
+
+            return label_order, matrix
 
     def eval_transition_matrix(self, rates_dict: dict) -> Tuple[List[str], sp.Matrix]:
-        """
-        Evaluate the transition matrix given values for each of the transition rates.
+        """Evaluate the transition matrix given values for each of the transition rates.
 
-        :param rates: A dictionary defining the value of each transition rate e.g rates['K1'] = 1.
+        :param rates: A dictionary defining the value of each transition
+            rate e.g rates['K1'] = 1.
         """
-
         if rates_dict is None:
             Exception('rate dictionary not provided')
 
-        l, Q = self.get_transition_matrix(use_parameters=True)
+        labels, Q = self.get_transition_matrix(use_parameters=True)
         Q_evaled = np.array(Q.evalf(subs=rates_dict)).astype(np.float64)
-        return l, Q_evaled
+        return labels, Q_evaled
 
     def eliminate_state_from_transition_matrix(self, labels: list = None,
                                                use_parameters: bool = False) -> Tuple[sp.Matrix, sp.Matrix]:
-        """Returns a matrix, A, and vector, B, corresponding to a linear ODE system describing the state probabilities.
-
-        Because the state occupancy probabilities must add up to zero, the
-        transition matrix is always singular. We can use this fact to remove
-        one state variable from the system of equations. The labels parameter
-        allows you to choose which variable is eliminated and also the ordering
-        of the states.
-
-        :param labels: A list of labels. The order of which determines the ordering of outputted system.
-        :param use_parameters: If true substitute in parameters of the transition rates
-        :return: A pair of symbolic matrices, A & B, defining a system of ODEs of the format dX/dt = AX + B.
+        """Return a matrix, A, and vector, B, corresponding to a linear ODE system describing the state probabilities.
+
+        Because the state occupancy probabilities must add up to zero,
+        the transition matrix is always singular. We can use this fact
+        to remove one state variable from the system of equations. The
+        labels parameter allows you to choose which variable is
+        eliminated and also the ordering of the states.
+
+        :param labels: A list of labels. The order of which determines
+            the ordering of outputted system.
+        :param use_parameters: If true substitute in parameters of the
+            transition rates
+        :return: A pair of symbolic matrices, A & B, defining a system
+            of ODEs of the format dX/dt = AX + B.
         """
-
         if labels is None:
             labels = list(self.graph.nodes)[:-1]
 
@@ -311,38 +377,25 @@ def eliminate_state_from_transition_matrix(self, labels: list = None,
             if label not in self.graph.nodes():
                 raise Exception(f"Provided label, {label} is not present in the graph")
 
-        _, matrix = self.get_transition_matrix()
-
         eliminated_states = [state for state in self.graph.nodes() if state not in labels]
         assert len(eliminated_states) == 1
         eliminated_state = eliminated_states[0]
 
+        _, matrix = self.get_transition_matrix(label_order=labels + [eliminated_state])
+
         matrix = matrix.T
         shape = sp.shape(matrix)
         assert shape[0] == shape[1]
 
-        # List describing the mapping from self.graph.nodes to labels.
-        # permutation[i] = j corresponds to a mapping which takes
-        # graph.nodes[i] to graph.nodes[j]. Map the row to be eliminated to the
-        # end.
-
-        permutation = [list(self.graph.nodes()).index(n) for n in labels + [eliminated_state]]
-
-        assert len(np.unique(permutation)) == shape[0]
-        matrix = matrix[permutation, permutation]
-
         M = sp.eye(shape[0])
         replacement_row = np.full(shape[0], -1)
-
         M[-1, :] = replacement_row[None, :]
-
         A_matrix = matrix @ M
-
         B_vec = matrix @ sp.Matrix([[0] * (shape[0] - 1) + [1]]).T
 
         if use_parameters:
             if len(self.rate_expressions) == 0:
-                raise Exception()
+                raise Exception("Tried substituting parameters but none found")
             else:
                 A_matrix = A_matrix.subs(self.rate_expressions)
                 B_vec = B_vec.subs(self.rate_expressions)
@@ -350,21 +403,20 @@ def eliminate_state_from_transition_matrix(self, labels: list = None,
         return A_matrix[0:-1, 0:-1], B_vec[0:-1, :]
 
     def get_embedded_chain(self, param_dict: dict = None) -> Tuple[List[str], np.ndarray, np.ndarray]:
-        """Compute the embedded DTMC and associated waiting times given values for each of the transition rates
-
-        :param rates: A dictionary defining the value of each transition rate e.g rates['K1'] = 1.
-
-        :return: 3-tuple: the state labels, the waiting times for each state, and the embedded Markov chain.
+        """Compute the embedded DTMC and associated waiting times given values for each of the transition rates.
 
+        :param rates: A dictionary defining the value of each transition
+            rate e.g rates['K1'] = 1.
+        :return: 3-tuple: the state labels, the waiting times for each
+            state, and the embedded Markov chain.
         """
-
         if param_dict is None:
             param_dict = self.default_values
 
         labs, Q = self.get_transition_matrix()
         _, Q_evaled = self.eval_transition_matrix(param_dict)
 
-        logging.debug("Q is {}".format(Q_evaled))
+        logging.debug("Q is %s", str(Q_evaled))
 
         mean_waiting_times = -1 / np.diagonal(Q_evaled)
 
@@ -385,34 +437,42 @@ def get_embedded_chain(self, param_dict: dict = None) -> Tuple[List[str], np.nda
     def sample_trajectories(self, no_trajectories: int, time_range: list = [0, 1],
                             param_dict: dict = None,
                             starting_distribution: list = None) -> pd.DataFrame:
-        """Samples trajectories of the Markov chain using a Gillespie algorithm.
-
-        :param no_trajectories: The number of simulations to run (number of channels)
-        :param time_range: A range of times durig which to simulate the model
-        :param param_dict: A dictionary defining the (constant) value of each transition rate
-        :param starting_distribution: The number of samples starting in each state. Defaults to an even distribution.
-        :return: A pandas dataframe describing the number of channels in each state for the times in time_range
-
+        """Sample trajectories of the Markov chain using a Gillespie algorithm.
+
+        :param no_trajectories: The number of simulations to run (number
+            of channels)
+        :param time_range: A range of times durig which to simulate the
+            model
+        :param param_dict: A dictionary defining the (constant) value of
+            each transition rate
+        :param starting_distribution: The number of samples starting in
+            each state. Defaults to an even distribution.
+        :return: A pandas dataframe describing the number of channels in
+            each state for the times in time_range
         """
-
         no_nodes = len(self.graph.nodes)
         logging.debug(f"There are {no_nodes} nodes")
 
         if param_dict is not None:
             param_list = self.get_parameter_list()
             # default missing values to those in self.default_values
+
+            _all_default_values = {**self.default_values, **self.shared_variables}
             param_dict = {param: param_dict[param]
                           if param in param_dict
-                          else self.default_values[param]
+                          else _all_default_values[param]
                           for param in param_list}
         else:
             param_dict = self.default_values
 
         if starting_distribution is None:
+            # If there is no user specified starting_distribution, create one
             starting_distribution = np.around(np.array([no_trajectories] * no_nodes) / no_nodes)
             starting_distribution[0] += no_trajectories - starting_distribution.sum()
+        else:
+            starting_distribution = np.array(starting_distribution)
 
-        distribution = starting_distribution
+        distribution = np.around(starting_distribution * no_trajectories / starting_distribution.sum())
 
         labels, mean_waiting_times, e_chain = self.get_embedded_chain(param_dict)
 
@@ -453,57 +513,71 @@ def sample_trajectories(self, no_trajectories: int, time_range: list = [0, 1],
         df = pd.DataFrame(data, columns=['time', *self.graph.nodes], dtype=np.float64)
         return df
 
-    def get_equilibrium_distribution(self, param_dict: dict = None) -> Tuple[List[str], np.array]:
-        """Compute the equilibrium distribution of the CTMC for the specified transition rate values
-
-        :param param_dict: A dictionary specifying the values of each transition rate
-
-        :return: A 2-tuple describing equilibrium distribution and labels defines which entry relates to which state
+    def get_equilibrium_distribution(self, param_dict: dict = {}) -> Tuple[List[str], np.array]:
+        """Compute the equilibrium distribution of the CTMC for the specified transition rate values.
 
+        :param param_dict: A dictionary specifying the values of each
+            transition rate
+        :return: A 2-tuple describing equilibrium distribution and
+            labels defines which entry relates to which state
+        :raises ValueError: If not every necessary parameter is defined
+            in param_dict
         """
         A, B = self.eliminate_state_from_transition_matrix(use_parameters=True)
-
         labels = self.graph.nodes()
-        try:
-            ss = -np.array(A.LUsolve(B).evalf(subs=param_dict)).astype(np.float64)
+        vars_used = [*A.free_symbols, *B.free_symbols]
 
-        except TypeError as exc:
-            logging.warning("Couldn't evaluate equilibrium distribution as float."
-                            "Is every parameter defined?"
-                            "%s" % str(exc))
-            raise exc
+        for _var in vars_used:
+            if str(_var) not in param_dict:
+                raise ValueError(f"Error evaluating equilibrium distribution {_var}\n"
+                                 f"A={A}\nB={B}\nparams={param_dict}\n")
 
-        logging.debug("ss is %s", ss)
+        ss = -np.array(A.LUsolve(B).evalf(subs=param_dict)).astype(np.float64)
         ss = np.append(ss, 1 - ss.sum())
         return labels, ss
 
-    def is_reversible(self) -> bool:
-        """Checks symbolically whether or not the Markov chain is reversible for any
-        set of non-zero transition rate values.
-
-        We assume that all transition rates are always non-zero and follow
-        Colquhoun et al. (2004) https://doi.org/10.1529/biophysj.103.
+    def is_connected(self) -> bool:
+        """Check if the graph is strongly connected.
 
-        :return: A bool which is true if Markov chain is reversible (assuming non-zero transition rates).
+        Returns True if each state can be reached from every other state assuming all transition rates are positive
 
+        :return: A bool which is true if the graph is strongly connected
+            and false otherwise
         """
+        return nx.algorithms.components.is_strongly_connected(self.graph)
 
+    def is_reversible(self) -> bool:
+        """Check symbolically if the Markov chain is reversible for any set of non-zero transition rate values.
+
+        We assume that all transition rates are always non-zero and
+        follow Colquhoun et al. (2004)
+        https://doi.org/10.1529/biophysj.103.
+
+        :return: A bool which is true if Markov chain is reversible
+            (assuming non-zero transition rates).
+        """
         # Digraph must be strongly connected in order for the chain to be
         # reversible. In other words it must be possible to transition from any
         # state to any other state in some finite number of transitions
-        if not nx.algorithms.components.is_strongly_connected(self.graph):
+        if not self.is_connected():
             return False
 
         undirected_graph = self.graph.to_undirected(reciprocal=False, as_view=True)
         cycle_basis = nx.cycle_basis(undirected_graph)
 
         for cycle in cycle_basis:
-            cycle.append(cycle[0])
             logging.debug("Checking cycle {}".format(cycle))
+            cycle.append(cycle[0])
+            iterator = list(zip(cycle[:-1], cycle[1:]))
+            forward_rate_list = [sp.sympify(self.graph.get_edge_data(frm, to)['rate']) for frm, to in iterator]
+            backward_rate_list = [sp.sympify(self.graph.get_edge_data(frm, to)['rate']) for to, frm in iterator]
 
-            iter = list(zip(cycle, itertools.islice(cycle, 1, None)))
-            forward_rate_list = [sp.sympify(self.graph.get_edge_data(frm, to)['rate']) for frm, to in iter]
-            backward_rate_list = [sp.sympify(self.graph.get_edge_data(frm, to)['rate']) for to, frm in iter]
+            logging.debug(self.rate_expressions)
+
+            # Substitute in expressions
+            forward_rate_list = [rate.subs(self.rate_expressions) for rate in forward_rate_list]
+            backward_rate_list = [rate.subs(self.rate_expressions) for rate in
+                                  backward_rate_list]
 
             logging.debug("Rates moving forwards around the cycle are: %s", forward_rate_list)
             logging.debug("Rates moving backwards around the cycle are: %s", backward_rate_list)
@@ -512,29 +586,33 @@ def is_reversible(self) -> bool:
                 logging.debug("Not all rates were specified.")
                 return False
 
-            forward_rate_product = sp.prod(forward_rate_list)
-            backward_rate_product = sp.prod(backward_rate_list)
+            forward_rate_product = sp.prod(forward_rate_list).subs(self.rate_expressions)
+            backward_rate_product = sp.prod(backward_rate_list).subs(self.rate_expressions)
+
             if (forward_rate_product - backward_rate_product).evalf() != 0:
                 return False
         return True
 
     def draw_graph(self, filepath: str = None, show_options: bool =
-                   False, show_rates: bool = False, show_parameters: bool = False):
+                   False, show_rates: bool = False, show_parameters: bool = False,
+                   show_html=False):
         """Visualise the graph as a webpage using pyvis.
 
-        :param filepath: An optional filepath to save the file to. If this is None, will be opened as a webpage instead.
-        :param show_options: Whether or not the options menu should be displayed on the webpage
-        :param show_parameters: Whether or not we should display the transition rates instead of their labels
+        :param filepath: An optional filepath to save the file to. If
+            this is None, will be opened as a webpage instead.
+        :param show_options: Whether or not the options menu should be
+            displayed on the webpage
+        :param show_parameters: Whether or not we should display the
+            transition rates instead of their labels
+        :param show_html: Whether or not to open the outputted html file
+            in the browser
         """
-
         for _, _, data in self.graph.edges(data=True):
             if 'label' not in data or show_rates:
                 data['label'] = data['rate']
             elif show_parameters:
-                if len(self.rate_expressions) == 0:
-                    raise Exception()
-                else:
-                    data['label'] = str(self.rate_expressions[data['rate']])
+                if len(self.rate_expressions) != 0:
+                    data['label'] = str(sp.sympify(data['rate']).subs(self.rate_expressions))
 
         nt = pyvis.network.Network(directed=True)
         nt.from_nx(self.graph)
@@ -543,7 +621,7 @@ def draw_graph(self, filepath: str = None, show_options: bool =
             nt.show_buttons()
         if filepath is not None:
             nt.save_graph(filepath)
-        else:
+        if show_html:
             nt.show('Markov_builder_graph.html')
 
     def substitute_rates(self, rates_dict: dict):
@@ -552,19 +630,20 @@ def substitute_rates(self, rates_dict: dict):
         This function modifies the `rate` attribute of edges in self.graph
 
         :param rates_dict: A dictionary of rates and their corresponding expressions.
-
         """
-
         for rate in rates_dict:
             if rate not in self.rates:
-                raise Exception()
+                raise Exception(f"Tried substituting for rate {rate} but it wasn't found in the model")
+            self.rate_expressions[rate] = rates_dict[rate]
+
         for _, _, d in self.graph.edges(data=True):
             if d['rate'] in rates_dict:
                 if 'label' not in d:
                     d['label'] = d['rate']
                 d['rate'] = str(sp.sympify(d['rate']).subs(rates_dict))
+                self.rate_expressions[d['label']] = sp.sympify(d['rate'])
 
-    def parameterise_rates(self, rate_dict: dict, shared_variables: list = []) -> None:
+    def parameterise_rates(self, rate_dict: dict, shared_variables: dict = {}) -> None:
         """Define a set of parameters for the transition rates.
 
         Parameters declared as
@@ -574,19 +653,9 @@ def parameterise_rates(self, rate_dict: dict, shared_variables: list = []) -> No
         exp(a + b*V) or k = exp(a - b*V) where a and b are dummy variables and
         V is the membrane voltage (a variable shared between transition rates).
 
-        :param rate_dict: A dictionary with a 2-tuple containing an expression and dummy variables for each rate.
-        :param shared_variables: A list of variables that may be shared between transition rates
-
+        :param rate_dict: A dictionary with a tuple containing an expression and optionally dummy variables and corresponding values for each transition rate in the model.
+        :param shared_variables: A dictionary of variables that may be shared between transition rates
         """
-
-        # Check that shared_variables is a list (and not a string!)
-        if isinstance(shared_variables, str):
-            raise TypeError("shared_variables is a string but must be a list")
-
-        for v in shared_variables:
-            if v in self.reserved_names:
-                raise Exception('name %s is reserved', v)
-
         # Validate rate dictionary
         for r in rate_dict:
             if r not in self.rates:
@@ -594,13 +663,23 @@ def parameterise_rates(self, rate_dict: dict, shared_variables: list = []) -> No
 
         rate_expressions = {}
         default_values_dict = {}
+
+        # Construct rate_expressions dict containing all transition rates and corresponding expressions, dummy variables and default values.
+        # Also construct default_values dict which contains all default values for every variable in the model
         for r in self.rates:
             if r in rate_dict:
-                if len(rate_dict[r]) == 2:
+                default_values = []
+                dummy_variables = []
+                if len(rate_dict[r]) == 1:
+                    expression = rate_dict[r][0]
+                elif len(rate_dict[r]) == 2:
                     expression, dummy_variables = rate_dict[r]
-                    default_values = []
-                else:
+                elif len(rate_dict[r]) == 3:
                     expression, dummy_variables, default_values = rate_dict[r]
+                else:
+                    raise ValueError(f"Rate dictionary was malformed. \
+                    Entry with key {r} ({rate_dict[r]}) should be a tuple/list of length 1-3")
+
                 if len(dummy_variables) < len(default_values):
                     raise ValueError("dummy variables list and default values list have mismatching lengths.\
                     Lengths {} and {}".format(len(dummy_variables), len(default_values)))
@@ -608,39 +687,48 @@ def parameterise_rates(self, rate_dict: dict, shared_variables: list = []) -> No
                 expression = sp.sympify(expression)
 
                 for symbol in expression.free_symbols:
-                    variables = list(dummy_variables) + list(shared_variables)
-                    if str(symbol) not in variables:
-                        raise Exception(
-                            f"Symbol, {symbol} was not found in dummy variables or shared_variables, {variables}.")
+                    symbol = str(symbol)
+                    variables = list(dummy_variables) + list(shared_variables.keys())
+                    if symbol not in variables:
+                        # Add symbol to shared variables dictionary
+                        shared_variables[symbol] = None
                 subs_dict = {u: f"{r}_{u}" for i, u in enumerate(dummy_variables)}
-                rate_expressions[r] = sp.sympify(expression).subs(subs_dict)
+                subs_dict = subs_dict | self.rate_expressions
+                rate_expressions[r] = sp.sympify(expression).subs(subs_dict).subs(subs_dict)
 
                 # Add default values to dictionary
                 for u, v in zip(dummy_variables, default_values):
                     new_var_name = f"{r}_{u}"
                     if new_var_name in default_values_dict:
-                        raise Exception(f"A parameter with label {new_var_name} is already present in the model.")
+                        raise ValueError(f"A parameter with label {new_var_name} is already present in the model.")
                     default_values_dict[new_var_name] = v
 
-        self.rate_expressions = {**self.rate_expressions, **rate_expressions}
-        self.default_values = {**self.default_values, **default_values_dict}
+        # TODO repeatedly substitute until no more substitutions possible
+        rate_expressions = {rate: expr.subs(rate_expressions).subs(rate_expressions) for rate, expr in
+                            rate_expressions.items()}
+
+        self.rate_expressions = rate_expressions
+        self.default_values = default_values_dict
+
+        for key in shared_variables.keys():
+            self.default_values[key] = shared_variables[key]
 
-        self.shared_rate_variables = self.shared_rate_variables.union(shared_variables)
+        for key in self.auxiliary_variables:
+            self.default_values[key] = self.auxiliary_variables[key]
 
     def get_parameter_list(self) -> List[str]:
-        """
-        Get a list describing every parameter in the model
+        """Get a list describing every parameter in the model.
 
-        :return: a list of strings corresponding the symbols in self.rate_expressions and self.shared_rate_variables.
+        :return: a list of strings corresponding the symbols in
+            self.rate_expressions and self.shared_rate_variables.
         """
-
         rates = set()
 
         for r in self.rate_expressions:
             for symbol in self.rate_expressions[r].free_symbols:
                 rates.add(str(symbol))
 
-        rates = rates.union([str(sym) for sym in self.shared_rate_variables])
+        rates = rates.union([str(sym) for sym in self.shared_variables])
 
         return sorted(rates)
 
@@ -649,23 +737,21 @@ def generate_myokit_model(self, name: str = "",
                               drug_binding=False, eliminate_state=None) -> myokit.Model:
         """Generate a myokit Model instance describing this Markov model.
 
-        Build a myokit model from this Markov chain using the parameterisation
-        defined by self.rate_expressions. If a rate does not have an entry in
-        self.rate_expressions, it is treated as a constant.
+        Build a myokit model from this Markov chain using the
+        parameterisation defined by self.rate_expressions. If a rate
+        does not have an entry in self.rate_expressions, it is treated
+        as a constant.
 
-        All initial conditions and parameter values should be set before the
-        model is run.
+        All initial conditions and parameter values should be set before
+        the model is run.
 
         :param name: A name to give to the model. Defaults to self.name.
-
-        :param membrane_voltage: A label defining which variable should be treated as the membrane potential.
-
-        :param eliminate_rate: Which rate (if any) to eliminate in order to reduce the number of ODEs in the system.
-
+        :param membrane_voltage: A label defining which variable should
+            be treated as the membrane potential.
+        :param eliminate_rate: Which rate (if any) to eliminate in order
+            to reduce the number of ODEs in the system.
         :return: A myokit.Model built using self
-
         """
-
         if name == "":
             name = self.name
 
@@ -694,8 +780,6 @@ def generate_myokit_model(self, name: str = "",
         model['engine']['time'].set_binding('time')
         model['engine']['time'].set_rhs(0)
 
-        drug_concentration = 'D' if drug_binding else None
-
         # Add parameters to the model
         for parameter in self.get_parameter_list():
             if parameter == membrane_potential:
@@ -704,18 +788,25 @@ def generate_myokit_model(self, name: str = "",
                 model['membrane']['V'].set_rhs(0)
                 model['membrane']['V'].set_binding('pace')
                 comp.add_alias(membrane_potential, model['membrane']['V'])
-            elif drug_binding and parameter == drug_concentration:
-                model.add_component('drug')
-                model['drug'].add_variable('D')
-                model['drug']['D'].set_rhs(0)
-                comp.add_alias(drug_concentration, model['drug']['D'])
             else:
                 comp.add_variable(parameter)
                 if parameter in self.default_values:
                     comp[parameter].set_rhs(self.default_values[parameter])
+                elif parameter in self.auxiliary_variables:
+                    comp[parameter].set_rhs(self.auxiliary_variables[parameter])
+                elif parameter in self.shared_variables:
+                    comp[parameter].set_rhs(self.shared_variables[parameter])
 
         for rate in self.rates:
-            comp.add_variable(rate)
+            free_symbols = sp.parse_expr(rate).free_symbols
+            for symb in free_symbols:
+                if symb not in comp.variables():
+                    try:
+                        comp.add_variable(str(symb))
+                    except myokit.DuplicateName:
+                        # Variable has already been added
+                        pass
+
             if rate in self.rate_expressions:
                 expr = self.rate_expressions[rate]
                 comp[rate].set_rhs(str(expr))
@@ -725,17 +816,21 @@ def generate_myokit_model(self, name: str = "",
             state = self.get_state_symbol(state)
             comp.add_variable(state)
 
+        connected_components = list(nx.connected_components(self.graph.to_undirected()))
+
         # Write down differential equations for the states (unless we chose to
         # eliminate it from the ODE system)
-        for state in states:
-            var = comp[state]
-            var.promote()
-            var.set_rhs(str(d_equations[state]))
-            var.set_state_value(0)
+        for state in self.graph.nodes():
+            state_symbol = self.get_state_symbol(state)
+            var = comp[state_symbol]
+
+            # Give all of the states equal occupancy
+            component = [c for c in connected_components if state in c][0]
+            initial_value = 1.0 / float(len(component))
 
-        # All but one states start with 0 occupancy. If they were all 0 nothing
-        # would happen!
-        comp[states[-1]].set_state_value(1)
+            if state != eliminate_state:
+                var.promote(initial_value)
+                var.set_rhs(str(d_equations[state_symbol]))
 
         # Write equation for eliminated state using the fact that the state
         # occupancies/probabilities must sum to 1.
@@ -749,19 +844,23 @@ def generate_myokit_model(self, name: str = "",
             comp[self.get_state_symbol(eliminate_state)].set_rhs(rhs_str)
 
         # Add auxiliary equation if required
-        if self.auxiliary_expression is not None:
+        if self.auxiliary_expression is not None and self.auxiliary_variable:
             comp.add_variable(self.auxiliary_variable)
-            comp[self.auxiliary_variable].set_rhs(str(self.auxiliary_expression))
+            comp[self.auxiliary_variable].set_rhs(str(self.auxiliary_expression).replace('**', '^'))
+
         return model
 
-    def define_auxiliary_expression(self, expression: sp.Expr, label: str = None, default_values: dict = {}) -> None:
+    def define_auxiliary_expression(self, expression: sp.Expr, label: str =
+                                    None, default_values: dict = {}) -> None:
         """Define an auxiliary output variable for the model.
 
-        :param expression: A sympy expression defining the auxiliary variable
+        :param expression: A sympy expression defining the auxiliary
+            variable
         :param label: A str naming the variable e.g IKr
-        :param default_values: A dictionary of the default values of any parameter used in the auxiliary expression
-
+        :param default_values: A dictionary of the default values of any
+            parameter used in the auxiliary expression
         """
+        self.auxiliary_variables = default_values
         if label in self.graph.nodes() or label in self.reserved_names:
             raise Exception('Name %s not available', label)
         else:
@@ -770,43 +869,51 @@ def define_auxiliary_expression(self, expression: sp.Expr, label: str = None, de
         self.auxiliary_variable = label
 
         if not isinstance(expression, sp.Expr):
-            raise Exception()
+            raise TypeError("Auxiliary expression must be a sympy expression")
 
         state_symbols = [self.get_state_symbol(state) for state in self.graph.nodes()]
-        for symbol in expression.free_symbols:
-            if str(symbol) not in state_symbols:
-                if self.shared_rate_variables is None:
-                    self.shared_rate_variables = set(str(symbol))
-                else:
-                    self.shared_rate_variables.add(str(symbol))
 
         for symbol in default_values:
             symbol = sp.sympify(symbol)
             if symbol not in expression.free_symbols:
-                raise Exception()
-            if symbol in self.default_values:
-                raise Exception()
+                raise Exception(f"Value provided for {symbol} but not found in auxiliarry expression")
+
+        for symbol in expression.free_symbols:
+            if str(symbol) not in state_symbols:
+                symbol = str(symbol)
+                if symbol in default_values.keys():
+                    self.shared_variables[symbol] = default_values[symbol]
+                elif symbol in self.shared_variables:
+                    continue
+                else:
+                    self.shared_variables[symbol] = np.nan
 
-        self.default_values = {**self.default_values, **default_values}
         self.auxiliary_expression = expression
 
+        # TODO check these variables are not elsewhere in the model
+        for key, val in default_values.items():
+            if key not in self.default_values:
+                self.default_values[key] = val
+
     def get_states(self):
+        """Construct a list of the states that are in the model.
+
+        :returns: a list of states
+        """
         return list(self.graph)
 
     def as_latex(self, state_to_remove: str = None, include_auxiliary_expression: bool = False,
                  column_vector=True, label_order: list = None) -> str:
-        """Creates a LaTeX expression describing the Markov chain, its parameters and
-        optionally, the auxiliary equation
-
-        :param state_to_remove: The name of the state (if any) to eliminate from the system.
-
-        :param include_auxiliary_expression: Whether or not to include the auxiliary expression in the output
-        :param column_vector: If False, write the system using row vectors instead of column vectors. Defaults to True
-
+        """Create a LaTeX expression describing the Markov chain, its parameters and optionally, the auxiliary equation.
+
+        :param state_to_remove: The name of the state (if any) to
+            eliminate from the system.
+        :param include_auxiliary_expression: Whether or not to include
+            the auxiliary expression in the output
+        :param column_vector: If False, write the system using row
+            vectors instead of column vectors. Defaults to True
         :returns: A python string containing the relevant LaTeX code.
-
         """
-
         if label_order is not None:
             if state_to_remove is None:
                 if len(label_order) != len(self.graph.nodes()):
@@ -826,7 +933,9 @@ def as_latex(self, state_to_remove: str = None, include_auxiliary_expression: bo
             if column_vector:
                 matrix_str = str(sp.latex(Q.T))
                 eqn = r'\begin{equation}\dfrac{dX}{dt} =  %s X \\end{equation}' % matrix_str
-                X_defn = r'\begin{equation} %s \end{equation}' % sp.latex(sp.Matrix(label_order))
+
+                x_vars = [[sp.sympify(self.get_state_symbol(label))] for label in label_order]
+                X_defn = r'\begin{equation} %s \end{equation}' % sp.latex(sp.Matrix(x_vars))
 
             else:
                 matrix_str = str(sp.latex(Q))
@@ -835,7 +944,7 @@ def as_latex(self, state_to_remove: str = None, include_auxiliary_expression: bo
 
         else:
             if state_to_remove not in map(str, self.graph.nodes()):
-                raise Exception("%s not in model", state_to_remove)
+                raise Exception(f"{state_to_remove} not in model")
             labels = [label for label in label_order]
 
             if len(labels) != len(self.graph.nodes()) - 1:
@@ -877,8 +986,13 @@ def as_latex(self, state_to_remove: str = None, include_auxiliary_expression: bo
                 r'\end{equation}' + "\n where" + return_str
         return return_str
 
-    def get_state_symbol(self, state):
+    def get_state_symbol(self, state: str):
+        """Convert a state found in the model to a suitable identifier for use with sympy.
+
+        :raises ValueError: When state is not present in the model
+        """
         if state in self.graph.nodes():
             return "state_" + state
         else:
-            raise Exception("State not present in model")
+            raise ValueError(f"State not present in model {self.name}: {state}")
+
diff --git a/markov_builder/MarkovStateAttributes.py b/markov_builder/MarkovStateAttributes.py
index 92ba4e8..94dea98 100644
--- a/markov_builder/MarkovStateAttributes.py
+++ b/markov_builder/MarkovStateAttributes.py
@@ -3,11 +3,10 @@
 
 @dataclass
 class MarkovStateAttributes:
-    """A dataclass defining what attributes each state in a MarkovChain should have, and their default values
-
-    In order to include additional attributes, you should create a subclass of
-    this class. Then you can create a MarkovChain with the keyword argument
-    state_attribute_class=my_attributes_subclass
+    """A dataclass defining what attributes each state in a MarkovChain should have, and their default values.
 
+    Additional attributes may be included by specifying another dataclass
     """
+
     open_state: bool = False
+    drug_bound: bool = False
diff --git a/markov_builder/__init__.py b/markov_builder/__init__.py
index bb8cd4d..bbff6f2 100644
--- a/markov_builder/__init__.py
+++ b/markov_builder/__init__.py
@@ -1,6 +1,4 @@
-"""
-Main module for markov_builder
-"""
+"""Main module for markov_builder."""
 import inspect
 import os
 
@@ -30,8 +28,9 @@
 # Expose version number
 #
 def version(formatted=False):
-    """
-    Returns the version number, as a 3-part integer (major, minor, revision).
+    """Returns the version number, as a 3-part integer (major, minor,
+    revision).
+
     If ``formatted=True``, it returns a string formatted version (e.g.
     "markov-builder 1.0.0").
     """
diff --git a/markov_builder/example_models.py b/markov_builder/example_models.py
index 744d9ec..7e2c254 100644
--- a/markov_builder/example_models.py
+++ b/markov_builder/example_models.py
@@ -6,8 +6,9 @@
 
 
 def construct_non_reversible_chain():
-    """Construct a model structure that is known to not satisfy microscopic
-    reversibiliy. This is used for testing.
+    """Construct a model structure that is known to not satisfy microscopic reversibiliy.
+
+    This is used for testing.
     """
     mc = MarkovChain(name='non_reversible_example')
 
@@ -25,10 +26,7 @@ def construct_non_reversible_chain():
 
 
 def construct_four_state_chain():
-    """Construct and parameterise the model introduced by Beattie et al. in
-    https://doi.org/10.1101/100677
-    """
-
+    """Construct and parameterise the model introduced by Beattie et al. in https://doi.org/10.1101/100677."""
     mc = MarkovChain(name='Beattie_model')
     states = ['C', 'I', 'IC']
 
@@ -49,7 +47,7 @@ def construct_four_state_chain():
                        'k_3': positive_rate_expr + ((0.0873, 8.91E-3),),
                        'k_4': negative_rate_expr + ((5.15E-3, 0.003158),)}
 
-    mc.parameterise_rates(rate_dictionary, shared_variables=('V',))
+    mc.parameterise_rates(rate_dictionary)
 
     open_state = mc.get_state_symbol('O')
     auxiliary_expression = sp.sympify(f"g_Kr * {open_state} * (V - E_Kr)")
@@ -62,32 +60,46 @@ def construct_four_state_chain():
 
 
 def construct_mazhari_chain():
-    """Construct the Mazhari model structure for hERG as described in
-    https://doi.org/10.1161/hh1301.093633
-    """
-
+    """Construct the Mazhari model structure for hERG as described in https://doi.org/10.1161/hh1301.093633."""
     mc = MarkovChain(name='Mazhari_model')
 
+    mc.add_state('O', open_state=True)
     for state in ('C1', 'C2', 'C3', 'I'):
         mc.add_state(state)
 
-    mc.add_state('O', open_state=True)
-
     rates = [('C1', 'C2', 'a0', 'b0'), ('C2', 'C3', 'kf', 'kb'), ('C3', 'O', 'a1', 'b1'),
              ('O', 'I', 'ai', 'bi'), ('I', 'C3', 'psi', 'ai3')]
 
     for r in rates:
         mc.add_both_transitions(*r)
 
+    constant_rate_expr = ('a', ('a',))
+    rate_dictionary = {
+        'a0': positive_rate_expr + ((0.0069, 0.0272),),
+        'a1': positive_rate_expr + ((0.0218, 0.0262),),
+        'ai': positive_rate_expr + ((0.662, 0.012),),
+        'ai3': positive_rate_expr + ((1.29e-5, 2.71e-6),),
+        'b0': negative_rate_expr + ((0.0227, 0.0431),),
+        'b1': negative_rate_expr + ((0.0009, 0.0269),),
+        'bi': negative_rate_expr + ((0.0059, 0.0443),),
+        'kf': constant_rate_expr + ((0.0266,),),
+        'kb': constant_rate_expr + ((0.1348,),)
+    }
+
+    open_state = mc.get_state_symbol('O')
+    auxiliary_expression = sp.sympify(f"g_Kr * {open_state} * (V + E_Kr)")
+    mc.define_auxiliary_expression(auxiliary_expression, 'I_kr',
+                                   {'g_Kr': 0.1524,
+                                    'E_Kr': -88})
+
     mc.substitute_rates({'psi': '(ai3*bi*b1)/(a1*ai)'})
+    mc.parameterise_rates(rate_dictionary)
 
     return mc
 
 
 def construct_wang_chain():
-    """Construct the Wang model structure for hERG as described in
-    https://doi.org/10.1111/j.1469-7793.1997.045bl.x
-    """
+    """Construct the Wang model structure for hERG as described in https://doi.org/10.1111/j.1469-7793.1997.045bl.x."""
     mc = MarkovChain(name='Wang_model')
 
     mc.add_state('O', open_state=True)
@@ -95,8 +107,10 @@ def construct_wang_chain():
     for state in ('C1', 'C2', 'C3', 'I'):
         mc.add_state(state)
 
-    rates = [('C1', 'C2', 'a_a0', 'b_a0'), ('C2', 'C3', 'k_f', 'k_b'), ('C3', 'O', 'a_a1', 'b_a1'),
-             ('O', 'I', 'a_1', 'b_1')]
+        rates = [('C1', 'C2', 'a_a0', 'b_a0'),
+                 ('C2', 'C3', 'k_f', 'k_b'),
+                 ('C3', 'O', 'a_a1', 'b_a1'),
+                 ('O', 'I', 'a_1', 'b_1')]
 
     for r in rates:
         mc.add_both_transitions(*r)
@@ -104,38 +118,34 @@ def construct_wang_chain():
     constant_rate_expr = ('a', ('a',))
 
     rate_dictionary = {'a_a0': positive_rate_expr + ((0.022348, 0.01176),),
-                       'b_a0': negative_rate_expr + ((0.047002, 0.0631),),
-                       'k_f': constant_rate_expr + ((0.023761,),),
-                       'k_b': constant_rate_expr + ((0.036778,),),
                        'a_a1': positive_rate_expr + ((0.013733, 0.038198),),
+                       'b_a0': negative_rate_expr + ((0.047002, 0.0631),),
                        'b_a1': negative_rate_expr + ((0.0000689, 0.04178),),
 
                        # Using 2mmol KCl values
                        'a_1': positive_rate_expr + ((0.090821, 0.023391),),
-                       'b_1': negative_rate_expr + ((0.006497, 0.03268),)
+                       'b_1': negative_rate_expr + ((0.006497, 0.03268),),
+                       'k_f': constant_rate_expr + ((0.023761,),),
+                       'k_b': constant_rate_expr + ((0.036778,),),
                        }
 
-    mc.parameterise_rates(rate_dictionary, shared_variables=('V',))
-
     open_state = mc.get_state_symbol('O')
-
     auxiliary_expression = sp.sympify(f"g_Kr * {open_state} * (V + E_Kr)")
     mc.define_auxiliary_expression(auxiliary_expression, 'I_kr',
                                    {'g_Kr': 0.1524,
                                     'E_Kr': -88})
+    mc.parameterise_rates(rate_dictionary)
+
     return mc
 
 
 def construct_HH_model(n: int, m: int, name: str = None):
-    """ Construct a Markov model equivalent to a Hodgkin-Huxley conductance models
+    """Construct a Markov model equivalent to a Hodgkin-Huxley conductance models.
 
     :param n: The number of activation gates in the model
     :param m: The number of inactivation gates in the model
-
     :return: A MarkovChain with n x m states
-
     """
-
     if n < 2 or m < 2:
         raise Exception()
 
@@ -179,23 +189,27 @@ def construct_HH_model(n: int, m: int, name: str = None):
 
 
 def construct_kemp_model():
-    """Construct and parameterise the model introduced by Kemp et al. in
+    """Construct and parameterise the model introduced by Kemp et al.
+
+    in
     https://doi.org/10.1085/jgp.202112923
     """
-
     mc = MarkovChain(name='Kemp_model')
 
     # Now the conducting state
-    mc.add_state('O', open_state=True)
+    mc.add_state('O1', open_state=True)
+
+    # Add a second conducting state
+    mc.add_state('O2', open_state=True)
 
     # First add the non-conducting states
-    for state in ('IO', 'C1', 'IC1', 'C2', 'IC2'):
+    for state in ('I', 'C1', 'C2'):
         mc.add_state(state)
 
     rates = [
-        ('O', 'IO', 'b_h', 'a_h'), ('C1', 'IC1', 'b_h', 'a_h'), ('C2', 'IC2', 'b_h', 'a_h'),
-        ('O', 'C1', 'b_2', 'a_2'), ('C1', 'C2', 'b_1', 'a_1'),
-        ('IO', 'IC1', 'b_2', 'a_2'), ('IC1', 'IC2', 'b_1', 'a_1')
+        ('O2', 'I', 'b_h', 'a_h'),
+        ('O1', 'C1', 'b_2', 'a_2'),
+        ('C1', 'C2', 'b_1', 'a_1')
     ]
 
     for r in rates:
@@ -217,14 +231,15 @@ def construct_kemp_model():
         'b_h': positive_rate_expr + ((2.70e-01, 1.58e-02),),
     }
 
-    mc.parameterise_rates(rate_dictionary, shared_variables=('V',))
-
-    open_state = mc.get_state_symbol('O')
+    open_state1 = mc.get_state_symbol('O1')
+    open_state2 = mc.get_state_symbol('O2')
 
-    auxiliary_expression = sp.sympify(f"g_Kr * {open_state} * (V + E_Kr)")
+    auxiliary_expression = sp.sympify(f"g_Kr * {open_state1} * {open_state2}\
+    * (V + E_Kr)")
     mc.define_auxiliary_expression(auxiliary_expression, 'I_kr',
                                    {
                                        'g_Kr': 7.05e-02,  # Use conductance from Cell 2
                                        'E_Kr': -88,  # -88mV chosen arbitrarily
                                    })
+    mc.parameterise_rates(rate_dictionary)
     return mc
diff --git a/markov_builder/models/__init__.py b/markov_builder/models/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/markov_builder/models/thirty_models/__init__.py b/markov_builder/models/thirty_models/__init__.py
new file mode 100644
index 0000000..665c79a
--- /dev/null
+++ b/markov_builder/models/thirty_models/__init__.py
@@ -0,0 +1,40 @@
+from .model0 import model_00
+from .model1 import model_01
+from .model2 import model_02
+from .model3 import model_03
+from .model4 import model_04
+from .model5 import model_05
+from .model6 import model_06
+from .model7 import model_07
+from .model8 import model_08
+from .model9 import model_09
+from .model10 import model_10
+from .model11 import model_11
+from .model12 import model_12
+from .model13 import model_13
+from .model14 import model_14
+from .model15 import model_15
+from .model16 import model_16
+from .model17 import model_17
+from .model18 import model_18
+from .model19 import model_19
+from .model20 import model_20
+from .model21 import model_21
+from .model22 import model_22
+from .model23 import model_23
+from .model24 import model_24
+from .model25 import model_25
+from .model26 import model_26
+from .model27 import model_27
+from .model28 import model_28
+from .model29 import model_29
+from .model30 import model_30
+
+
+__all__ = ['model_00', 'model_01', 'model_02', 'model_03', 'model_04',
+           'model_05', 'model_06', 'model_07', 'model_08', 'model_09',
+           'model_10', 'model_11', 'model_12', 'model_13', 'model_14',
+           'model_15', 'model_16', 'model_17', 'model_18', 'model_19',
+           'model_20', 'model_21', 'model_22', 'model_23', 'model_24',
+           'model_25', 'model_26', 'model_27', 'model_28', 'model_29',
+           'model_30']
diff --git a/markov_builder/models/thirty_models/model0.py b/markov_builder/models/thirty_models/model0.py
new file mode 100644
index 0000000..69e4157
--- /dev/null
+++ b/markov_builder/models/thirty_models/model0.py
@@ -0,0 +1,40 @@
+import numpy as np
+
+from markov_builder import MarkovChain
+from markov_builder.rate_expressions import (negative_rate_expr,
+                                             positive_rate_expr)
+
+
+class model_00(MarkovChain):
+    description = ""
+    states = ('O', 'C', 'I', 'IC')
+    rates = [('O', 'C', 'k_2', 'k_1'),
+             ('I', 'IC', 'k_2', 'k_1'),
+             ('O', 'I', 'k_3', 'k_4'),
+             ('C', 'IC', 'k_3', 'k_4')]
+
+    shared_variables_dict = {'V': np.nan}
+
+    shared_variables_dict = {
+        'V': np.nan,
+        'g_Kr': 0.1524,
+    }
+
+    rate_dictionary = {'k_1': positive_rate_expr + ((2.26E-4, 6.99E-2),),
+                       'k_2': negative_rate_expr + ((3.44E-5, 5.460E-2),),
+                       'k_3': positive_rate_expr + ((0.0873, 8.91E-3),),
+                       'k_4': negative_rate_expr + ((5.15E-3, 0.03158),)}
+    auxiliary_expression = "g_Kr * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+    auxiliary_parameters_dict = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters_dict
+                         )
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model1.py b/markov_builder/models/thirty_models/model1.py
new file mode 100644
index 0000000..7eb662c
--- /dev/null
+++ b/markov_builder/models/thirty_models/model1.py
@@ -0,0 +1,37 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_01(MarkovChain):
+    description = ""
+    states = ('O', 'C')
+    rates = [('O', 'C', 'k_21', 'k_12')]
+
+    rate_dictionary = {
+        'k_12': ('p1 * exp(p2*V)',),
+        'k_21': ('p3 * exp(-p4*V)',)
+    }
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'g_Kr': 0.1524,
+                             }
+
+    auxiliary_expression = "g_Kr * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model10.py b/markov_builder/models/thirty_models/model10.py
new file mode 100644
index 0000000..85cd517
--- /dev/null
+++ b/markov_builder/models/thirty_models/model10.py
@@ -0,0 +1,56 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_10(MarkovChain):
+    description = ""
+    states = ('O', 'C1', 'C2', 'O2', 'I')
+    rates = [
+        ('C1', 'O', 'a2', 'b2'),
+        ('C2', 'C1', 'a1', 'b1'),
+        ('I', 'O2', 'ah', 'bh')
+    ]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 0.08730,
+                             'p10': 8.91e-3,
+                             'p11': 5.15e-3,
+                             'p12': 0.03158,
+                             'p13': 0.15240,
+                             }
+
+    rate_dictionary = {
+        'a1': ('p1 * exp(p2*V)',),
+        'b1': ('p3 * exp(-p4*V)',),
+        'bh': ('p5 * exp(p6*V)',),
+        'ah': ('p7 * exp(-p8*V)',),
+        'a2': ('p9 * exp(p10*V)',),
+        'b2': ('p11 * exp(-p12*V)',),
+    }
+
+    auxiliary_expression = "p13 * state_O * state_O2 * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
+        self.set_state_attribute('O2', open_state=True)
diff --git a/markov_builder/models/thirty_models/model11.py b/markov_builder/models/thirty_models/model11.py
new file mode 100644
index 0000000..3a670ed
--- /dev/null
+++ b/markov_builder/models/thirty_models/model11.py
@@ -0,0 +1,67 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_11(MarkovChain):
+    description = ""
+    states = ('O', 'I', 'IC1', 'IC2', 'C2', 'C1')
+    rates = [
+        ('O', 'I', 'ah', 'bh'),
+        ('C1', 'O', 'a1', 'b1'),
+        ('IC1', 'I', 'a3', 'b3'),
+        ('IC2', 'IC1', 'a4', 'b4'),
+        ('C2', 'C1', 'a2', 'b2'),
+        ('C2', 'IC2', 'ah', 'bh'),
+        ('C1', 'IC1', 'ah', 'bh')
+    ]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 2.26e-4,
+                             'p10': 0.06990,
+                             'p11': 3.45e-5,
+                             'p12': 0.05462,
+                             'p13': 0.08730,
+                             'p14': 8.91e-3,
+                             'p15': 5.15e-3,
+                             'p16': 0.03158,
+                             'p17': 0.15240
+                             }
+
+    rate_dictionary = {
+        'a1': ('p1 * exp( p2 * V)',),
+        'b1': ('p3 * exp(-p4 * V)',),
+        'ah': ('p5 * exp( p6 * V)',),
+        'bh': ('p7 * exp(-p8 * V)',),
+        'a2': ('p9 * exp( p10 * V)',),
+        'b2': ('p11 * exp(-p12 * V)',),
+        'a3': ('p13 * exp( p14 * V)',),
+        'b3': ('(a3*b1)/a1',),
+        'a4': ('p15 * exp(p16* V)',),
+        'b4': ('(a4*b2)/a2',)
+    }
+
+    auxiliary_expression = "p17 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model12.py b/markov_builder/models/thirty_models/model12.py
new file mode 100644
index 0000000..b2fd465
--- /dev/null
+++ b/markov_builder/models/thirty_models/model12.py
@@ -0,0 +1,61 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_12(MarkovChain):
+    description = ""
+    states = ('O', 'I', 'IC1', 'IC2', 'C2', 'C1')
+    rates = [
+        ('O', 'I', 'a1', 'b1'),
+        ('C1', 'O', 'am', '2*bm'),
+        ('IC1', 'I', 'am', '2*bm'),
+        ('IC2', 'IC1', '2*am', 'bm'),
+        ('C2', 'C1', '2*am', 'bm'),
+        ('C2', 'IC2', 'a3', 'b3'),
+        ('C1', 'IC1', 'a2', 'b2')
+    ]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 2.26e-4,
+                             'p10': 0.06990,
+                             'p11': 3.45e-5,
+                             'p12': 0.05462,
+                             'p13': 0.15240
+                             }
+
+    rate_dictionary = {
+        'am': ('p1 * exp( p2 * V)',),
+        'bm': ('p3 * exp(-p4 * V)',),
+        'a1': ('p5 * exp( p6 * V)',),
+        'b1': ('p7 * exp(-p8 * V)',),
+        'a2': ('p9 * exp( p10 * V)',),
+        'a3': ('p11 * exp( p12 * V)',),
+        'b2': ('a2 * b1 / a1',),
+        'b3': ('a3 * b2 / a2',),
+    }
+
+    auxiliary_expression = "p13 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -85}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model13.py b/markov_builder/models/thirty_models/model13.py
new file mode 100644
index 0000000..0961258
--- /dev/null
+++ b/markov_builder/models/thirty_models/model13.py
@@ -0,0 +1,79 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_13(MarkovChain):
+    description = ""
+    states = ('O', 'I', 'IC1', 'IC2', 'C2', 'C1')
+    rates = [
+        ('O', 'I', 'a3', 'b3'),
+        ('C1', 'O', 'a2', 'b2'),
+        ('IC1', 'I', 'a7', 'b7'),
+        ('IC2', 'IC1', 'a6', 'b6'),
+        ('C2', 'C1', 'a1', 'b1'),
+        ('C2', 'IC2', 'a5', 'b5'),
+        ('C1', 'IC1', 'a4', 'b4')
+    ]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 2.26e-4,
+                             'p10': 0.06990,
+                             'p11': 3.45e-5,
+                             'p12': 0.05462,
+                             'p13': 2.26e-4,
+                             'p14': 0.06990,
+                             'p15': 3.45e-5,
+                             'p16': 0.05462,
+                             'p17': 0.08730,
+                             'p18': 8.91e-3,
+                             'p19': 5.15e-3,
+                             'p20': 0.03158,
+                             'p21': 2.26e-4,
+                             'p22': 0.06990,
+                             'p23': 3.45e-5,
+                             'p24': 0.05462,
+                             'p25': 0.15240,
+                             }
+
+    rate_dictionary = {
+        'a1': ('p1 * exp( p2 * V)',),
+        'b1': ('p3 * exp(-p4 * V)',),
+        'a2': ('p5 * exp( p6 * V)',),
+        'b2': ('p7 * exp(-p8 * V)',),
+        'a3': ('p9 * exp( p10 * V)',),
+        'b3': ('p11 * exp(-p12 * V)',),
+        'a4': ('p13 * exp( p14 * V)',),
+        'b4': ('(a7*b3*b2*a4)/(a2*a3*b7)',),
+        'a5': ('p15 * exp( p16 * V)',),
+        'b5': ('(a5*a6*b4*b1)/(a1*a4*b6)',),
+        'a6': ('p17 * exp( p18 * V)',),
+        'b6': ('p19 * exp(-p20 * V)',),
+        'a7': ('p21 * exp(p22 * V)',),
+        'b7': ('p23 * exp(-p24 * V)',),
+    }
+
+    auxiliary_expression = "p25 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -85}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model14.py b/markov_builder/models/thirty_models/model14.py
new file mode 100644
index 0000000..27e39f7
--- /dev/null
+++ b/markov_builder/models/thirty_models/model14.py
@@ -0,0 +1,50 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_14(MarkovChain):
+    description = ""
+    states = ('O', 'C2', 'C1', 'C3', 'I')
+    rates = [
+        ('C1', 'O', 'am', '3*bm'),
+        ('C2', 'C1', '2*am', '2*bm'),
+        ('C3', 'C2', '3*am', 'bm'),
+        ('O', 'I', 'b1', 'a1'),
+    ]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 0.1524
+                             }
+
+    rate_dictionary = {
+        'am': ('p1 * exp( p2 * V)',),
+        'bm': ('p3 * exp(-p4 * V)',),
+        'b1': ('p5 * exp(p6 * V)',),
+        'a1': ('p7 * exp(-p8 * V)',),
+    }
+
+    auxiliary_expression = "p9 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -85}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model15.py b/markov_builder/models/thirty_models/model15.py
new file mode 100644
index 0000000..490c445
--- /dev/null
+++ b/markov_builder/models/thirty_models/model15.py
@@ -0,0 +1,81 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_15(MarkovChain):
+    description = "IKr Markov model 15 (Sanmitra Ghosh)"
+
+    # States
+    states = ('C3', 'C2', 'C1', 'O', 'I')
+
+    # Transitions derived directly from dot() equations
+    rates = [
+        ('C3', 'C2', 'a1', 'b1'),
+        ('C2', 'C1', 'a2', 'b2'),
+        ('C1', 'O', 'a3', 'b3'),
+        ('O', 'I', 'a4', 'b4'),
+    ]
+
+    # Rate expressions (verbatim from .mmt)
+    rate_dictionary = {
+        'a1': ('p1 * exp(p2 * V)',),
+        'b1': ('p3 * exp(-p4 * V)',),
+
+        'a2': ('p5 * exp(p6 * V)',),
+        'b2': ('p7 * exp(-p8 * V)',),
+
+        'a3': ('p9 * exp(p10 * V)',),
+        'b3': ('p11 * exp(-p12 * V)',),
+
+        'a4': ('p13 * exp(p14 * V)',),
+        'b4': ('p15 * exp(-p16 * V)',),
+    }
+
+    # Shared variables and parameters
+    shared_variables_dict = {
+        'V': nan,
+
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+
+        'p9': 2.26e-4,
+        'p10': 0.06990,
+        'p11': 3.45e-5,
+        'p12': 0.05462,
+
+        'p13': 0.08730,
+        'p14': 8.91e-3,
+        'p15': 5.15e-3,
+        'p16': 0.03158,
+
+        'p17': 0.15240,
+    }
+
+    # Current expression
+    auxiliary_expression = 'p17 * state_O * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {
+        'E_Kr': -88,
+    }
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters,
+        )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model16.py b/markov_builder/models/thirty_models/model16.py
new file mode 100644
index 0000000..f9404e2
--- /dev/null
+++ b/markov_builder/models/thirty_models/model16.py
@@ -0,0 +1,80 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_16(MarkovChain):
+    description = "IKr Markov model 16 (Wang et al. 1997, Sanmitra Ghosh)"
+
+    # States
+    states = ('C3', 'C2', 'C1', 'O', 'I')
+
+    # Transitions derived from dot() equations
+    rates = [
+        ('C3', 'C2', 'a1', 'b1'),
+        ('C2', 'C1', 'a2', 'b2'),
+        ('C1', 'O', 'a3', 'b3'),
+        ('O', 'I', 'a4', 'b4'),
+    ]
+
+    # Rate expressions (verbatim from .mmt)
+    rate_dictionary = {
+        'a1': ('p1 * exp(p2 * V)',),
+        'b1': ('p3 * exp(-p4 * V)',),
+
+        'a3': ('p5 * exp(p6 * V)',),
+        'b3': ('p7 * exp(-p8 * V)',),
+
+        'a4': ('p9 * exp(p10 * V)',),
+        'b4': ('p11 * exp(-p12 * V)',),
+
+        # Voltage-independent rates
+        'a2': ('p13',),
+        'b2': ('p14',),
+    }
+
+    # Shared variables and parameters
+    shared_variables_dict = {
+        'V': nan,
+
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+
+        'p9': 0.08730,
+        'p10': 8.91e-3,
+        'p11': 5.15e-3,
+        'p12': 0.03158,
+
+        'p13': 0.08730,
+        'p14': 8.91e-3,
+
+        'p15': 0.15240,
+    }
+
+    # Current expression
+    auxiliary_expression = 'p15 * state_O * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {
+        'E_Kr': -85,
+    }
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters,
+        )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model17.py b/markov_builder/models/thirty_models/model17.py
new file mode 100644
index 0000000..a2cb54a
--- /dev/null
+++ b/markov_builder/models/thirty_models/model17.py
@@ -0,0 +1,76 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_17(MarkovChain):
+    description = "IKr Markov model (Sanmitra Ghosh)"
+
+    # States
+    states = ('C3', 'C2', 'C1', 'O', 'I')
+
+    # Transitions inferred directly from dot() equations
+    rates = [
+        # C3 <-> C2
+        ('C3', 'C2', '3*am', 'bm'),
+
+        # C2 <-> C1
+        ('C2', 'C1', '2*am', '2*bm'),
+
+        # C1 <-> O
+        ('C1', 'O', 'am', '3*bm'),
+
+        # C1 <-> I
+        ('C1', 'I', 'a2', 'b2'),
+
+        # O <-> I
+        ('O', 'I', 'a1', 'b1'),
+    ]
+
+    # Rate expressions (verbatim from .mmt)
+    rate_dictionary = {
+        'am': ('p1 * exp(p2 * V)',),
+        'bm': ('p3 * exp(-p4 * V)',),
+        'a1': ('p5 * exp(p6 * V)',),
+        'b1': ('p7 * exp(-p8 * V)',),
+        'a2': ('p9 * exp(p10 * V)',),
+        'b2': ('(a2 * 3 * bm * b1) / (am * a1)',),
+    }
+
+    # Shared variables and parameters
+    shared_variables_dict = {
+        'V': nan,
+
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+        'p9': 5.15e-3,
+        'p10': 0.03158,
+        'p11': 0.15240,
+    }
+
+    # Current expression
+    auxiliary_expression = 'p11 * state_O * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {
+        'E_Kr': -85,
+    }
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters,
+        )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model18.py b/markov_builder/models/thirty_models/model18.py
new file mode 100644
index 0000000..fa1c251
--- /dev/null
+++ b/markov_builder/models/thirty_models/model18.py
@@ -0,0 +1,97 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_18(MarkovChain):
+    description = "IKr Markov model 18 (Sanmitra Ghosh)"
+
+    # States
+    states = ('C3', 'C2', 'C1', 'O', 'I')
+
+    # Transitions derived directly from dot() equations
+    rates = [
+        # C3 <-> C2
+        ('C3', 'C2', 'a1', 'b1'),
+
+        # C2 <-> C1
+        ('C2', 'C1', 'a2', 'b2'),
+
+        # C1 <-> O
+        ('C1', 'O', 'a3', 'b3'),
+
+        # O <-> I
+        ('O', 'I', 'a4', 'b4'),
+
+        # C1 <-> I
+        ('C1', 'I', 'a5', 'b5'),
+    ]
+
+    # Rate expressions (verbatim from .mmt)
+    rate_dictionary = {
+        'a1': ('p1 * exp(p2 * V)',),
+        'b1': ('p3 * exp(-p4 * V)',),
+
+        'a2': ('p5 * exp(p6 * V)',),
+        'b2': ('p7 * exp(-p8 * V)',),
+
+        'a3': ('p9 * exp(p10 * V)',),
+        'b3': ('p11 * exp(-p12 * V)',),
+
+        'a4': ('p13 * exp(p14 * V)',),
+        'b4': ('p15 * exp(-p16 * V)',),
+
+        'a5': ('p17 * exp(p18 * V)',),
+        'b5': ('(a5 * b3 * b4) / (a3 * a4)',),
+    }
+
+    # Shared variables and parameters
+    shared_variables_dict = {
+        'V': nan,
+
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+
+        'p9': 5.15e-3,
+        'p10': 0.03158,
+        'p11': 2.26e-4,
+        'p12': 0.06990,
+
+        'p13': 3.45e-5,
+        'p14': 0.05462,
+        'p15': 0.08730,
+        'p16': 8.91e-3,
+
+        'p17': 5.15e-3,
+        'p18': 0.03158,
+
+        'p19': 0.15240,
+    }
+
+    # Current expression
+    auxiliary_expression = 'p19 * state_O * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {
+        'E_Kr': -85,
+    }
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters,
+        )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model19.py b/markov_builder/models/thirty_models/model19.py
new file mode 100644
index 0000000..f0e33c4
--- /dev/null
+++ b/markov_builder/models/thirty_models/model19.py
@@ -0,0 +1,61 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_19(MarkovChain):
+    description = "IKr Markov model 19 (Sanmitra Ghosh)"
+
+    # States
+    states = ('C', 'O_1', 'I', 'O_2')
+
+    # Transitions derived directly from dot() equations
+    rates = [
+        ('C', 'O_1', 'am', 'bm'),
+        ('I', 'O_2', 'bh', 'ah')
+    ]
+
+    # Rate expressions (verbatim from .mmt)
+    rate_dictionary = {
+        'am': ('p1 * exp(p2 * V)',),
+        'bm': ('p3 * exp(-p4 * V)',),
+        'ah': ('p5 * exp(p6 * V)',),
+        'bh': ('p7 * exp(-p8 * V)',),
+    }
+
+    # Shared variables and parameters
+    shared_variables_dict = {
+        'V': nan,
+
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+
+        'p9': 0.15240,
+    }
+
+    # Current expression
+    auxiliary_expression = 'p9 * state_O_1**3 * state_O_2 * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {
+        'E_Kr': -85,
+    }
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            # open_state=self.open_state,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters,
+        )
diff --git a/markov_builder/models/thirty_models/model2.py b/markov_builder/models/thirty_models/model2.py
new file mode 100644
index 0000000..c76c605
--- /dev/null
+++ b/markov_builder/models/thirty_models/model2.py
@@ -0,0 +1,44 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_02(MarkovChain):
+    description = ""
+    states = ('O', 'C', 'I')
+    rates = [('O', 'C', 'bm', 'am'),
+             ('O', 'I', 'ah', 'bh')]
+
+    rate_dictionary = {'am': ('p1 * exp(p2 * V)',),
+                       'bm': ('p3 * exp(-p4 * V)',),
+                       'ah': ('p5 * exp(p6 * V)',),
+                       'bh': ('p7 * exp(-p8 * V)',)
+                       }
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'g_Kr': 0.1524,
+                             }
+
+    auxiliary_expression = "g_Kr * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model20.py b/markov_builder/models/thirty_models/model20.py
new file mode 100644
index 0000000..2e9f436
--- /dev/null
+++ b/markov_builder/models/thirty_models/model20.py
@@ -0,0 +1,62 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_20(MarkovChain):
+    description = ""
+    states = ('O', 'O2', 'C2', 'C1', 'C3', 'I')
+    rates = [
+        ('C1', 'O', 'a3', 'b3'),
+        ('C2', 'C1', 'a2', 'b2'),
+        ('C3', 'C2', 'a1', 'b1'),
+        ('O2', 'I', 'bh', 'ah'),
+    ]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 0.08730,
+                             'p10': 8.91e-3,
+                             'p11': 5.15e-3,
+                             'p12': 0.03158,
+                             'p13': 0.08730,
+                             'p14': 8.91e-3,
+                             'p15': 5.15e-3,
+                             'p16': 0.03158,
+                             'p17': 0.15240,
+                             }
+
+    rate_dictionary = {
+        'a1': ('p1 * exp(p2*V)',),
+        'b1': ('p3 * exp(-p4*V)',),
+        'bh': ('p5 * exp(p6*V)',),
+        'ah': ('p7 * exp(-p8*V)',),
+        'a2': ('p9 * exp(p10*V)',),
+        'b2': ('p11 * exp(-p12*V)',),
+        'a3': ('p13 * exp(p14*V)',),
+        'b3': ('p15 * exp(-p16*V)',)
+    }
+
+    auxiliary_expression = "p17 * state_O * state_O2 * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -85}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model21.py b/markov_builder/models/thirty_models/model21.py
new file mode 100644
index 0000000..f016073
--- /dev/null
+++ b/markov_builder/models/thirty_models/model21.py
@@ -0,0 +1,82 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_21(MarkovChain):
+    description = ""
+    states = ('O', 'C1', 'C2', 'C3', 'I', 'IC1', 'IC2', 'IC3')
+    rates = [
+        ('C1', 'O', 'a3', 'b3'),
+        ('C2', 'C1', 'a2', 'b2'),
+        ('C3', 'C2', 'a1', 'b1'),
+
+        ('IC1', 'I', 'a6', 'b6'),
+        ('IC2', 'IC1', 'a5', 'b5'),
+        ('IC3', 'IC2', 'a4', 'b4'),
+
+        ('C3', 'IC3', 'ah', 'bh'),
+        ('C2', 'IC2', 'ah', 'bh'),
+        ('C1', 'IC1', 'ah', 'bh'),
+        ('O', 'I', 'ah', 'bh'),
+    ]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 2.26e-4,
+                             'p10': 0.0699,
+                             'p11': 3.45e-5,
+                             'p12': 0.05462,
+                             'p13': 0.0873,
+                             'p14': 8.91e-3,
+                             'p15': 5.15e-3,
+                             'p16': 0.03158,
+                             'p17': 0.08730,
+                             'p18': 8.91e-3,
+                             'p19': 5.15e-3,
+                             'p20': 0.03158,
+                             'p21': 5.15e-3,
+                             'p22': 0.03158,
+                             'p23': 0.15240,
+                             }
+
+    rate_dictionary = {
+        'ah': ('p1 * exp(p2*V)',),
+        'bh': ('p3 * exp(-p4*V)',),
+        'a1': ('p5 * exp(p6*V)',),
+        'b1': ('p7 * exp(-p8*V)',),
+        'a2': ('p9 * exp(p10*V)',),
+        'b2': ('p11 * exp(-p12*V)',),
+        'a3': ('p13 * exp(p14*V)',),
+        'b3': ('p15 * exp(-p16*V)',),
+        'a4': ('p17 * exp(p18*V)',),
+        'b4': ('(a4*b1/a1)',),
+        'a5': ('p19 * exp(p20*V)',),
+        'b5': ('(a5*b2/a2)',),
+        'a6': ('p21 * exp(p22*V)',),
+        'b6': ('(a6*b3/a3)',)
+    }
+
+    auxiliary_expression = "p23 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -85}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model22.py b/markov_builder/models/thirty_models/model22.py
new file mode 100644
index 0000000..ed8a2b6
--- /dev/null
+++ b/markov_builder/models/thirty_models/model22.py
@@ -0,0 +1,90 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_22(MarkovChain):
+    description = "IKr Markov model 22 (Sanmitra Ghosh)"
+
+    # States
+    states = (
+        'C3', 'C2', 'C1', 'O',
+        'I', 'IC1', 'IC2', 'IC3'
+    )
+
+    rates = [
+        ('C3', 'C2', '3*am', 'bm'),
+        ('C2', 'C1', '2*am', '2*bm'),
+        ('C1', 'O', 'am', '3*bm'),
+
+        ('O', 'I', 'a1', 'b1'),
+        ('C1', 'IC1', 'a2', 'b2'),
+        ('C2', 'IC2', 'a3', 'b3'),
+        ('C3', 'IC3', 'a4', 'b4'),
+
+        ('IC1', 'I', 'am', '3*bm'),
+        ('IC2', 'IC1', '2*am', '2*bm'),
+        ('IC3', 'IC2', '3*am', 'bm'),
+    ]
+
+    # Rate expressions
+    rate_dictionary = {
+        'am': ('p1 * exp(p2 * V)',),
+        'bm': ('p3 * exp(-p4 * V)',),
+
+        'a1': ('p5 * exp(p6 * V)',),
+        'b1': ('p7 * exp(-p8 * V)',),
+
+        'a2': ('p9 * exp(p10 * V)',),
+        'b2': ('(a2 * b1) / a1',),
+
+        'a3': ('p11 * exp(p12 * V)',),
+        'b3': ('(a3 * b2) / a2',),
+
+        'a4': ('p13 * exp(p14 * V)',),
+        'b4': ('(a4 * b3) / a3',),
+    }
+
+    shared_variables_dict = {
+        'V': nan,
+
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+
+        'p9': 2.26e-4,
+        'p10': 0.06990,
+        'p11': 3.45e-5,
+        'p12': 0.05462,
+
+        'p13': 0.08730,
+        'p14': 8.91e-3,
+
+        'p15': 0.15240,
+    }
+
+    auxiliary_expression = 'p15 * state_O * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {
+        'E_Kr': -85,
+    }
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters,
+        )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model23.py b/markov_builder/models/thirty_models/model23.py
new file mode 100644
index 0000000..23dda56
--- /dev/null
+++ b/markov_builder/models/thirty_models/model23.py
@@ -0,0 +1,112 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_23(MarkovChain):
+    description = "IKr Markov model 23 (Sanmitra Ghosh)"
+
+    states = (
+        'C3', 'C2', 'C1', 'O',
+        'I', 'IC1', 'IC2', 'IC3'
+    )
+
+    rates = [
+        ('C3', 'C2', 'a1', 'b1'),
+        ('C2', 'C1', 'a2', 'b2'),
+        ('C1', 'O', 'a3', 'b3'),
+        ('O', 'I', 'a4', 'b4'),
+
+        ('C3', 'IC3', 'a7', 'b7'),
+        ('C2', 'IC2', 'a6', 'b6'),
+        ('C1', 'IC1', 'a5', 'b5'),
+
+        ('IC3', 'IC2', 'a8', 'b8'),
+        ('IC2', 'IC1', 'a9', 'b9'),
+        ('IC1', 'I', 'a10', 'b10'),
+    ]
+
+    rate_dictionary = {
+        'a1': ('p1 * exp(p2 * V)',),
+        'b1': ('p3 * exp(-p4 * V)',),
+        'a2': ('p5 * exp(p6 * V)',),
+        'b2': ('p7 * exp(-p8 * V)',),
+        'a3': ('p9 * exp(p10 * V)',),
+        'b3': ('p11 * exp(-p12 * V)',),
+        'a4': ('p13 * exp(p14 * V)',),
+        'b4': ('p15 * exp(-p16 * V)',),
+
+        'a8': ('p17 * exp(p18 * V)',),
+        'b8': ('p19 * exp(-p20 * V)',),
+        'a9': ('p21 * exp(p22 * V)',),
+        'b9': ('p23 * exp(-p24 * V)',),
+        'a10': ('p25 * exp(p26 * V)',),
+        'b10': ('p27 * exp(-p28 * V)',),
+
+        'a5': ('p29 * exp(p30 * V)',),
+        'b5': ('(a10 * a5 * b4 * b3) / (a3 * a4 * b10)',),
+        'a6': ('p31 * exp(p32 * V)',),
+        'b6': ('(a9 * a6 * b5 * b2) / (a2 * a5 * b9)',),
+        'a7': ('p33 * exp(p34 * V)',),
+        'b7': ('(a8 * a7 * b1 * b6) / (a1 * a6 * b8)',),
+    }
+
+    shared_variables_dict = {
+        'V': nan,
+
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+        'p9': 2.26e-4,
+        'p10': 0.06990,
+        'p11': 3.45e-5,
+        'p12': 0.05462,
+        'p13': 0.08730,
+        'p14': 8.91e-3,
+        'p15': 5.15e-3,
+        'p16': 0.03158,
+        'p17': 0.08730,
+        'p18': 8.91e-3,
+        'p19': 5.15e-3,
+        'p20': 0.03158,
+        'p21': 0.06990,
+        'p22': 3.45e-5,
+        'p23': 0.05462,
+        'p24': 0.08730,
+        'p25': 8.91e-3,
+        'p26': 5.15e-3,
+        'p27': 0.03158,
+        'p28': 0.08730,
+        'p29': 8.91e-3,
+        'p30': 5.15e-3,
+        'p31': 0.03158,
+        'p32': 5.15e-3,
+        'p33': 0.03158,
+        'p34': 0.03158,
+        'p35': 0.15240,
+    }
+
+    auxiliary_expression = 'p35 * state_O * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {
+        'E_Kr': -85,
+    }
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters,
+        )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model24.py b/markov_builder/models/thirty_models/model24.py
new file mode 100644
index 0000000..6d44328
--- /dev/null
+++ b/markov_builder/models/thirty_models/model24.py
@@ -0,0 +1,58 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_24(MarkovChain):
+    description = "IKr Markov model 24 (Sanmitra Ghosh)"
+
+    states = ('C4', 'C3', 'C2', 'C1', 'O', 'I')
+
+    rates = [
+        ('C4', 'C3', '4*am', 'bm'),
+        ('C3', 'C2', '3*am', '2*bm'),
+        ('C2', 'C1', '2*am', '3*bm'),
+        ('C1', 'O', 'am', '4*bm'),
+        ('O', 'I', 'a1', 'b1'),
+    ]
+
+    rate_dictionary = {
+        'am': ('p1 * exp(p2 * V)',),
+        'bm': ('p3 * exp(-p4 * V)',),
+        'a1': ('p5 * exp(p6 * V)',),
+        'b1': ('p7 * exp(-p8 * V)',),
+    }
+
+    shared_variables_dict = {
+        'V': nan,
+
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+        'p9': 0.15240,
+    }
+
+    auxiliary_expression = 'p9 * state_O * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {
+        'E_Kr': -85,
+    }
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters,
+        )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model25.py b/markov_builder/models/thirty_models/model25.py
new file mode 100644
index 0000000..b398703
--- /dev/null
+++ b/markov_builder/models/thirty_models/model25.py
@@ -0,0 +1,76 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_25(MarkovChain):
+    description = "IKr Markov model 25 (Sanmitra Ghosh)"
+
+    states = ('C4', 'C3', 'C2', 'C1', 'O', 'I')
+
+    rates = [
+        ('C4', 'C3', 'a1', 'b1'),
+        ('C3', 'C2', 'a2', 'b2'),
+        ('C2', 'C1', 'a3', 'b3'),
+        ('C1', 'O', 'a4', 'b4'),
+        ('O', 'I', 'a5', 'b5'),
+    ]
+
+    rate_dictionary = {
+        'a1': ('p1 * exp(p2 * V)',),
+        'b1': ('p3 * exp(-p4 * V)',),
+        'a2': ('p5 * exp(p6 * V)',),
+        'b2': ('p7 * exp(-p8 * V)',),
+        'a3': ('p9 * exp(p10 * V)',),
+        'b3': ('p11 * exp(-p12 * V)',),
+        'a4': ('p13 * exp(p14 * V)',),
+        'b4': ('p15 * exp(-p16 * V)',),
+        'a5': ('p17 * exp(p18 * V)',),
+        'b5': ('p19 * exp(-p20 * V)',),
+    }
+
+    shared_variables_dict = {
+        'V': nan,
+
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+        'p9': 2.26e-4,
+        'p10': 0.06990,
+        'p11': 3.45e-5,
+        'p12': 0.05462,
+        'p13': 0.08730,
+        'p14': 8.91e-3,
+        'p15': 5.15e-3,
+        'p16': 0.03158,
+        'p17': 5.15e-3,
+        'p18': 0.03158,
+        'p19': 5.15e-3,
+        'p20': 0.03158,
+        'p21': 0.15240,
+    }
+
+    auxiliary_expression = 'p21 * state_O * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {
+        'E_Kr': -85,
+    }
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters,
+        )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model26.py b/markov_builder/models/thirty_models/model26.py
new file mode 100644
index 0000000..acecfe6
--- /dev/null
+++ b/markov_builder/models/thirty_models/model26.py
@@ -0,0 +1,52 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_26(MarkovChain):
+    description = "IKr Markov model 26 (Sanmitra Ghosh) with two open states"
+
+    states = ('O_m', 'C_m', 'O_h', 'C_h')
+    rates = [
+        ('C_m', 'O_m', 'am', 'bm'),   # m gate
+        ('C_h', 'O_h', 'ah', 'bh')    # h gate
+    ]
+
+    shared_variables_dict = {
+        'V': nan,
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+        'p9': 0.15240,
+    }
+
+    rate_dictionary = {
+        'am': ('p1 * exp(p2*V)',),
+        'bm': ('p3 * exp(-p4*V)',),
+        'ah': ('p7 * exp(-p8*V)',),
+        'bh': ('p5 * exp(p6*V)',),
+    }
+
+    auxiliary_expression = 'p9 * state_O_m**4 * state_O_h * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -85}
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters
+        )
+
+        self.set_state_attribute('O_h', open_state=True)
+        self.set_state_attribute('O_m', open_state=True)
diff --git a/markov_builder/models/thirty_models/model27.py b/markov_builder/models/thirty_models/model27.py
new file mode 100644
index 0000000..4f8e7b0
--- /dev/null
+++ b/markov_builder/models/thirty_models/model27.py
@@ -0,0 +1,73 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_27(MarkovChain):
+    description = "IKr Markov model 27 (Sanmitra Ghosh) with Om and h gates"
+
+    states = ('C1', 'C2', 'C3', 'C4', 'Om', 'Oh', 'I')
+    rates = [
+        ('C1', 'Om', 'a1', 'b1'),
+        ('C2', 'C1', 'a2', 'b2'),
+        ('C3', 'C2', 'a3', 'b3'),
+        ('C3', 'C4', 'a4', 'b4'),
+        ('I', 'Oh', 'ah', 'bh')
+    ]
+
+    rate_dictionary = {
+        'a1': ('p1 * exp(p2*V)',),
+        'b1': ('p3 * exp(-p4*V)',),
+        'ah': ('p7 * exp(-p8*V)',),
+        'bh': ('p5 * exp(p6*V)',),
+        'a2': ('p9 * exp(p10*V)',),
+        'b2': ('p11 * exp(-p12*V)',),
+        'a3': ('p13 * exp(p14*V)',),
+        'b3': ('p15 * exp(-p16*V)',),
+        'a4': ('p17 * exp(p18*V)',),
+        'b4': ('p19 * exp(-p20*V)',),
+    }
+
+    shared_variables_dict = {
+        'V': nan,
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+        'p9': 0.08730,
+        'p10': 8.91e-3,
+        'p11': 5.15e-3,
+        'p12': 0.03158,
+        'p13': 0.08730,
+        'p14': 8.91e-3,
+        'p15': 5.15e-3,
+        'p16': 0.03158,
+        'p17': 5.15e-3,
+        'p18': 0.03158,
+        'p19': 5.15e-3,
+        'p20': 0.03158,
+        'p21': 0.15240,
+    }
+
+    auxiliary_expression = 'p21 * state_Om * state_Oh * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -85}
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters
+        )
+
+        self.set_state_attribute('Om', open_state=True)
+        self.set_state_attribute('Oh', open_state=True)
diff --git a/markov_builder/models/thirty_models/model28.py b/markov_builder/models/thirty_models/model28.py
new file mode 100644
index 0000000..c59d49e
--- /dev/null
+++ b/markov_builder/models/thirty_models/model28.py
@@ -0,0 +1,101 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_28(MarkovChain):
+    description = "Model 28: IKr Markov model from Sanmitra Ghosh"
+
+    # States
+    states = ('C1', 'C2', 'C3', 'C4', 'O', 'I', 'IC1', 'IC2', 'IC3', 'IC4')
+
+    # Transitions: (from, to, forward_rate, backward_rate)
+    rates = [
+        ('C1', 'O', 'a3', 'b3'),
+        ('C2', 'C1', 'a2', 'b2'),
+        ('C3', 'C2', 'a1', 'b1'),
+        ('C4', 'C3', 'a7', 'b7'),
+        ('IC1', 'I', 'a6', 'b6'),
+        ('IC2', 'IC1', 'a5', 'b5'),
+        ('IC3', 'IC2', 'a4', 'b4'),
+        ('IC4', 'IC3', 'a8', 'b8'),
+        ('C1', 'IC1', 'ah', 'bh'),
+        ('C2', 'IC2', 'ah', 'bh'),
+        ('C3', 'IC3', 'ah', 'bh'),
+        ('C4', 'IC4', 'ah', 'bh'),
+        ('O', 'I', 'ah', 'bh'),
+    ]
+
+    # Shared variables and parameters
+    shared_variables_dict = {
+        'V': nan,
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+        'p9': 2.26e-4,
+        'p10': 0.06990,
+        'p11': 3.45e-5,
+        'p12': 0.05462,
+        'p13': 0.08730,
+        'p14': 8.91e-3,
+        'p15': 5.15e-3,
+        'p16': 0.03158,
+        'p17': 0.08730,
+        'p18': 8.91e-3,
+        'p19': 5.15e-3,
+        'p20': 0.03158,
+        'p21': 5.15e-3,
+        'p22': 0.03158,
+        'p23': 5.15e-3,
+        'p24': 0.03158,
+        'p25': 5.15e-3,
+        'p26': 0.03158,
+        'p27': 5.15e-3,
+        'p28': 0.03158,
+        'p29': 0.15240,
+    }
+
+    # Rate expressions
+    rate_dictionary = {
+        'a1': ('p5 * exp(p6 * V)',),
+        'b1': ('p7 * exp(-p8 * V)',),
+        'a2': ('p9 * exp(p10 * V)',),
+        'b2': ('p11 * exp(-p12 * V)',),
+        'a3': ('p13 * exp(p14 * V)',),
+        'b3': ('p15 * exp(-p16 * V)',),
+        'a4': ('p17 * exp(p18 * V)',),
+        'b4': ('(a4 * b1) / a1',),
+        'a5': ('p19 * exp(p20 * V)',),
+        'b5': ('(a5 * b2) / a2',),
+        'a6': ('p21 * exp(p22 * V)',),
+        'b6': ('(a6 * b3) / a3',),
+        'a7': ('p23 * exp(p24 * V)',),
+        'b7': ('p25 * exp(-p26 * V)',),
+        'a8': ('p27 * exp(p28 * V)',),
+        'b8': ('(a8 * b7) / a7',),
+        'ah': ('p1 * exp(p2 * V)',),
+        'bh': ('p3 * exp(-p4 * V)',),
+    }
+
+    # Auxiliary (current) expression
+    auxiliary_expression = "p29 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+    auxiliary_parameters = {'E_Kr': -85}
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters,
+        )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model29.py b/markov_builder/models/thirty_models/model29.py
new file mode 100644
index 0000000..1bc1d35
--- /dev/null
+++ b/markov_builder/models/thirty_models/model29.py
@@ -0,0 +1,83 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_29(MarkovChain):
+    description = "IKr Markov model 29 (Sanmitra Ghosh, HH-derived, no duplicate reactions)"
+
+    # States
+    states = ('C1', 'C2', 'C3', 'C4', 'O', 'I', 'IC1', 'IC2', 'IC3', 'IC4')
+
+    # Transitions with HH-derived multiplicities (forward/backward in one tuple)
+    rates = [
+        ('C4', 'C3', '4*am', 'bm'),
+        ('C3', 'C2', '3*am', '2*bm'),
+        ('C2', 'C1', '2*am', '3*bm'),
+        ('C1', 'O', 'am', '4*bm'),
+        ('O', 'I', 'a1', 'b1'),
+        ('C1', 'IC1', 'a2', 'b2'),
+        ('C2', 'IC2', 'a3', 'b3'),
+        ('C3', 'IC3', 'a4', 'b4'),
+        ('C4', 'IC4', 'a5', 'b5'),
+        ('I', 'IC1', '4*bm', 'am'),
+        ('IC1', 'IC2', '3*bm', '2*am'),
+        ('IC2', 'IC3', '2*bm', '3*am'),
+        ('IC3', 'IC4', 'bm', '4*am'),
+    ]
+
+    # Rate expressions
+    rate_dictionary = {
+        'am': ('p1 * exp(p2 * V)',),
+        'bm': ('p3 * exp(-p4 * V)',),
+        'a1': ('p5 * exp(p6 * V)',),
+        'b1': ('p7 * exp(-p8 * V)',),
+        'a2': ('p9 * exp(p10 * V)',),
+        'b2': ('(a2 * b1) / a1',),
+        'a3': ('p11 * exp(p12 * V)',),
+        'b3': ('(a3 * b2) / a2',),
+        'a4': ('p13 * exp(p14 * V)',),
+        'b4': ('(a4 * b3) / a3',),
+        'a5': ('p15 * exp(p16 * V)',),
+        'b5': ('(a5 * b4) / a4',),
+    }
+
+    # Shared variables / parameters
+    shared_variables_dict = {
+        'V': nan,
+        'p1': 2.26e-4,
+        'p2': 0.06990,
+        'p3': 3.45e-5,
+        'p4': 0.05462,
+        'p5': 0.08730,
+        'p6': 8.91e-3,
+        'p7': 5.15e-3,
+        'p8': 0.03158,
+        'p9': 2.26e-4,
+        'p10': 0.06990,
+        'p11': 3.45e-5,
+        'p12': 0.05462,
+        'p13': 0.08730,
+        'p14': 8.91e-3,
+        'p15': 0.08730,
+        'p16': 8.91e-3,
+        'p17': 0.15240,
+    }
+
+    # Auxiliary current
+    auxiliary_expression = 'p17 * state_O * (V - E_Kr)'
+    auxiliary_symbol = 'I_Kr'
+    auxiliary_parameters = {'E_Kr': -85}
+
+    def __init__(self):
+        super().__init__(
+            states=self.states,
+            transition_rates=self.rates,
+            rate_expressions=self.rate_dictionary,
+            auxiliary_expression=self.auxiliary_expression,
+            auxiliary_symbol=self.auxiliary_symbol,
+            shared_variables=self.shared_variables_dict,
+            auxiliary_parameters=self.auxiliary_parameters
+        )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model3.py b/markov_builder/models/thirty_models/model3.py
new file mode 100644
index 0000000..f565f21
--- /dev/null
+++ b/markov_builder/models/thirty_models/model3.py
@@ -0,0 +1,45 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_03(MarkovChain):
+    description = ""
+    states = ('O1', 'C', 'O2', 'I')
+    rates = [('O1', 'C', 'bm', 'am'),
+             ('O2', 'I', 'ah', 'bh')]
+
+    rate_dictionary = {'am': ('p1 * exp(V * p2)',),
+                       'bm': ('p3 * exp(-V * p4)',),
+                       'ah': ('p5 * exp(V * p6)',),
+                       'bh': ('p7 * exp(-V * p8)',),
+                       }
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 0.15240,
+                             }
+
+    auxiliary_expression = "p9 * state_O1 * state_O2 * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O1', open_state=True)
+        self.set_state_attribute('O2', open_state=True)
diff --git a/markov_builder/models/thirty_models/model30.py b/markov_builder/models/thirty_models/model30.py
new file mode 100644
index 0000000..3120f74
--- /dev/null
+++ b/markov_builder/models/thirty_models/model30.py
@@ -0,0 +1,117 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_30(MarkovChain):
+    description = ""
+    states = ('O', 'I', 'IC1', 'IC2', 'C2', 'C1', 'C3', 'IC3', 'IC4', 'C4')
+    rates = [
+        ('C1', 'O', 'a3', 'b3'),
+        ('C4', 'C3', 'a12', 'b12'),
+        ('C3', 'C2', 'a1', 'b1'),
+        ('IC1', 'I', 'a10', 'b10'),
+        ('IC2', 'IC1', 'a9', 'b9'),
+        ('IC3', 'IC2', 'a8', 'b8'),
+        ('IC4', 'IC3', 'a11', 'b11'),
+        ('C2', 'C1', 'a2', 'b2'),
+        ('O', 'I', 'a4', 'b4'),
+        ('C2', 'IC2', 'a6', 'b6'),
+        ('C1', 'IC1', 'a5', 'b5'),
+        ('C4', 'IC4', 'a13', 'b13'),
+        ('C3', 'IC3', 'a7', 'b7')
+    ]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 2.26e-4,
+                             'p10': 0.06990,
+                             'p11': 3.45e-5,
+                             'p12': 0.05462,
+                             'p13': 0.08730,
+                             'p14': 8.91e-3,
+                             'p15': 5.15e-3,
+                             'p16': 0.03158,
+                             'p17': 0.08730,
+                             'p18': 8.91e-3,
+                             'p19': 5.15e-3,
+                             'p20': 0.03158,
+                             'p21': 0.06990,
+                             'p22': 3.45e-5,
+                             'p23': 0.05462,
+                             'p24': 0.08730,
+                             'p25': 8.91e-3,
+                             'p26': 5.15e-3,
+                             'p27': 0.03158,
+                             'p28': 0.08730,
+                             'p29': 8.91e-3,
+                             'p30': 5.15e-3,
+                             'p31': 0.03158,
+                             'p32': 5.15e-3,
+                             'p33': 0.03158,
+                             'p34': 0.03158,
+                             'p35': 8.91e-3,
+                             'p36': 5.15e-3,
+                             'p37': 0.03158,
+                             'p38': 0.08730,
+                             'p39': 8.91e-3,
+                             'p40': 5.15e-3,
+                             'p41': 0.03158,
+                             'p42': 5.15e-3,
+                             'p43': 0.03158,
+                             'p44': 0.03158,
+                             'p45': 0.15240
+                             }
+
+    rate_dictionary = {
+        'a1': ('p1 * exp( p2 * V)',),
+        'b1': ('p3 * exp(-p4 * V)',),
+        'a2': ('p5 * exp( p6 * V)',),
+        'b2': ('p7 * exp(-p8 * V)',),
+        'a3': ('p9 * exp( p10 * V)',),
+        'b3': ('p11 * exp(-p12* V)',),
+        'a4': ('p13 * exp( p14 * V)',),
+        'b4': ('p15 * exp(-p16 * V)',),
+        'a8': ('p17 * exp(p18* V)',),
+        'b8': ('p19 * exp(-p20 * V)',),
+        'a9': ('p21 * exp(p22* V)',),
+        'b9': ('p23 * exp(-p24 * V)',),
+        'a10': ('p25 * exp(p26* V)',),
+        'b10': ('p27 * exp(-p28 * V)',),
+        'a5': ('p29 * exp(p30 * V)',),
+        'b5': ('a10*a5*b4*(b3)/(a3*a4*b10)',),
+        'a6': ('p31 * exp(p32 * V)',),
+        'b6': ('a9*a6*b5*(b2)/(a2*a5*b9)',),
+        'a7': ('p33 * exp(p34 * V)',),
+        'b7': ('a8*a7*b1*(b6)/(a1*a6*b8)',),
+        'a11': ('p35 * exp(p36* V)',),
+        'b11': ('p37 * exp(-p38 * V)',),
+        'a12': ('p39 * exp(p40* V)',),
+        'b12': ('p41 * exp(-p42 * V)',),
+        'a13': ('p43 * exp(p44 * V)',),
+        'b13': ('a13*a11*b7*b12/(a12*a7*b11)',),
+    }
+
+    auxiliary_expression = "p45 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -85}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model4.py b/markov_builder/models/thirty_models/model4.py
new file mode 100644
index 0000000..9c0d580
--- /dev/null
+++ b/markov_builder/models/thirty_models/model4.py
@@ -0,0 +1,50 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_04(MarkovChain):
+    description = ""
+    states = ('O', 'C', 'I', 'IC')
+    rates = [('O', 'C', 'b1', 'a1'),
+             ('I', 'IC', 'a3', 'b3'),
+             ('O', 'I', 'a4', 'b4'),
+             ('C', 'IC', 'a4', 'b4')]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 5.15e-3,
+                             'p10': 0.03158,
+                             'p11': 0.15240,
+                             }
+
+    rate_dictionary = {'a1': ('p1 * exp(p2*V)',),
+                       'b1': ('p3 * exp(-p4*V)',),
+                       'a4': ('p5 * exp(p6*V)',),
+                       'b4': ('p7 * exp(-p8*V)',),
+                       'a3': ('p9 * exp(-p10*V)',),
+                       'b3': ('(a3*a1)/b1',)
+                       }
+
+    auxiliary_expression = "p11 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model5.py b/markov_builder/models/thirty_models/model5.py
new file mode 100644
index 0000000..160773a
--- /dev/null
+++ b/markov_builder/models/thirty_models/model5.py
@@ -0,0 +1,51 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_05(MarkovChain):
+    description = ""
+    states = ('O', 'C', 'I', 'IC')
+    rates = [('O', 'C', 'bm', 'am'),
+             ('I', 'IC', 'bm', 'am'),
+             ('O', 'I', 'b1', 'a1'),
+             ('C', 'IC', 'b2', 'a2')]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 5.15e-3,
+                             'p10': 0.03158,
+                             'p11': 0.15240,
+                             }
+
+    rate_dictionary = {'am': ('p1 * exp(p2*V)',),
+                       'bm': ('p3 * exp(-p4*V)',),
+                       'b1': ('p5 * exp(p6*V)',),
+                       'a1': ('p7 * exp(-p8*V)',),
+                       'b2': ('p9 * exp(p10*V)',),
+                       'a2': ('(b2*a1)/b1',)
+                       }
+
+    auxiliary_expression = "p11 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model6.py b/markov_builder/models/thirty_models/model6.py
new file mode 100644
index 0000000..cc9e96b
--- /dev/null
+++ b/markov_builder/models/thirty_models/model6.py
@@ -0,0 +1,58 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_06(MarkovChain):
+    description = ""
+    states = ('O', 'C', 'I', 'IC')
+    rates = [('O', 'C', 'a1', 'b1'),
+             ('O', 'I', 'a2', 'b2'),
+             ('I', 'IC', 'a3', 'b3'),
+             ('C', 'IC', 'a4', 'b4')]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 5.15e-3,
+                             'p10': 0.03158,
+                             'p11': 5.15e-3,
+                             'p12': 0.03158,
+                             'p13': 5.15e-3,
+                             'p14': 0.03158,
+                             'p15': 0.15240
+                             }
+
+    rate_dictionary = {
+        'b1': ('p1 * exp(p2*V)',),
+        'a1': ('p3 * exp(-p4*V)',),
+        'a2': ('p5 * exp(p6*V)',),
+        'b2': ('p7 * exp(-p8*V)',),
+        'b3': ('p9 * exp(p10*V)',),
+        'a3': ('p11 * exp(-p12*V)',),
+        'a4': ('p13 * exp(p14*V)',),
+        'b4': ('(b3*b2*a1*a4)/(b1*a2*a3)',)
+    }
+
+    auxiliary_expression = "p15 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model7.py b/markov_builder/models/thirty_models/model7.py
new file mode 100644
index 0000000..5a05936
--- /dev/null
+++ b/markov_builder/models/thirty_models/model7.py
@@ -0,0 +1,49 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_07(MarkovChain):
+    description = ""
+    states = ('O', 'C1', 'C2', 'I')
+    rates = [
+        ('O', 'C1', 'bm*2', 'am'),
+        ('C1', 'C2', 'bm', 'am*2'),
+        ('O', 'I', 'a1', 'b1')
+    ]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 0.15240,
+                             }
+
+    rate_dictionary = {
+        'am': ('p1 * exp(p2*V)',),
+        'bm': ('p3 * exp(-p4*V)',),
+        'a1': ('p5 * exp(p6*V)',),
+        'b1': ('p7 * exp(-p8*V)',),
+    }
+
+    auxiliary_expression = "p9 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model8.py b/markov_builder/models/thirty_models/model8.py
new file mode 100644
index 0000000..86f65f6
--- /dev/null
+++ b/markov_builder/models/thirty_models/model8.py
@@ -0,0 +1,55 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_08(MarkovChain):
+    description = ""
+    states = ('O', 'Y1', 'Y2', 'Y4')
+    rates = [
+        ('O', 'Y2', 'k43', 'k34'),
+        ('O', 'Y4', 'k56', 'k65'),
+        ('Y1', 'Y2', 'k12', 'k21')
+    ]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 0.08730,
+                             'p10': 8.91e-3,
+                             'p11': 5.15e-3,
+                             'p12': 0.03158,
+                             'p13': 0.15240
+                             }
+
+    rate_dictionary = {
+        'k12': ('p1 * exp( p2 * V)',),
+        'k21': ('p3 * exp(-p4 * V)',),
+        'k34': ('p5 * exp( p6 * V)',),
+        'k43': ('p7 * exp(-p8 * V)',),
+        'k56': ('p9 * exp( p10 * V)',),
+        'k65': ('p11 * exp(-p12 * V)',)
+    }
+
+    auxiliary_expression = "p13 * state_O * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+
+        self.set_state_attribute('O', open_state=True)
diff --git a/markov_builder/models/thirty_models/model9.py b/markov_builder/models/thirty_models/model9.py
new file mode 100644
index 0000000..5464a71
--- /dev/null
+++ b/markov_builder/models/thirty_models/model9.py
@@ -0,0 +1,46 @@
+from numpy import nan
+
+from markov_builder.MarkovChain import MarkovChain
+
+
+class model_09(MarkovChain):
+    description = ""
+    states = ('O', 'C', 'O2', 'I')
+    rates = [('C', 'O', 'am', 'bm'),
+             ('I', 'O2', 'ah', 'bh')]
+
+    shared_variables_dict = {'V': nan,
+                             'p1': 2.26e-4,
+                             'p2': 0.06990,
+                             'p3': 3.45e-5,
+                             'p4': 0.05462,
+                             'p5': 0.08730,
+                             'p6': 8.91e-3,
+                             'p7': 5.15e-3,
+                             'p8': 0.03158,
+                             'p9': 0.15240,
+                             }
+
+    rate_dictionary = {
+        'am': ('p1 * exp(p2*V)',),
+        'bm': ('p3 * exp(-p4*V)',),
+        'bh': ('p5 * exp(p6*V)',),
+        'ah': ('p7 * exp(-p8*V)',),
+    }
+
+    auxiliary_expression = "p9 * state_O ** 2 * state_O2 * (V - E_Kr)"
+    auxiliary_symbol = 'I_Kr'
+
+    auxiliary_parameters = {'E_Kr': -88}
+
+    def __init__(self):
+        super().__init__(states=self.states,
+                         transition_rates=self.rates,
+                         rate_expressions=self.rate_dictionary,
+                         auxiliary_expression=self.auxiliary_expression,
+                         auxiliary_symbol=self.auxiliary_symbol,
+                         shared_variables=self.shared_variables_dict,
+                         auxiliary_parameters=self.auxiliary_parameters
+                         )
+        self.set_state_attribute('O', open_state=True)
+        self.set_state_attribute('O2', open_state=True)
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..f83093f
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,102 @@
+[build-system]
+
+requires = [
+         "setuptools>=61",
+         "wheel",
+         "setuptools_scm[toml]>=8.0"
+         ]
+
+build-backend = "setuptools.build_meta"
+
+
+[tool.setuptools_scm]
+local_scheme = "dirty-tag"
+write_to = "markov_builder/_version.py"
+
+[project]
+requires-python = ">=3.10"
+dynamic = ["version"]
+readme = { file = "README.md", content-type = "text/markdown" }
+
+name = "markov_builder"
+
+description = "Programmatically generate Markov models of ion-channel currents"
+
+authors = [
+            {name="Joseph Shuttleworth", email="joey.shuttleworth@nottingham.ac.uk"},
+            {name="Dominic Whittaker"},
+            {name="Michael Clerx"},
+            {name="Maurice Hendrix"},
+            {name="Gary Mirams"}
+           ]
+
+
+classifiers = [
+    'Programming Language :: Python :: 3',
+    'License :: OSI Approved :: BSD License',
+    'Operating System :: OS Independent',
+]
+
+dependencies = [
+        'scipy>=1.7',
+        'numpy>=1.21',
+        'sympy>=1.10',
+        'networkx>=3.4.2',
+        'numba>=0.61.2',
+        'pygraphviz>=1.14'
+]
+
+[project.optional-dependencies]
+test = [
+          'pytest-cov>=2.10',     # For coverage checking
+          'pytest>=4.6',          # For unit tests
+          'flake8>=3.8',            # For code style checking
+          'isort>=7.0.0',
+          'codecov>=2.1.3',
+          'matplotlib>=3.7.0',
+          'myokit>=1.38.0',
+          'seaborn>=0.13.0',
+	        'pyvis>=0.3.2'
+      ]
+
+docs = [
+	"sphinx>=9.0.0",
+	"ruff>=0.14.7"
+	]
+
+
+[tool.setuptools.packages.find]
+include = ["markov_builder"]
+
+[project.urls]
+Homepage = "https://github.com/CardiacModelling/markov-builder"
+Source = "https://github.com/CardiacModelling/markov-builder"
+
+[tool.isort]
+skip_glob = ["**/_version.py", "**.conf.py", "**/docs", ".*"]
+
+lines_after_imports = 2
+
+[tool.ruff]
+extend-select = ["D"]
+fix = true
+exclude = ["markov_builder/models/thirty_models/", "*/__init__.py"]
+
+# Ignore empty docstring on init
+ignore = ["D100", "D107"]
+
+[tool.flake8]
+max-line-length = 120
+ignore = ["W391", "W503", "W504", "E226", "E501"]
+
+exclude = [
+    ".git",
+    ".*",
+    ".venv",
+    "docs",
+    "build",
+    "doctrees",
+    "test_output",
+    "**/__pycache__",
+    "*.egg_info"
+]
\ No newline at end of file
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 30e87da..0000000
--- a/setup.py
+++ /dev/null
@@ -1,72 +0,0 @@
-#
-# markov_builder setuptools script
-#
-import os
-
-from setuptools import find_packages, setup
-
-
-# Load text for description
-with open('README.md') as f:
-    readme = f.read()
-
-# Load version number
-with open(os.path.join('markov_builder', 'version.txt'), 'r') as f:
-    version = f.read()
-
-# Go!
-setup(
-    # Module name (lowercase)
-    name='markov_builder',
-
-    version=version,
-    description='markov builder for cardiac modelling',
-    long_description=readme,
-    long_description_content_type="text/markdown",
-    author='Joseph Shuttleworth, Dominic Whittaker, Michael Clerx, Maurice Hendrix, Gary Mirams',
-    author_email='joseph.shuttleworth@nottingham.ac.uk',
-    maintainer='Maurice Hendrix',
-    maintainer_email='joseph.shuttleworth@nottingham.ac.uk',
-    url='https://github.com/CardiacModelling/markov-builder',
-    classifiers=[
-        "Programming Language :: Python :: 3",
-        "License :: OSI Approved :: BSD License",
-        "Operating System :: OS Independent",
-    ],
-
-    # Packages to include
-    packages=find_packages(
-        include=('markov_builder', 'markov_builder.*')),
-
-    # Include non-python files (via MANIFEST.in)
-    include_package_data=True,
-
-    # Required Python version
-    python_requires='>=3.6',
-
-    # List of dependencies
-    install_requires=[
-        'pints>=0.3',
-        'scipy>=1.7',
-        'numpy>=1.21',
-        'matplotlib>=3.4',
-        'pandas>=1.3',
-        'networkx>=2.6',
-        'plotly>=5.3',
-        'symengine>=0.8',
-        'sympy>=1.8',
-        'pyvis>=0.1.9',
-        'pygraphviz>=1.0',
-        'myokit>=1.33.0'
-    ],
-    extras_require={
-        'test': [
-            'pytest-cov>=2.10',     # For coverage checking
-            'pytest>=4.6',          # For unit tests
-            'flake8>=3',            # For code style checking
-            'isort',
-            'mock>=3.0.5',         # For mocking command line args etc.
-            'codecov>=2.1.3',
-        ],
-    },
-)
diff --git a/tests/models_myokit/model-0.mmt b/tests/models_myokit/model-0.mmt
new file mode 100644
index 0000000..e717a94
--- /dev/null
+++ b/tests/models_myokit/model-0.mmt
@@ -0,0 +1,54 @@
+[[model]]
+name: Beattie-2017-IKr
+author: Michael Clerx
+desc: Kylie's model as imported from the Mex file
+# Initial values
+ikr.y1 = 0
+ikr.y2 = 0
+ikr.y3 = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+# Model from Kylie's mex file (HH)
+#
+[ikr]
+use membrane.V
+IKr = p9 * y3 * (V - nernst.EK)
+y4 = 1 - y1 - y2 - y3
+k12 = p1 * exp( p2 * V)
+k21 = p3 * exp(-p4 * V)
+k41 = p5 * exp( p6 * V)
+k14 = p7 * exp(-p8 * V)
+dot(y1) = -(k12 + k14) * y1 + k21 * y2 + k41 * y4 
+dot(y2) = -(k14 + k21) * y2 + k41 * y3 + k12 * y1 
+dot(y3) = -(k21 + k41) * y3 + k12 * y4 + k14 * y2 
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.44e-5 [1/ms]
+p4 = 0.0546 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.15240 [mS/uF]
+
+n_params = 9
+
diff --git a/tests/models_myokit/model-1.mmt b/tests/models_myokit/model-1.mmt
new file mode 100644
index 0000000..39d1ce6
--- /dev/null
+++ b/tests/models_myokit/model-1.mmt
@@ -0,0 +1,44 @@
+[[model]]
+name: model-1
+author: Sanmitra Ghosh
+desc: Check associated model deifintion document
+# Initial values
+ikr.O = 0
+
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+# Model from Kylie's mex file (HH)
+#
+[ikr]
+use membrane.V
+IKr = p5 * O * (V - nernst.EK)
+k12 = p1 * exp( p2 * V)
+k21 = p3 * exp(-p4 * V)
+dot(O) = k12 * (1-O) - k21*O
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+
+p5 = 0.15240 [mS/uF]
+
+n_params = 5
\ No newline at end of file
diff --git a/tests/models_myokit/model-10.mmt b/tests/models_myokit/model-10.mmt
new file mode 100644
index 0000000..bdcb019
--- /dev/null
+++ b/tests/models_myokit/model-10.mmt
@@ -0,0 +1,64 @@
+[[model]]
+name: model-10
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C1 = 0
+ikr.Om = 0
+ikr.h = 0#5.202995292628385e-06
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p13 * O * (V - nernst.EK)
+C2 = 1 -(Om+C1)
+O = Om*h
+a1 = p1 * exp(p2*V)
+b1 = p3 * exp(-p4*V)
+bh = p5 * exp(p6*V)
+ah = p7 * exp(-p8*V)
+a2 = p9 * exp(p10*V)
+b2 = p11 * exp(-p12*V)
+
+h_inf = ah/(ah + bh)
+tauh = 1/(ah + bh)
+
+dot(C1) = b2*Om + a1*C2 - C1*(a2 +b1)
+dot(Om) = a2*C1 - b2*Om
+dot(h) = (h_inf - h)/tauh
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.08730 [1/ms]
+p10 = 8.91e-3 [1/mV]
+p11 = 5.15e-3 [1/ms]
+p12 = 0.03158 [1/mV]
+p13 = 0.15240 [mS/uF]
+
+n_params = 13
+
diff --git a/tests/models_myokit/model-11.mmt b/tests/models_myokit/model-11.mmt
new file mode 100644
index 0000000..4f6a92a
--- /dev/null
+++ b/tests/models_myokit/model-11.mmt
@@ -0,0 +1,82 @@
+[[model]]
+name: model-11
+author: Sanmitra Ghosh
+desc: Check associated model deifintion document
+# Initial values
+ikr.C2 = 1
+ikr.C1 = 0
+ikr.O = 0
+ikr.I = 0
+ikr.IC1 = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+# MODEL
+#          a4       a3
+#     IC2 <==> IC1  <==>  I
+#     ||   b4*  ||   b3* || 
+#   bh||ah    bh||ah   bh||ah
+#     ||   a2   ||   a1  || 
+#     C2  <==>  C1  <==>  O
+#          b2        b1
+#
+[ikr]
+use membrane.V
+IKr = p17 * O * (V - nernst.EK)
+
+IC2 = 1 - (C1 + C2 + O + I + IC1)
+
+a1 = p1 * exp( p2 * V)
+b1 = p3 * exp(-p4 * V)
+bh = p5 * exp( p6 * V)
+ah = p7 * exp(-p8 * V)
+a2 = p9 * exp( p10 * V)
+b2 = p11 * exp(-p12 * V)
+a3 = p13 * exp( p14 * V)
+b3 = (a3*b1)/a1
+a4 = p15 * exp(p16* V)
+b4 = (a4*b2)/a2
+
+dot(C2) = ah*IC2 + b2*C1 - C2*(a2 + bh)
+dot(C1) = ah*IC1 + b1*O + a2*C2 -C1*(a1+bh+b2)
+dot(O) = ah*I + a1*C1 - O*(bh+b1)
+dot(I) = bh*O + a3* IC1 - I*(ah+b3)
+dot(IC1) = b3*I + a4*IC2 + bh*C1 - IC1*(a3 + ah + b4)
+
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 5.15e-3 [1/ms]
+p16 = 0.03158 [1/mV]
+p17 = 0.15240 [mS/uF]
+
+n_params = 17
+
diff --git a/tests/models_myokit/model-12.mmt b/tests/models_myokit/model-12.mmt
new file mode 100644
index 0000000..4051862
--- /dev/null
+++ b/tests/models_myokit/model-12.mmt
@@ -0,0 +1,77 @@
+[[model]]
+name: model-12
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C2 = 1
+ikr.C1 = 0
+ikr.O = 0
+ikr.I = 0
+ikr.IC1 = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+# MODEL
+#          2am       am
+#     IC2 <==> IC1  <==>  I
+#     ||   bm   ||   2bm || 
+#  b3*||a3   b2*||a2   b1||a1
+#     ||   2am  ||   am  || 
+#     C2  <==>  C1  <==>  O
+#          bm        2bm
+#
+[ikr]
+use membrane.V
+IKr = p13 * O * (V - nernst.EK)
+
+IC2 = 1 - (C1 + C2 + O + I + IC1)
+
+am = p1 * exp( p2 * V)
+bm = p3 * exp(-p4 * V)
+a1 = p5 * exp( p6 * V)
+b1 = p7 * exp(-p8 * V)
+a2 = p9 * exp( p10 * V)
+a3 = p11 * exp(p12* V)
+
+b2 = (a2*b1)/a1
+b3 = (a3*b2)/a2
+
+dot(C2) = b3*IC2 + bm*C1 - C2*(a3 + 2*am)
+dot(C1) = b2*IC1 + 2*bm*O + 2*am*C2 -C1*(am + bm + a2)
+dot(O) = b1*I + am*C1 -O*(a1+2*bm)
+dot(I) = a1*O + am* IC1 - I*(b1 + 2*bm)
+dot(IC1) = 2*bm*I + 2*am*IC2 + a2*C1 - IC1*(am + bm + b2)
+
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 0.15240 [mS/uF]
+
+n_params = 13
+
diff --git a/tests/models_myokit/model-13.mmt b/tests/models_myokit/model-13.mmt
new file mode 100644
index 0000000..08d125d
--- /dev/null
+++ b/tests/models_myokit/model-13.mmt
@@ -0,0 +1,88 @@
+[[model]]
+name: model-13
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C2 = 1
+ikr.C1 = 0
+ikr.O = 0
+ikr.I = 0
+ikr.IC1 = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+
+#
+[ikr]
+use membrane.V
+IKr = p25 * O * (V - nernst.EK)
+
+IC2 = 1 - (C1 + C2 + O + I + IC1)
+
+a1 = p1 * exp( p2 * V)
+b1 = p3 * exp(-p4 * V)
+a2 = p5 * exp( p6 * V)
+b2 = p7 * exp(-p8 * V)
+a3 = p9 * exp( p10 * V)
+b3 = p11 * exp(-p12* V)
+a4 = p13 * exp(p14*V)
+b4 = (a7*b3*b2*a4)/(a2*a3*b7)
+a5 = p15 * exp(p16*V)
+b5 = (a5*a6*b4*b1)/(a1*a4*b6)
+a6 = p17 * exp( p18 * V)
+b6 = p19 * exp( -p20 * V)
+a7 = p21 * exp( p22 * V)
+b7 = p23 * exp( -p24 * V)
+
+dot(C2) = b5*IC2 + b1*C1 - C2*(a1 + a5)
+dot(C1) = b4*IC1 + b2*O + a1*C2 -C1*(a4 + b1 + a2)
+dot(O) = b3*I + a2*C1 - O*(a3+b2)
+dot(I) = a3*O + a7* IC1 - I*(b3 + b7)
+dot(IC1) = b7*I + a6*IC2 + a4*C1 - IC1*(a7 + b4 + b6)
+
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 2.26e-4 [1/ms]
+p14 = 0.06990 [1/mV]
+p15 = 3.45e-5 [1/ms]
+p16 = 0.05462 [1/mV]
+p17 = 0.08730 [1/ms]
+p18 = 8.91e-3 [1/mV]
+p19 = 5.15e-3 [1/ms]
+p20 = 0.03158 [1/mV]
+p21 = 2.26e-4 [1/ms]
+p22 = 0.06990 [1/mV]
+p23 = 3.45e-5 [1/ms]
+p24 = 0.05462 [1/mV]
+p25 = 0.15240 [mS/uF]
+
+n_params = 25
+
+
diff --git a/tests/models_myokit/model-14.mmt b/tests/models_myokit/model-14.mmt
new file mode 100644
index 0000000..f23ce83
--- /dev/null
+++ b/tests/models_myokit/model-14.mmt
@@ -0,0 +1,57 @@
+[[model]]
+name: model-14
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C3 = 1
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p9 * O * (V - nernst.EK)
+I = 1 - (O + C1 + C2 + C3)
+
+am = p1 * exp( p2 * V)
+bm = p3 * exp(-p4 * V)
+b1 = p5 * exp( p6 * V)
+a1 = p7 * exp(-p8 * V)
+
+
+dot(C1) = 3*bm  * O + 2*am * C2 -C1*(am + 2*bm)
+dot(C2) = 2*bm * C1 + 3*am * C3 - C2*(2*am + bm)
+dot(C3) = bm * C2 - 3*am * C3
+dot(O) = a1 * I + am * C1 - O*(b1 + 3*bm)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.15240 [mS/uF]
+
+n_params = 9
diff --git a/tests/models_myokit/model-15.mmt b/tests/models_myokit/model-15.mmt
new file mode 100644
index 0000000..3b121f6
--- /dev/null
+++ b/tests/models_myokit/model-15.mmt
@@ -0,0 +1,69 @@
+[[model]]
+name: model-15
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C3 = 1
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p17 * O * (V - nernst.EK)
+I = 1 - (O + C1 + C2 + C3)
+
+a1 = p1 * exp( p2 * V)
+b1 = p3 * exp(-p4 * V)
+a2 = p5 * exp( p6 * V)
+b2 = p7 * exp(-p8 * V)
+a3 = p9 * exp( p10 * V)
+b3 = p11 * exp(-p12 * V)
+a4 = p13 * exp( p14 * V)
+b4 = p15 * exp(-p16 * V)
+
+
+dot(C1) = b3  * O + a2 * C2 -C1*(a3 + b2)
+dot(C2) = b2 * C1 + a1 * C3 - C2*(a2 + b1)
+dot(C3) = b1 * C2 - a1 * C3
+dot(O) = b4 * I + a3 * C1 - O*(a4 + b3)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 5.15e-3 [1/ms]
+p16 = 0.03158 [1/mV]
+p17 = 0.15240 [mS/uF]
+
+n_params = 17
diff --git a/tests/models_myokit/model-16.mmt b/tests/models_myokit/model-16.mmt
new file mode 100644
index 0000000..2489967
--- /dev/null
+++ b/tests/models_myokit/model-16.mmt
@@ -0,0 +1,66 @@
+[[model]]
+name: model-16
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C3 = 1
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+# Wang et al. 1997 Model from Kylie's mex file (HH)
+#
+[ikr]
+use membrane.V
+IKr = p15 * O * (V - nernst.EK)
+I = 1 - (O + C1 + C2 + C3)
+a1 = p1 * exp( p2 * V)
+b1 = p3 * exp(-p4 * V)
+a3 = p5 * exp( p6 * V)
+b3 = p7 * exp(-p8 * V)
+a4 = p9 * exp( p10 * V)
+b4 = p11 * exp(-p12 * V)
+a2 = p13 
+b2 = p14 
+
+dot(C1) = b3  * O + a2 * C2 -(a3 + b2) * C1
+dot(C2) = b2 * C1 + a1 * C3 - (a2 +b1) * C2
+dot(C3) = b1 * C2 - a1 * C3
+dot(O) = b4 * I + a3 * C1 - (a4 + b3) * O
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.08730 [1/ms]
+p10 = 8.91e-3 [1/mV]
+p11 = 5.15e-3 [1/ms]
+p12 = 0.03158 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 0.15240 [mS/uF]
+
+n_params = 15
diff --git a/tests/models_myokit/model-17.mmt b/tests/models_myokit/model-17.mmt
new file mode 100644
index 0000000..4ea03b2
--- /dev/null
+++ b/tests/models_myokit/model-17.mmt
@@ -0,0 +1,60 @@
+[[model]]
+name: model-17
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C3 = 1
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p11 * O * (V - nernst.EK)
+I = 1 - (O + C1 + C2 + C3)
+
+am = p1 * exp( p2 * V)
+bm = p3 * exp(-p4 * V)
+a1 = p5 * exp( p6 * V)
+b1 = p7 * exp(-p8 * V)
+a2 = p9 * exp(p10 * V)
+b2 = (a2*3*bm*b1)/(am*a1)
+
+dot(C3) = bm * C2 - 3*am * C3
+dot(C2) = 2*bm * C1 + 3*am * C3 - C2*(2*am + bm)
+dot(C1) = 3*bm  * O + 2*am * C2 + b2*I -C1*(am + 2*bm + a2)
+dot(O) = b1 * I + am * C1 - O*(a1 + 3*bm)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 5.15e-3 [1/ms]
+p10 = 0.03158 [1/mV]
+p11 = 0.15240 [mS/uF]
+
+n_params = 11
diff --git a/tests/models_myokit/model-18.mmt b/tests/models_myokit/model-18.mmt
new file mode 100644
index 0000000..dc12305
--- /dev/null
+++ b/tests/models_myokit/model-18.mmt
@@ -0,0 +1,72 @@
+[[model]]
+name: model-18
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C3 = 1
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p19 * O * (V - nernst.EK)
+I = 1 - (O + C1 + C2 + C3)
+
+a1 = p1 * exp( p2 * V)
+b1 = p3 * exp(-p4 * V)
+a2 = p5 * exp( p6 * V)
+b2 = p7 * exp(-p8 * V)
+a3 = p9 * exp(p10 * V)
+b3 = p11 * exp(-p12 * V)
+a4 = p13 * exp(p14 * V)
+b4 = p15 * exp(-p16 * V)
+a5 = p17 * exp(p18 * V)
+b5 = (a5*b3*b4)/(a3*a4)
+
+dot(C3) = +b1*C2 - a1*C3
+dot(C2) = +b2*C1 + a1*C3 - C2*(a2 + b1)
+dot(C1) = +b3*O  + a2*C2 + b5*I - C1*(a3 + b2 + a5)
+dot(O)  = +b4*I  + a3*C1 - O*(b3 + a4)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 5.15e-3 [1/ms]
+p10 = 0.03158 [1/mV]
+p11 = 2.26e-4 [1/ms]
+p12 = 0.06990 [1/mV]
+p13 = 3.45e-5 [1/ms]
+p14 = 0.05462 [1/mV]
+p15 = 0.08730 [1/ms]
+p16 = 8.91e-3 [1/mV]
+p17 = 5.15e-3 [1/ms]
+p18 = 0.03158 [1/mV]
+p19 = 0.15240 [mS/uF]
+
+n_params = 19
diff --git a/tests/models_myokit/model-19.mmt b/tests/models_myokit/model-19.mmt
new file mode 100644
index 0000000..6fb0c19
--- /dev/null
+++ b/tests/models_myokit/model-19.mmt
@@ -0,0 +1,60 @@
+[[model]]
+name: model-19
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.m = 0#7.421666396065612e-06
+ikr.h = 0#5.202995292628385e-06
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p9 * O * (V - nernst.EK)
+
+O = m*m*m*h
+am = p1 * exp(p2*V)
+bm = p3 * exp(-p4*V)
+bh = p5 * exp(p6*V)
+ah = p7 * exp(-p8*V)
+
+m_inf = am/(am + bm)
+taum = 1/(am + bm)
+h_inf = ah/(ah + bh)
+tauh = 1/(ah + bh)
+
+dot(m) = (m_inf - m)/taum
+dot(h) = (h_inf - h)/tauh
+
+
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.15240 [mS/uF]
+
+n_params = 9
+
diff --git a/tests/models_myokit/model-2.mmt b/tests/models_myokit/model-2.mmt
new file mode 100644
index 0000000..10bfd09
--- /dev/null
+++ b/tests/models_myokit/model-2.mmt
@@ -0,0 +1,52 @@
+[[model]]
+name: model-2
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.O = 0
+ikr.C = 1
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+
+[ikr]
+use membrane.V
+IKr = p9 * O * (V - nernst.EK)
+I = 1 - (O + C)
+am = p1 * exp( p2 * V)
+bm = p3 * exp(-p4 * V)
+b = p5 * exp( p6 * V)
+a = p7 * exp(-p8 * V)
+
+dot(O) = -(bm+b)*O + am*C + a*I
+dot(C) = -am*C + bm*O
+
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.15240 [mS/uF]
+
+n_params = 9
+
diff --git a/tests/models_myokit/model-20.mmt b/tests/models_myokit/model-20.mmt
new file mode 100644
index 0000000..87b7a20
--- /dev/null
+++ b/tests/models_myokit/model-20.mmt
@@ -0,0 +1,74 @@
+[[model]]
+name: model-20
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.Om = 0
+ikr.h = 0           #5.202995292628385e-06
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p17 * O * (V - nernst.EK)
+O = Om*h
+
+
+C3 = 1 - (Om + C1 + C2)
+a1 = p1 * exp(p2*V)
+b1 = p3 * exp(-p4*V)
+bh = p5 * exp(p6*V)
+ah = p7 * exp(-p8*V)
+a2 = p9 * exp(p10*V)
+b2 = p11 * exp(-p12*V)
+a3 = p13 * exp(p14*V)
+b3 = p15 * exp(-p16*V)
+
+h_inf = ah/(ah + bh)
+tauh = 1/(ah + bh)
+
+dot(C2) = b2*C1 + a1*C3 - C2*(a2 +b1)
+dot(C1) = b3*Om + a2*C2 - C1*(a3 +b2)
+dot(Om) = a3*C1 - b3*Om
+dot(h) = (h_inf - h)/tauh
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.08730 [1/ms]
+p10 = 8.91e-3 [1/mV]
+p11 = 5.15e-3 [1/ms]
+p12 = 0.03158 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 5.15e-3 [1/ms]
+p16 = 0.03158 [1/mV]
+p17 = 0.15240 [mS/uF]
+
+n_params = 17
+
diff --git a/tests/models_myokit/model-21.mmt b/tests/models_myokit/model-21.mmt
new file mode 100644
index 0000000..4e4c6eb
--- /dev/null
+++ b/tests/models_myokit/model-21.mmt
@@ -0,0 +1,88 @@
+[[model]]
+name: model-21
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C3 = 1
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+ikr.I = 0
+ikr.IC1 = 0
+ikr.IC2 = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p23 * O * (V - nernst.EK)
+
+IC3 = 1 - (C1 + C2 + C3 + O + I + IC1 + IC2)
+
+ah = p1 * exp( p2 * V)
+bh = p3 * exp(-p4 * V)
+a1 = p5 * exp( p6 * V)
+b1 = p7 * exp(-p8 * V)
+a2 = p9 * exp( p10 * V)
+b2 = p11 * exp(-p12* V)
+a3 = p13 * exp( p14 * V)
+b3 = p15 * exp(-p16 * V)
+a4 = p17 * exp(p18* V)
+b4 = (a4*b1)/a1
+a5 = p19 * exp(p20 * V)
+b5 = (a5*b2)/a2
+a6 = p21 * exp(p22 * V)
+b6 = (a6*b3)/a3
+
+dot(C3) = b1*C2 + bh*IC3 - C3*(a1 + ah)
+dot(C2) = b2*C1 + a1*C3 +bh*IC2 - C2*(a2 + b1 + ah)
+dot(C1) = b3*O + a2*C2 + bh*IC1 - C1*(a3 + b2 + ah)
+dot(O) = bh*I + a3*C1 - O*(b3 + ah)
+dot(I) = ah*O + a6*IC1 - I*(bh + b6)
+dot(IC1) = b6*I + a5*IC2 + ah*C1 - IC1*(a6 + bh + b5)
+dot(IC2) = b5*IC1 + a4*IC3 + ah*C2 - IC2*(a5 + bh + b4)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 5.15e-3 [1/ms]
+p16 = 0.03158 [1/mV]
+p17 = 0.08730 [1/ms]
+p18 = 8.91e-3 [1/mV]
+p19 = 5.15e-3 [1/ms]
+p20 = 0.03158 [1/mV]
+p21 = 5.15e-3 [1/ms]
+p22 = 0.03158 [1/mV]
+p23 = 0.15240 [mS/uF]
+
+n_params = 23
+
diff --git a/tests/models_myokit/model-22.mmt b/tests/models_myokit/model-22.mmt
new file mode 100644
index 0000000..4df1a60
--- /dev/null
+++ b/tests/models_myokit/model-22.mmt
@@ -0,0 +1,77 @@
+[[model]]
+name: model-22
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C3 = 1
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+ikr.I = 0
+ikr.IC1 = 0
+ikr.IC2 = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+# Model from Kylie's mex file (HH)
+#
+[ikr]
+use membrane.V
+IKr = p15 * O * (V - nernst.EK)
+
+IC3 = 1 - (C1 + C2 + C3 + O + I + IC1 + IC2)
+
+am = p1 * exp( p2 * V)
+bm = p3 * exp(-p4 * V)
+a1 = p5 * exp( p6 * V)
+b1 = p7 * exp(-p8 * V)
+a2 = p9 * exp( p10 * V)
+b2 = (a2*b1)/a1
+a3 = p11 * exp( p12* V)
+b3 = (a3*b2)/a2
+a4 = p13 * exp( p14 * V)
+b4 = (a4*b3)/a3
+
+dot(C3) = bm*C2 + b4*IC3 - C3*(a4 + 3*am)
+dot(C2) = 2*bm*C1 + 3*am*C3 + b3*IC2 - C2*(a3 + 2*am + bm)
+dot(C1) = 3*bm*O + 2*am*C2 + b2*IC1 - C1*(am + 2*bm + a2)
+dot(O) = b1*I + am*C1 - O*(3*bm + a1)
+dot(I) = a1*O + am*IC1 - I*(b1 + 3*bm)
+dot(IC1) = 3*bm*I + 2*am*IC2 + a2*C1 - IC1*(2*bm + am + b2)
+dot(IC2) = 2*bm*IC1 + 3*am*IC3 + a3*C2 - IC2*(2*am + b3 + bm)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 0.15240 [mS/uF]
+
+n_params = 15
+
diff --git a/tests/models_myokit/model-23.mmt b/tests/models_myokit/model-23.mmt
new file mode 100644
index 0000000..501ee72
--- /dev/null
+++ b/tests/models_myokit/model-23.mmt
@@ -0,0 +1,107 @@
+[[model]]
+name: model-23
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C3 = 1
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+ikr.I = 0
+ikr.IC1 = 0
+ikr.IC2 = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p35 * O * (V - nernst.EK)
+
+IC3 = 1 - (C1 + C2 + C3 + O + I + IC1 + IC2)
+
+a1 = p1 * exp( p2 * V)
+b1 = p3 * exp(-p4 * V)
+a2 = p5 * exp( p6 * V)
+b2 = p7 * exp(-p8 * V)
+a3 = p9 * exp( p10 * V)
+b3 = p11 * exp(-p12* V)
+a4 = p13 * exp( p14 * V)
+b4 = p15 * exp(-p16 * V)
+a8 = p17 * exp(p18* V)
+b8 = p19 * exp(-p20 * V)
+a9 = p21 * exp(p22* V)
+b9 = p23 * exp(-p24 * V) 
+a10 = p25 * exp(p26* V)
+b10 = p27 * exp(-p28 * V) 
+a5 = p29 * exp(p30 * V)
+b5 = (a10*a5*b4*b3)/(a3*a4*b10)
+a6 = p31 * exp(p32 * V)
+b6 = (a9*a6*b5*b2)/(a2*a5*b9)
+a7 = p33 * exp(p34 * V)
+b7 = (a8*a7*b1*b6)/(a1*a6*b8)
+
+dot(C3) = b1*C2 + b7*IC3 - C3*(a1 + a7)
+dot(C2) = b2*C1 + a1*C3 +b6*IC2 - C2*(a2 + b1 + a6)
+dot(C1) = b3*O + a2*C2 + b5*IC1 - C1*(a3 + b2 + a5)
+dot(O) = b4*I + a3*C1 - O*(b3 + a4)
+dot(I) = a4*O + a10*IC1 - I*(b4 + b10)
+dot(IC1) = b10*I + a9*IC2 + a5*C1 - IC1*(b9 + a10 + b5)
+dot(IC2) = b9*IC1 + a8*IC3 + a6*C2 - IC2*(a9 + b8 + b6)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 5.15e-3 [1/ms]
+p16 = 0.03158 [1/mV]
+p17 = 0.08730 [1/ms]
+p18 = 8.91e-3 [1/mV]
+p19 = 5.15e-3 [1/ms]
+p20 = 0.03158 [1/mV]
+p21 = 0.06990 [1/mV]
+p22 = 3.45e-5 [1/ms]
+p23 = 0.05462 [1/mV]
+p24 = 0.08730 [1/ms]
+p25 = 8.91e-3 [1/mV]
+p26 = 5.15e-3 [1/ms]
+p27 = 0.03158 [1/mV]
+p28 = 0.08730 [1/ms]
+p29 = 8.91e-3 [1/mV]
+p30 = 5.15e-3 [1/ms]
+p31 = 0.03158 [1/mV]
+p32 = 5.15e-3 [1/ms]
+p33 = 0.03158 [1/mV]
+p34 = 0.03158 [1/mV]
+p35 = 0.15240 [mS/uF]
+
+n_params = 35
+
+
diff --git a/tests/models_myokit/model-24.mmt b/tests/models_myokit/model-24.mmt
new file mode 100644
index 0000000..a194f35
--- /dev/null
+++ b/tests/models_myokit/model-24.mmt
@@ -0,0 +1,58 @@
+[[model]]
+name: model-24
+author: Sanmitra Ghosh
+desc: Check associated model definition document 
+# Initial values
+ikr.C4 = 1
+ikr.C3 = 0
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p9 * O * (V - nernst.EK)
+I = 1 - (O + C1 + C2 + C3 + C4)
+
+am = p1 * exp( p2 * V)
+bm = p3 * exp(-p4 * V)
+a1 = p5 * exp( p6 * V)
+b1 = p7 * exp(-p8 * V)
+
+dot(C4) = bm*C3 - 4*am*C4
+dot(C3) = 2*bm * C2 + 4*am * C4 - C3*(3*am + bm)
+dot(C2) = 3*bm * C1 + 3*am * C3 - C2*(2*am + 2*bm)
+dot(C1) = 4*bm  * O + 2*am * C2 -C1*(am + 3*bm)
+dot(O) = b1 * I + am * C1 - O*(a1 + 4*bm)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.15240 [mS/uF]
+
+n_params = 9
diff --git a/tests/models_myokit/model-25.mmt b/tests/models_myokit/model-25.mmt
new file mode 100644
index 0000000..8c46d4d
--- /dev/null
+++ b/tests/models_myokit/model-25.mmt
@@ -0,0 +1,76 @@
+[[model]]
+name: model-25
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C4 = 1
+ikr.C3 = 0
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p21 * O * (V - nernst.EK)
+I = 1 - (O + C1 + C2 + C3 + C4)
+
+a1 = p1 * exp( p2 * V)
+b1 = p3 * exp(-p4 * V)
+a2 = p5 * exp( p6 * V)
+b2 = p7 * exp(-p8 * V)
+a3 = p9 * exp( p10 * V)
+b3 = p11 * exp(-p12 * V)
+a4 = p13 * exp( p14 * V)
+b4 = p15 * exp(-p16 * V)
+a5 = p17 * exp( p18 * V)
+b5 = p19 * exp(-p20 * V)
+
+dot(C4) = b1 * C3 - a1 * C4
+dot(C3) = b2 * C2 + a1 * C4 - C3*(a2 + b1)
+dot(C2) = b3 * C1 + a2 * C3 - C2*(a3 + b2)
+dot(C1) = b4 * O  + a3 * C2 - C1*(a4 + b3)
+dot(O) =  b5 * I  + a4 * C1 -  O*(a5 + b4)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 5.15e-3 [1/ms]
+p16 = 0.03158 [1/mV]
+p17 = 5.15e-3 [1/ms]
+p18 = 0.03158 [1/mV]
+p19 = 5.15e-3 [1/ms]
+p20 = 0.03158 [1/mV]
+p21 = 0.15240 [mS/uF]
+
+n_params = 21
diff --git a/tests/models_myokit/model-26.mmt b/tests/models_myokit/model-26.mmt
new file mode 100644
index 0000000..d82f40d
--- /dev/null
+++ b/tests/models_myokit/model-26.mmt
@@ -0,0 +1,58 @@
+[[model]]
+name: model-26
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.m = 0 #7.421666396065612e-06
+ikr.h = 0 #5.202995292628385e-06
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p9 * O * (V - nernst.EK)
+
+O = m*m*m*m*h
+am = p1 * exp(p2*V)
+bm = p3 * exp(-p4*V)
+bh = p5 * exp(p6*V)
+ah = p7 * exp(-p8*V)
+
+m_inf = am/(am + bm)
+taum = 1/(am + bm)
+h_inf = ah/(ah + bh)
+tauh = 1/(ah + bh)
+
+dot(m) = (m_inf - m)/taum
+dot(h) = (h_inf - h)/tauh
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.15240 [mS/uF]
+
+n_params = 9
+
diff --git a/tests/models_myokit/model-27.mmt b/tests/models_myokit/model-27.mmt
new file mode 100644
index 0000000..17bbc4c
--- /dev/null
+++ b/tests/models_myokit/model-27.mmt
@@ -0,0 +1,81 @@
+[[model]]
+name: model-27
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C3 = 0
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.Om = 0
+ikr.h = 0           #5.202995292628385e-06
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p21 * O * (V - nernst.EK)
+
+O = Om*h
+C4 = 1 - (Om + C1 + C2 + C3)
+a1 = p1 * exp(p2*V)
+b1 = p3 * exp(-p4*V)
+bh = p5 * exp(p6*V)
+ah = p7 * exp(-p8*V)
+a2 = p9 * exp(p10*V)
+b2 = p11 * exp(-p12*V)
+a3 = p13 * exp(p14*V)
+b3 = p15 * exp(-p16*V)
+a4 = p17 * exp(p18*V)
+b4 = p19 * exp(-p20*V)
+
+h_inf = ah/(ah + bh)
+tauh = 1/(ah + bh)
+
+dot(C3) = b3*C2 + a4*C4 - C3*(a3 +b4)
+dot(C2) = b2*C1 + a3*C3 - C2*(a2 +b3)
+dot(C1) = b1*Om + a2*C2 - C1*(a1 +b2)
+dot(Om) = a1*C1 - b1*Om
+dot(h) = (h_inf - h)/tauh
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.08730 [1/ms]
+p10 = 8.91e-3 [1/mV]
+p11 = 5.15e-3 [1/ms]
+p12 = 0.03158 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 5.15e-3 [1/ms]
+p16 = 0.03158 [1/mV]
+p17 = 5.15e-3 [1/ms]
+p18 = 0.03158 [1/mV]
+p19 = 5.15e-3 [1/ms]
+p20 = 0.03158 [1/mV]
+p21 = 0.15240 [mS/uF]
+
+n_params = 21
+
diff --git a/tests/models_myokit/model-28.mmt b/tests/models_myokit/model-28.mmt
new file mode 100644
index 0000000..d26bb27
--- /dev/null
+++ b/tests/models_myokit/model-28.mmt
@@ -0,0 +1,102 @@
+[[model]]
+name: model-28
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C4 = 1
+ikr.C3 = 0
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+ikr.I = 0
+ikr.IC1 = 0
+ikr.IC2 = 0
+ikr.IC3 = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p29 * O * (V - nernst.EK)
+
+IC4 = 1 - (C1 + C2 + C3 + C4 + O + I + IC1 + IC2 + IC3)
+
+bh = p1 * exp( p2 * V)
+ah = p3 * exp(-p4 * V)
+a1 = p5 * exp( p6 * V)
+b1 = p7 * exp(-p8 * V)
+a2 = p9 * exp( p10 * V)
+b2 = p11 * exp(-p12* V)
+a3 = p13 * exp( p14 * V)
+b3 = p15 * exp(-p16 * V)
+a4 = p17 * exp(p18* V)
+b4 = (a4*b1)/a1
+a5 = p19 * exp(p20 * V)
+b5 = (a5*b2)/a2
+a6 = p21 * exp(p22 * V)
+b6 = (a6*b3)/a3
+a7 = p23 * exp(p24 * V)
+b7 = p25 * exp(-p26 * V)
+a8 = p27 * exp(p28 * V)
+b8 = (a8*b7)/a7
+
+dot(C4) = ah*IC4 +b7*C3 - C4*(a7+bh)
+dot(C3) = b1*C2 + ah*IC3 + a7*C4 - C3*(a1 + bh + b7)
+dot(C2) = b2*C1 + a1*C3 +ah*IC2 - C2*(a2 + b1 + bh)
+dot(C1) = b3*O + a2*C2 + ah*IC1 - C1*(a3 + b2 + bh)
+dot(O) = ah*I + a3*C1 - O*(b3 + bh)
+dot(I) = bh*O + a6*IC1 - I*(ah + b6)
+dot(IC1) = b6*I + a5*IC2 + bh*C1 - IC1*(a6 + ah + b5)
+dot(IC2) = b5*IC1 + a4*IC3 + bh*C2 - IC2*(a5 + ah + b4)
+dot(IC3) = a8*IC4 + bh*C3 + b4*IC2 - IC3*(ah + a4 + b8)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 5.15e-3 [1/ms]
+p16 = 0.03158 [1/mV]
+p17 = 0.08730 [1/ms]
+p18 = 8.91e-3 [1/mV]
+p19 = 5.15e-3 [1/ms]
+p20 = 0.03158 [1/mV]
+p21 = 5.15e-3 [1/ms]
+p22 = 0.03158 [1/mV]
+p23 = 5.15e-3 [1/ms]
+p24 = 0.03158 [1/mV]
+p25 = 5.15e-3 [1/ms]
+p26 = 0.03158 [1/mV]
+p27 = 5.15e-3 [1/ms]
+p28 = 0.03158 [1/mV]
+p29 = 0.15240 [mS/uF]
+
+n_params = 29
+
diff --git a/tests/models_myokit/model-29.mmt b/tests/models_myokit/model-29.mmt
new file mode 100644
index 0000000..51769a3
--- /dev/null
+++ b/tests/models_myokit/model-29.mmt
@@ -0,0 +1,85 @@
+[[model]]
+name: model-29
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C4 = 1
+ikr.C3 = 0
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+ikr.I = 0
+ikr.IC1 = 0
+ikr.IC2 = 0
+ikr.IC3 = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+# Model from Kylie's mex file (HH)
+#
+[ikr]
+use membrane.V
+IKr = p17 * O * (V - nernst.EK)
+
+IC4 = 1 - (C1 + C2 + C3 + C4 + O + I + IC1 + IC2 + IC3)
+
+am = p1 * exp( p2 * V)
+bm = p3 * exp(-p4 * V)
+a1 = p5 * exp( p6 * V)
+b1 = p7 * exp(-p8 * V)
+a2 = p9 * exp( p10 * V)
+b2 = (a2*b1)/a1
+a3 = p11 * exp( p12* V)
+b3 = (a3*b2)/a2
+a4 = p13 * exp( p14 * V)
+b4 = (a4*b3)/a3
+a5 = p15 * exp( p16 * V)
+b5 = (a5*b4)/a4
+
+dot(C4) = bm*C3 + b5*IC4 - C4*(4*am+a5) 
+dot(C3) = 2*bm*C2 + 4*am*C4 + b4*IC3 - C3*(a4 + 3*am + bm)
+dot(C2) = 3*bm*C1 + 3*am*C3 + b3*IC2 - C2*(a3 + 2*am + 2*bm)
+dot(C1) = 4*bm*O + 2*am*C2 + b2*IC1 - C1*(am + 3*bm + a2)
+dot(O) = b1*I + am*C1 - O*(4*bm + a1)
+dot(I) = a1*O + am*IC1 - I*(b1 + 4*bm)
+dot(IC1) = 4*bm*I + 2*am*IC2 + a2*C1 - IC1*(3*bm + am + b2)
+dot(IC2) = 3*bm*IC1 + 3*am*IC3 + a3*C2 - IC2*(2*am + b3 + 2*bm)
+dot(IC3) = 2*bm*IC2 + 4*am*IC4 + a4*C3 - IC3*(3*am + b4 + bm)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 0.08730 [1/ms]
+p16 = 8.91e-3 [1/mV]
+p17 = 0.15240 [mS/uF]
+
+n_params = 17
+
diff --git a/tests/models_myokit/model-3.mmt b/tests/models_myokit/model-3.mmt
new file mode 100644
index 0000000..fd7ead2
--- /dev/null
+++ b/tests/models_myokit/model-3.mmt
@@ -0,0 +1,57 @@
+[[model]]
+name: model-3
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.m = 0#7.421666396065612e-06
+ikr.h = 0#5.202995292628385e-06
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -88 [mV]
+
+[ikr]
+use membrane.V
+IKr = p9 * O * (V - nernst.EK)
+
+am = p1 * exp(p2*V)
+bm = p3 * exp(-p4*V)
+ah = p5 * exp(p6*V)
+bh = p7 * exp(-p8*V)
+
+m_inf = am/(am + bm)
+taum = 1/(am + bm)
+h_inf = bh/(ah + bh)
+tauh = 1/(ah + bh)
+
+O = m*h
+
+dot(m) = (m_inf - m)/taum
+dot(h) = (h_inf - h)/tauh
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.15240 [mS/uF]
+
+n_params = 9
+
diff --git a/tests/models_myokit/model-30.mmt b/tests/models_myokit/model-30.mmt
new file mode 100644
index 0000000..b7f0036
--- /dev/null
+++ b/tests/models_myokit/model-30.mmt
@@ -0,0 +1,127 @@
+[[model]]
+name: model-30
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C4 = 1
+ikr.C3 = 0
+ikr.C2 = 0
+ikr.C1 = 0
+ikr.O = 0
+ikr.I = 0
+ikr.IC1 = 0
+ikr.IC2 = 0
+ikr.IC3 = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p45 * O * (V - nernst.EK)
+
+IC4 = 1 - (C1 + C2 + C3 + C4 + O + I + IC1 + IC2 + IC3)
+
+a1 = p1 * exp( p2 * V)
+b1 = p3 * exp(-p4 * V)
+a2 = p5 * exp( p6 * V)
+b2 = p7 * exp(-p8 * V)
+a3 = p9 * exp( p10 * V)
+b3 = p11 * exp(-p12* V)
+a4 = p13 * exp( p14 * V)
+b4 = p15 * exp(-p16 * V)
+a8 = p17 * exp(p18* V)
+b8 = p19 * exp(-p20 * V)
+a9 = p21 * exp(p22* V)
+b9 = p23 * exp(-p24 * V) 
+a10 = p25 * exp(p26* V)
+b10 = p27 * exp(-p28 * V) 
+a5 = p29 * exp(p30 * V)
+b5 = (a10*a5*b4*b3)/(a3*a4*b10)
+a6 = p31 * exp(p32 * V)
+b6 = (a9*a6*b5*b2)/(a2*a5*b9)
+a7 = p33 * exp(p34 * V)
+b7 = (a8*a7*b1*b6)/(a1*a6*b8)
+a11 = p35 * exp(p36* V)
+b11 = p37 * exp(-p38 * V) 
+a12 = p39 * exp(p40* V)
+b12 = p41 * exp(-p42 * V) 
+a13 = p43 * exp(p44 * V)
+b13 = (a13*a11*b7*b12)/(a12*a7*b11)
+
+dot(C4) = b13*IC4 + b12*C3 - C4*(a13+a12)
+dot(C3) = b1*C2 + a12*C4 + b7*IC3 - C3*(a1 + a7 + b12)
+dot(C2) = b2*C1 + a1*C3  + b6*IC2 - C2*(a2 + b1 + a6)
+dot(C1) = b3*O  + a2*C2  + b5*IC1 - C1*(a3 + b2 + a5)
+dot(O) = b4*I + a3*C1 - O*(b3 + a4)
+dot(I) = a4*O + a10*IC1 - I*(b4 + b10)
+dot(IC1) = b10*I + a9*IC2 + a5*C1 - IC1*(b9 + a10 + b5)
+dot(IC2) = b9*IC1 + a8*IC3 + a6*C2 - IC2*(a9 + b8 + b6)
+dot(IC3) = b8*IC2 + a11*IC4 + a7*C3 - IC3*(a8 + b7 + b11)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 2.26e-4 [1/ms]
+p10 = 0.06990 [1/mV]
+p11 = 3.45e-5 [1/ms]
+p12 = 0.05462 [1/mV]
+p13 = 0.08730 [1/ms]
+p14 = 8.91e-3 [1/mV]
+p15 = 5.15e-3 [1/ms]
+p16 = 0.03158 [1/mV]
+p17 = 0.08730 [1/ms]
+p18 = 8.91e-3 [1/mV]
+p19 = 5.15e-3 [1/ms]
+p20 = 0.03158 [1/mV]
+p21 = 0.06990 [1/mV]
+p22 = 3.45e-5 [1/ms]
+p23 = 0.05462 [1/mV]
+p24 = 0.08730 [1/ms]
+p25 = 8.91e-3 [1/mV]
+p26 = 5.15e-3 [1/ms]
+p27 = 0.03158 [1/mV]
+p28 = 0.08730 [1/ms]
+p29 = 8.91e-3 [1/mV]
+p30 = 5.15e-3 [1/ms]
+p31 = 0.03158 [1/mV]
+p32 = 5.15e-3 [1/ms]
+p33 = 0.03158 [1/mV]
+p34 = 0.03158 [1/mV]
+p35 = 8.91e-3 [1/mV]
+p36 = 5.15e-3 [1/ms]
+p37 = 0.03158 [1/mV]
+p38 = 0.08730 [1/ms]
+p39 = 8.91e-3 [1/mV]
+p40 = 5.15e-3 [1/ms]
+p41 = 0.03158 [1/mV]
+p42 = 5.15e-3 [1/ms]
+p43 = 0.03158 [1/mV]
+p44 = 0.03158 [1/mV]
+p45 = 0.15240 [mS/uF]
+
+n_params = 45
+
+
diff --git a/tests/models_myokit/model-4.mmt b/tests/models_myokit/model-4.mmt
new file mode 100644
index 0000000..f289199
--- /dev/null
+++ b/tests/models_myokit/model-4.mmt
@@ -0,0 +1,60 @@
+[[model]]
+name: model-4
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C = 1
+ikr.O = 0
+ikr.I = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p11 * O * (V - nernst.EK)
+IC = 1 - (I + O + C)
+
+a1 = p1 * exp(p2*V)
+b1 = p3 * exp(-p4*V)
+a4 = p5 * exp(p6*V)
+b4 = p7 * exp(-p8*V)
+a3 = p9 * exp(-p10*V)
+b3 = (a3*a1)/b1
+
+
+dot(C) = b1*O + b4*IC - C*(a1+a4)
+dot(O) = a1*C + b4*I - O*(a4+b1)
+dot(I) = a4*O + b3*IC -I*(b4 + a3)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 5.15e-3 [1/ms]
+p10 = 0.03158 [1/mV]
+p11 = 0.15240 [mS/uF]
+
+n_params = 11
+
diff --git a/tests/models_myokit/model-5.mmt b/tests/models_myokit/model-5.mmt
new file mode 100644
index 0000000..9820657
--- /dev/null
+++ b/tests/models_myokit/model-5.mmt
@@ -0,0 +1,60 @@
+[[model]]
+name: model-5
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.C = 1
+ikr.O = 0
+ikr.I = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+#
+[ikr]
+use membrane.V
+IKr = p11 * O * (V - nernst.EK)
+IC = 1 - (I + O + C)
+
+am = p1 * exp(p2*V)
+bm = p3 * exp(-p4*V)
+b1 = p5 * exp(p6*V)
+a1 = p7 * exp(-p8*V)
+b2 = p9 * exp(p10*V)
+a2 = (b2*a1)/b1
+
+dot(C) = bm*O + a2*IC - C*(am+b2)
+dot(O) = a1*I + am*C - O*(bm+b1)
+dot(I) = b1*O + am*IC -I*(a1 + bm)
+
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 5.15e-3 [1/ms]
+p10 = 0.03158 [1/mV]
+p11 = 0.15240 [mS/uF]
+
+n_params = 11
+
diff --git a/tests/models_myokit/model-6.mmt b/tests/models_myokit/model-6.mmt
new file mode 100644
index 0000000..9b8c069
--- /dev/null
+++ b/tests/models_myokit/model-6.mmt
@@ -0,0 +1,66 @@
+[[model]]
+name: model-6
+author: Sanmitra Ghosh
+desc: Check associated model deifintion document
+# Initial values
+ikr.C = 1
+ikr.O = 0
+ikr.I = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membane potential
+    in [mV]
+
+[nernst]
+EK = -88 [mV]
+
+#
+# Model from Kylie's mex file (HH)
+#
+[ikr]
+use membrane.V
+IKr = p15 * O * (V - nernst.EK)
+IC = 1 - (I + O + C)
+
+b1 = p1 * exp(p2*V)
+a1 = p3 * exp(-p4*V)
+a2 = p5 * exp(p6*V)
+b2 = p7 * exp(-p8*V)
+b3 = p9 * exp(p10*V)
+a3 = p11* exp(-p12*V)
+a4 = p13* exp(p14*V)
+b4 = (b3*b2*a1*a4)/(b1*a2*a3)
+
+dot(C) = a1*O + b4*IC - C*(b1+a4)
+dot(O) = b1*C + b2*I - O*(a2+a1)
+dot(I) = a2*O + b3*IC -I*(b2 + a3)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 5.15e-3 [1/ms]
+p10 = 0.03158 [1/mV]
+p11 = 5.15e-3 [1/ms]
+p12 = 0.03158 [1/mV]
+p13 = 5.15e-3 [1/ms]
+p14 = 0.03158 [1/mV]
+p15 = 0.15240 [mS/uF]
+
+n_params = 15
+
diff --git a/tests/models_myokit/model-7.mmt b/tests/models_myokit/model-7.mmt
new file mode 100644
index 0000000..7e4e9fc
--- /dev/null
+++ b/tests/models_myokit/model-7.mmt
@@ -0,0 +1,57 @@
+[[model]]
+name: model-7
+author: Sanmitra Ghosh
+desc: Check associated model defintion document
+# Initial values
+ikr.C2 = 1
+ikr.C1 = 0
+ikr.O = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -88 [mV]
+
+#
+# Model from Kylie's mex file (HH)
+#
+[ikr]
+use membrane.V
+IKr = p9 * O * (V - nernst.EK)
+
+I = 1 - (C1 + C2 + O)
+
+am = p1 * exp( p2 * V)
+bm = p3 * exp(-p4 * V)
+a1 = p5 * exp( p6 * V)
+b1 = p7 * exp(-p8 * V)
+
+
+dot(C2) = bm*C1 - 2*am*C2
+dot(C1) = 2*am*C2 + 2*bm*O - C1*(bm+am)
+dot(O) = am*C1 + b1*I -O*(a1 + 2*bm)
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.15240 [mS/uF]
+
+n_params =9
diff --git a/tests/models_myokit/model-8.mmt b/tests/models_myokit/model-8.mmt
new file mode 100644
index 0000000..dc014b2
--- /dev/null
+++ b/tests/models_myokit/model-8.mmt
@@ -0,0 +1,61 @@
+[[model]]
+name: model-8
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.y1 = 1
+ikr.y2 = 0
+ikr.O = 0
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+# Model from Kylie's mex file (HH)
+#
+[ikr]
+use membrane.V
+IKr = p13 * O * (V - nernst.EK)
+y4 = 1 - (y1 + y2 + O)
+k12 = p1 * exp( p2 * V)
+k21 = p3 * exp(-p4 * V)
+k34 = p5 * exp( p6 * V)
+k43 = p7 * exp(-p8 * V)
+k56 = p9 * exp( p10 * V)
+k65 = p11 * exp(-p12 * V)
+
+dot(y1) = -k12 * y1 + k21 * y2
+dot(y2) = -(k21 + k34) * y2 + k12 * y1 + k43 * O
+dot(O) = -(k43 + k56) * O + k34 * y2 + k65 * y4
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.08730 [1/ms]
+p10 = 8.91e-3 [1/mV]
+p11 = 5.15e-3 [1/ms]
+p12 = 0.03158 [1/mV]
+p13 = 0.15240 [mS/uF]
+
+
+n_params = 13
diff --git a/tests/models_myokit/model-9.mmt b/tests/models_myokit/model-9.mmt
new file mode 100644
index 0000000..6d95851
--- /dev/null
+++ b/tests/models_myokit/model-9.mmt
@@ -0,0 +1,61 @@
+[[model]]
+name: model-9
+author: Sanmitra Ghosh
+desc: Check associated model definition document
+# Initial values
+ikr.m = 0#7.421666396065612e-06
+ikr.h = 0#5.202995292628385e-06
+
+#
+# Simulation engine variables
+#
+[engine]
+time = 0 bind time
+pace = 0 bind pace
+
+#
+# Membrane potential
+#
+[membrane]
+V = engine.pace
+    desc: membrane potential
+    in [mV]
+
+[nernst]
+EK = -85 [mV]
+
+#
+# Model from Kylie's mex file (HH)
+#
+[ikr]
+use membrane.V
+IKr = p9 * O * (V - nernst.EK)
+
+O = m*m*h
+am = p1 * exp(p2*V)
+bm = p3 * exp(-p4*V)
+bh = p5 * exp(p6*V)
+ah = p7 * exp(-p8*V)
+
+m_inf = am/(am + bm)
+taum = 1/(am + bm)
+h_inf = ah/(ah + bh)
+tauh = 1/(ah + bh)
+
+dot(m) = (m_inf - m)/taum
+dot(h) = (h_inf - h)/tauh
+
+
+
+p1 = 2.26e-4 [1/ms]
+p2 = 0.06990 [1/mV]
+p3 = 3.45e-5 [1/ms]
+p4 = 0.05462 [1/mV]
+p5 = 0.08730 [1/ms]
+p6 = 8.91e-3 [1/mV]
+p7 = 5.15e-3 [1/ms]
+p8 = 0.03158 [1/mV]
+p9 = 0.15240 [mS/uF]
+
+n_params = 9
+
diff --git a/tests/models_myokit/simplified-staircase.mmt b/tests/models_myokit/simplified-staircase.mmt
new file mode 100644
index 0000000..f15b398
--- /dev/null
+++ b/tests/models_myokit/simplified-staircase.mmt
@@ -0,0 +1,35 @@
+[[protocol]]
+# V     start   duration period repeats
+-80     0       250     0       0
+-120    next    50      0       0
+-120    next    400     0       0
+-80     next    200     0       0
+40      next    1000    0       0
+-120    next    500     0       0
+-80     next    1000    0       0
+-40     next    500     0       0
+-60     next    500     0       0
+-20     next    500     0       0
+-40     next    500     0       0
+0       next    500     0       0
+-20     next    500     0       0
+20      next    500     0       0
+0       next    500     0       0
+40      next    500     0       0
+20      next    500     0       0
+40      next    500     0       0
+0       next    500     0       0
+20      next    500     0       0
+-20     next    500     0       0
+0       next    500     0       0
+-40     next    500     0       0
+-20     next    500     0       0
+-60     next    500     0       0
+-40     next    500     0       0
+-80     next    1000    0       0
+40      next    500     0       0
+-70     next    10      0       0
+-70     next    100     0       0
+-120    next    390     0       0
+-80     next    500     0       0
+
diff --git a/tests/test_30_models.py b/tests/test_30_models.py
new file mode 100644
index 0000000..02a65a3
--- /dev/null
+++ b/tests/test_30_models.py
@@ -0,0 +1,182 @@
+#!/usr/bin/env python3
+
+import itertools
+import logging
+import os
+import unittest
+
+import matplotlib
+import matplotlib.pyplot as plt
+import myokit
+import myokit as mk
+import networkx as nx
+import numpy as np
+import pytest
+import sympy as sp
+
+import markov_builder
+import markov_builder.models.thirty_models
+
+
+matplotlib.use('pdf')
+
+models = [getattr(markov_builder.models.thirty_models, f"model_{i:02}") for i in range(31)]
+disconnected_models = [models[i] for i in [3, 9, 10, 19, 20, 26, 27]]
+
+output_dir = os.environ.get('MARKOVBUILDER_TEST_OUTPUT', 'test_output')
+os.makedirs(output_dir, exist_ok=True)
+if not os.path.exists(output_dir):
+    os.makedirs(output_dir)
+logging.info("outputting to " + output_dir)
+
+
+@pytest.mark.parametrize("model", models)
+def test_generate_myokit(model):
+    name = model.__name__
+    logging.debug(f"initiating {name}")
+    mc = model()
+    myokit_model = mc.generate_myokit_model()
+    myokit.save(os.path.join(output_dir, f"{model.__name__}.mmt"), myokit_model)
+
+
+@pytest.mark.parametrize("model", models)
+def test_visualise_graphs(model):
+    name = model.__name__
+    logging.debug(f"initiating {name}")
+    mc = model()
+
+    mc.draw_graph(os.path.join(output_dir, f"{name}_graph.html"),
+                  show_parameters=False)
+    mc.draw_graph(os.path.join(output_dir, f"{name}_graph_with_parameters.html"),
+                  show_parameters=True)
+
+    nx.drawing.nx_agraph.write_dot(mc.graph, os.path.join(output_dir,
+                                                          "%s_dotfile.dot" % name))
+
+
+@pytest.mark.parametrize("model", models)
+def test_connected(model):
+    name = model.__name__
+    logging.debug(f"initiating {name}")
+
+    mc = model()
+    assert mc.is_connected() ^ (model in disconnected_models), f"model {model} is not connected"
+
+
+@pytest.mark.parametrize("model", models)
+def test_reversible(model):
+    name = model.__name__
+    logging.debug(f"initiating {name}")
+
+    if model in disconnected_models:
+        return
+
+    mc = model()
+    if not mc.is_reversible():
+        undirected_graph = mc.graph.to_undirected(reciprocal=False, as_view=True)
+        cycle_basis = nx.cycle_basis(undirected_graph)
+
+        for cycle in cycle_basis:
+            iterator = list(zip(cycle, itertools.islice(cycle, 1, None)))
+            forward_rate_list = [sp.sympify(mc.graph.get_edge_data(frm, to)['rate']) for frm, to in iterator]
+            backward_rate_list = [sp.sympify(mc.graph.get_edge_data(frm, to)['rate']) for to, frm in iterator]
+
+            logging.debug(mc.rate_expressions)
+
+            # Substitute in expressions
+            forward_rate_list = [rate.subs(mc.rate_expressions) for
+                                 rate in forward_rate_list]
+            backward_rate_list = [rate.subs(mc.rate_expressions) for rate in
+                                  backward_rate_list]
+
+            forward_rate_product = sp.prod(forward_rate_list)
+            backward_rate_product = sp.prod(backward_rate_list)
+            if (forward_rate_product - backward_rate_product).evalf() != 0:
+                logging.error("%s: rates moving forwards around the cycle are: %s", name, forward_rate_list)
+                logging.error("%s: rates moving backwards around the cycle are: %s", name, backward_rate_list)
+                logging.error(f"States are {cycle}")
+
+    assert mc.is_reversible(), "MC is not reversible"
+
+
+@pytest.mark.parametrize("model", models)
+def test_myokit_simulation_output(model):
+    mmt_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)),
+                           'models_myokit')
+    Erev = -88
+
+    # Extract number from model name i.e. model_02 -> 02
+    index = int(str(model.__name__)[-2:])
+
+    comparison_plot_dir = os.path.join(output_dir, 'comparison_plots')
+
+    if not os.path.exists(comparison_plot_dir):
+        os.makedirs(comparison_plot_dir)
+
+    name = model.__name__
+
+    mc = model()
+
+    mk_protocol_filename = os.path.join(mmt_dir,
+                                        'simplified-staircase.mmt')
+
+    mk_protocol = mk.load_protocol(mk_protocol_filename)
+
+    mk_model = mk.load_model(os.path.join(mmt_dir,
+                                          f"model-{index}.mmt"))
+
+    # logging.debug(f"Loaded model-{i}.mmt")
+
+    sim = mk.Simulation(mk_model, mk_protocol)
+    sim.set_tolerance(1e-9, 1e-9)
+    sim.set_constant('nernst.EK', Erev)
+
+    tmax = mk_protocol.characteristic_time()
+    times = np.linspace(0, tmax, int(tmax) + 1)
+    sim.pre(10000)
+
+    log = sim.run(tmax, log_times=times, log=['ikr.IKr'])
+
+    mk_IKr = np.array(log['ikr.IKr'])
+
+    remove = None
+    if model not in disconnected_models + [models[0], models[1], models[8]]:
+        # Find removed state
+        for state_var in mk_model["ikr"]:
+            str_state_var = str(state_var).split(".")[1]
+            print(state_var)
+            if (not state_var.is_state()) and (str_state_var in mc.graph):
+                remove = str_state_var
+                break
+        if remove is None:
+            assert remove, "Couldn't find redundant state to remove"
+
+    generated_mk_model = mc.generate_myokit_model(eliminate_state=remove)
+
+    sim = mk.Simulation(generated_mk_model, mk_protocol)
+    sim.set_tolerance(1e-9, 1e-9)
+    sim.set_constant('markov_chain.E_Kr', Erev)
+    sim.pre(10000)
+
+    log = sim.run(tmax, log_times=times, log=['markov_chain.I_Kr'])
+    gen_mk_IKr = np.array(log['markov_chain.I_Kr'])
+
+    fig = plt.figure(figsize=(12, 9))
+    fig.gca().plot(times[:-1], mk_IKr, label='original myokit simulation')
+    fig.gca().plot(times[:-1], gen_mk_IKr, label='generated markov_builder simulation')
+
+    fig.gca().legend()
+    fig.savefig(os.path.join(comparison_plot_dir,
+                             f"{name}_myokit_comparison"))
+    plt.close(fig)
+
+    error = np.sqrt(np.mean((gen_mk_IKr - mk_IKr)**2))
+
+    test_threshold = 0.02 if model in disconnected_models else 1e-5
+
+    assert error < test_threshold, f"Excessive error in model {index}: {name} {error} >= {test_threshold}"
+
+
+if __name__ == "__main__":
+    logging.getLogger().setLevel(logging.DEBUG)
+    unittest.main()
diff --git a/tests/test_MarkovChain.py b/tests/test_MarkovChain.py
index e4bef52..2b2edb4 100644
--- a/tests/test_MarkovChain.py
+++ b/tests/test_MarkovChain.py
@@ -4,284 +4,332 @@
 import os
 import unittest
 
+import matplotlib
 import matplotlib.pyplot as plt
 import myokit
 import networkx as nx
+import pytest
 import sympy as sp
 
 import markov_builder.example_models as example_models
-from markov_builder.rate_expressions import negative_rate_expr, positive_rate_expr
+from markov_builder.rate_expressions import (negative_rate_expr,
+                                             positive_rate_expr)
 
 
-class TestMarkovChain(unittest.TestCase):
+matplotlib.use('pdf')
 
-    def setUp(self):
-        """Run by unittest before the tests in this class are performed.
+models = [example_models.construct_four_state_chain(),
+          example_models.construct_non_reversible_chain(),
+          example_models.construct_mazhari_chain(),
+          example_models.construct_wang_chain(),
+          example_models.construct_kemp_model()]
 
-        Create an output directory (if it doesn't already exist). An
-        alternative output path can be used by setting the
-        MARKOVBUILDER_TEST_OUTPUT environment variable
+parameterised_models = [example_models.construct_four_state_chain(),
+                        example_models.construct_mazhari_chain(),
+                        example_models.construct_wang_chain(),
+                        example_models.construct_kemp_model()]
 
-        """
-        test_output_dir = os.environ.get('MARKOVBUILDER_TEST_OUTPUT', os.path.join(
-            'test_output', os.path.dirname(os.path.abspath(__file__)), self.__class__.__name__))
+current_file = os.path.basename(__file__)
+output_dir = os.environ.get('MARKOVBUILDER_TEST_OUTPUT', os.path.join(
+    'test_output', current_file))
 
-        if not os.path.exists(test_output_dir):
-            os.makedirs(test_output_dir)
-        self.output_dir = test_output_dir
-        logging.info("outputting to " + test_output_dir)
+if not os.path.exists(output_dir):
+    os.makedirs(output_dir)
 
-        self.models = [example_models.construct_four_state_chain(),
-                       example_models.construct_non_reversible_chain(),
-                       example_models.construct_mazhari_chain(),
-                       example_models.construct_wang_chain(),
-                       example_models.construct_kemp_model()]
+logging.info("outputting to " + output_dir)
 
-    def test_transition_matrix(self):
-        """Construct various examples of Markov models.
-        Output dot files of these graphs in the test output directory.
 
-        Check that our transition rate matrix for the Beattie model is correct
-        by comparing it against a reference solution.
+def test_transition_matrix():
+    """Construct various examples of Markov models.
+    Output dot files of these graphs in the test output directory.
 
-        """
+    Check that our transition rate matrix for the Beattie model is correct
+    by comparing it against a reference solution.
 
-        # Construct Beattie model
-        logging.info("Constructing four-state Beattie model")
+    """
 
-        mc = example_models.construct_four_state_chain()
+    # Construct Beattie model
+    logging.info("Constructing four-state Beattie model")
 
-        # Output DOT file
-        nx.drawing.nx_agraph.write_dot(mc.graph, "Beattie_dotfile.dot")
+    mc = example_models.construct_four_state_chain()
 
-        # Draw graph using pyvis
-        mc.draw_graph(os.path.join(self.output_dir, "BeattieModel.html"), show_options=True)
-        mc.draw_graph(os.path.join(self.output_dir, "BeattieModel.html"))
-        logging.debug(mc.graph)
+    # Output DOT file
+    nx.drawing.nx_agraph.write_dot(mc.graph, "Beattie_dotfile.dot")
+
+    # Draw graph using pyvis
+    mc.draw_graph(os.path.join(output_dir, "BeattieModel.html"), show_options=True)
+    mc.draw_graph(os.path.join(output_dir, "BeattieModel.html"))
+    logging.debug(mc.graph)
+
+    labels, Q = mc.get_transition_matrix(use_parameters=True)
+    labels, Q = mc.get_transition_matrix()
+
+    logging.debug("Q^T matrix is {}, labels are {}".format(Q.T, labels))
 
-        labels, Q = mc.get_transition_matrix(use_parameters=True)
-        labels, Q = mc.get_transition_matrix()
+    system = mc.eliminate_state_from_transition_matrix(['C', 'O', 'I'])
 
-        logging.debug("Q^T matrix is {}, labels are {}".format(Q.T, labels))
+    pen_and_paper_A = sp.Matrix([['-k_1 - k_3 - k_4', 'k_2 - k_4', '-k_4'],
+                                 ['k_1', '-k_2 - k_3', 'k_4'],
+                                 ['-k_1', 'k_3 - k_1', '-k_2 - k_4 - k_1']])
+    pen_and_paper_B = sp.Matrix(['k_4', 0, 'k_1'])
 
-        system = mc.eliminate_state_from_transition_matrix(['C', 'O', 'I'])
+    assert pen_and_paper_A == system[0]
+    assert pen_and_paper_B == system[1]
 
-        pen_and_paper_A = sp.Matrix([['-k_1 - k_3 - k_4', 'k_2 - k_4', '-k_4'],
-                                     ['k_1', '-k_2 - k_3', 'k_4'],
-                                     ['-k_1', 'k_3 - k_1', '-k_2 - k_4 - k_1']])
-        pen_and_paper_B = sp.Matrix(['k_4', 0, 'k_1'])
+    system2 = mc.eliminate_state_from_transition_matrix(['C', 'I', 'O'])
 
-        self.assertEqual(pen_and_paper_A, system[0])
-        self.assertEqual(pen_and_paper_B, system[1])
+    assert pen_and_paper_A[[0, 2, 1], [0, 2, 1]] == system2[0]
+    assert pen_and_paper_B[[0, 2, 1], 0] == system2[1]
 
-    def test_construct_examples(self):
-        """ Output the example models as pyvis files and DOT files """
-        for mc in self.models:
-            name = mc.name
-            mc.draw_graph(os.path.join(self.output_dir, "%s_graph.html" % name))
 
-            if len(mc.rate_expressions) > 0:
-                mc.draw_graph(os.path.join(self.output_dir, "%s_graph_parameters.html" %
-                                           name), show_parameters=True)
+@pytest.mark.parametrize("mc", models)
+def test_construct_examples(mc):
+    """ Output the example models as pyvis files and DOT files """
+    name = mc.name
+    mc.draw_graph(os.path.join(output_dir, "%s_graph.html" % name))
 
-            nx.drawing.nx_agraph.write_dot(mc.graph, "%s_dotfile.dot" % name)
-            nx.drawing.nx_agraph.write_dot(mc.graph, "%s_dotfile.dot" % name)
+    if len(mc.rate_expressions) > 0:
+        mc.draw_graph(os.path.join(output_dir, "%s_graph_parameters.html" %
+                                   name), show_parameters=True)
 
-    def test_parameterise_rates_no_default(self):
-        """Test parameterise rates using a dictionary with no default parameter values.
+    nx.drawing.nx_agraph.write_dot(mc.graph, "%s_dotfile.dot" % name)
 
-        Using the Beattie model
-        """
 
-        mc = example_models.construct_four_state_chain()
+def test_parameterise_rates_no_default():
+    """Test parameterise rates using a dictionary with no default parameter values.
 
-        rate_dictionary = {'k_1': positive_rate_expr,
-                           'k_2': negative_rate_expr,
-                           'k_3': positive_rate_expr,
-                           'k_4': negative_rate_expr,
-                           }
+    Using the Beattie model
+    """
 
-        mc.parameterise_rates(rate_dictionary, shared_variables=('V',))
+    mc = example_models.construct_four_state_chain()
 
-    def test_myokit_output(self):
-        """
-        Test the MarkovChain.parameterise_rates function.
-        """
+    rate_dictionary = {'k_1': positive_rate_expr,
+                       'k_2': negative_rate_expr,
+                       'k_3': positive_rate_expr,
+                       'k_4': negative_rate_expr,
+                       }
 
-        for mc in (example_models.construct_four_state_chain(), example_models.construct_kemp_model()):
-            # Expressions to be used for the rates. The variable V (membrane
-            # voltage) is shared across expressions and so it should only appear
-            # once in the parameter list.
+    mc.parameterise_rates(rate_dictionary, shared_variables={'V': 'V'})
 
-            # Output system of equations
-            logging.debug("ODE system is %s", str(mc.get_transition_matrix(use_parameters=True)))
 
-            # Output reduced system of equations
-            logging.debug("Reduced ODE system is :%s",
-                          str(mc.eliminate_state_from_transition_matrix(list(mc.graph.nodes)[:-1],
-                                                                        use_parameters=True)))
+@pytest.mark.parametrize("mc", parameterised_models)
+def test_myokit_output(mc):
+    """
+    Test the MarkovChain.parameterise_rates function.
+    """
 
-            # Output list of parameters
-            param_list = mc.get_parameter_list()
-            logging.debug("parameters are %s", mc.get_parameter_list())
+    # Expressions to be used for the rates. The variable V (membrane
+    # voltage) is shared across expressions and so it should only appear
+    # once in the parameter list.
 
-            self.assertEqual(param_list.count('V'), 1)
+    # Output system of equations
+    logging.debug("ODE system is %s", str(mc.get_transition_matrix(use_parameters=True)))
 
-            # Generate myokit code
-            myokit_model = mc.generate_myokit_model(eliminate_state='O')
-            myokit_model = mc.generate_myokit_model()
-            myokit_model = mc.generate_myokit_model(drug_binding=True)
-            myokit_model = mc.generate_myokit_model(drug_binding=True, eliminate_state='O')
+    # Output reduced system of equations
+    logging.debug("Reduced ODE system is :%s",
+                  str(mc.eliminate_state_from_transition_matrix(list(mc.graph.nodes)[:-1],
+                                                                use_parameters=True)))
 
-            myokit.save(os.path.join(self.output_dir, 'beattie_model.mmt'), myokit_model)
+    # Output list of parameters
+    param_list = mc.get_parameter_list()
+    logging.debug("parameters are %s", mc.get_parameter_list())
 
-            # Eliminate last node
-            myokit_model = mc.generate_myokit_model(list(mc.graph)[-1])
-            myokit.save(os.path.join(self.output_dir, 'beattie_model_reduced.mmt'), myokit_model)
+    assert param_list.count('V') == 1
 
-    def test_construct_open_trapping_model(self):
-        """
+    if 'O' in mc.graph.nodes():
+        eliminate_state = 'O'
+    else:
+        eliminate_state = None
 
-        Construct a model then run the add_open_trapping function to mirror the
-        model and add drug trapping to the open channel.
+    # Generate myokit code
+    myokit_model = mc.generate_myokit_model(eliminate_state=eliminate_state)
+    myokit_model = mc.generate_myokit_model()
 
-        """
+    if mc.is_connected():
+        myokit_model = mc.generate_myokit_model(drug_binding=True)
+        myokit_model = mc.generate_myokit_model(drug_binding=True,
+                                                eliminate_state=eliminate_state)
+        # Eliminate last node
+        myokit_model = mc.generate_myokit_model(list(mc.graph)[-1])
+        myokit.save(os.path.join(output_dir,
+                                 f"{mc.name}_model_reduced.mmt"), myokit_model)
 
-        for mc in self.models:
-            mc.add_open_trapping(prefix="d_", new_rates=True)
+    myokit.save(os.path.join(output_dir, f"{mc.name}_model.mmt"),
+                myokit_model)
 
-            # Save dotfile
-            nx.drawing.nx_agraph.write_dot(mc.graph, os.path.join(self.output_dir, "open_trapping.dot"))
 
-            mc.draw_graph(os.path.join(self.output_dir, "%s_open_trapping.html" % mc.name))
-            logging.debug(mc.graph)
+def test_construct_open_trapping_model():
+    """
 
-    def test_assert_reversibility_using_cycles(self):
-        """Test that MarkovChain().is_reversible correctly identifies if markov
-        chains are reversible or not.
+    Construct a model then run the add_open_trapping function to mirror the
+    model and add drug trapping to the open channel.
 
-        MarkovChain().is_reversible checks reversibility using the method
-        presented in Colquhoun et al. (2014)
-        https://doi.org/10.1529/biophysj.103.038679
+    """
 
-        TODO add more test cases
+    for mc in models:
 
-        """
+        if not mc.is_connected():
+            continue
 
-        # Test function on models that we know are reversible
-        reversible_models = [example_models.construct_four_state_chain(),
-                             example_models.construct_mazhari_chain(),
-                             example_models.construct_kemp_model()]
+        mc.add_open_trapping(prefix="d_", new_rates=True)
+
+        # Save dotfile
+        nx.drawing.nx_agraph.write_dot(mc.graph, os.path.join(output_dir, "open_trapping.dot"))
+
+        mc.draw_graph(os.path.join(output_dir, "%s_open_trapping.html" % mc.name))
+        logging.debug(mc.graph)
 
-        for mc in reversible_models:
-            logging.info("Checking reversibility")
-            self.assertTrue(mc.is_reversible())
-            logging.info("Checking reversibility with open trapping")
-            mc.add_open_trapping(new_rates=True)
-            self.assertTrue(mc.is_reversible())
 
-        # Test is_reversible on a model we know is not reversible. This example
-        # is a simple three state model with 6 independent transition rates
+def test_assert_reversibility_using_cycles():
+    """Test that MarkovChain().is_reversible correctly identifies if markov
+    chains are reversible or not.
 
-        logging.info("Checking reversibility of non-reversible chain")
-        mc = example_models.construct_non_reversible_chain()
-        logging.debug("graph is %s", mc.graph)
+    MarkovChain().is_reversible checks reversibility using the method
+    presented in Colquhoun et al. (2014)
+    https://doi.org/10.1529/biophysj.103.038679
 
-        self.assertFalse(mc.is_reversible())
-        logging.info("Checking reversibility of non-reversible chain")
-        mc.add_open_trapping()
-        logging.debug("graph is %s", mc.graph)
-        self.assertFalse(mc.is_reversible())
+    TODO add more test cases
 
-    def test_equate_rates(self):
-        """
-        Test that the MarkovChain.substitute_rates function performs the expected substitution
-
-        TODO: add more test cases
-        """
-        mc = example_models.construct_four_state_chain()
-        rates_dict = {'k_1': 'k_2*k_4'}
-        mc.substitute_rates(rates_dict)
-        mc.draw_graph(os.path.join(self.output_dir, '%s_rates_substitution.html' % mc.name), show_rates=True)
-        label_found = False
-        for _, _, d in mc.graph.edges(data=True):
-            if 'label' in d:
-                label_found = True
-                if d['label'] == 'k1':
-                    self.assertEqual(rates_dict['k1'], d['rate'])
-        self.assertTrue(label_found)
+    """
 
-        mc = example_models.construct_four_state_chain()
-        mc.substitute_rates(rates_dict)
+    # Test function on models that we know are reversible
+    reversible_models = [example_models.construct_four_state_chain(),
+                         example_models.construct_mazhari_chain(),
+                         example_models.construct_kemp_model()]
+
+    for mc in reversible_models:
+        # Skip models with multiple components
+        if not mc.is_connected():
+            continue
+
+        logging.info("Checking reversibility")
+        assert mc.is_reversible()
+        logging.info("Checking reversibility with open trapping")
         mc.add_open_trapping(new_rates=True)
-        mc.draw_graph(os.path.join(self.output_dir,
-                                   '%s_open_trapping_rates_substitution.html' %
-                                   mc.name), show_rates=True)
-        transition_rates = [d['rate'] for _, _, d in mc.graph.edges(data=True)]
+        assert mc.is_reversible()
+
+    # Test is_reversible on a model we know is not reversible. This example
+    # is a simple three state model with 6 independent transition rates
+
+    logging.info("Checking reversibility of non-reversible chain")
+    mc = example_models.construct_non_reversible_chain()
+    logging.debug("graph is %s", mc.graph)
+
+    assert not mc.is_reversible()
+    logging.info("Checking reversibility of non-reversible chain")
+    mc.add_open_trapping()
+    logging.debug("graph is %s", mc.graph)
+    assert not mc.is_reversible()
+
+
+def test_equate_rates():
+    """
+    Test that the MarkovChain.substitute_rates function performs the expected substitution
+
+    TODO: add more test cases
+    """
+    mc = example_models.construct_four_state_chain()
+    rates_dict = {'k_1': 'k_2*k_4'}
+    mc.substitute_rates(rates_dict)
+    mc.draw_graph(os.path.join(output_dir, '%s_rates_substitution.html' % mc.name), show_rates=True)
+    label_found = False
+    for _, _, d in mc.graph.edges(data=True):
+        if 'label' in d:
+            label_found = True
+            if d['label'] == 'k1':
+                assert rates_dict['k1'] == d['rate']
+    assert label_found
+
+    mc = example_models.construct_four_state_chain()
+    mc.substitute_rates(rates_dict)
+    mc.add_open_trapping(new_rates=True)
+    mc.draw_graph(os.path.join(output_dir,
+                               '%s_open_trapping_rates_substitution.html' %
+                               mc.name), show_rates=True)
+    transition_rates = [d['rate'] for _, _, d in mc.graph.edges(data=True)]
+
+    # Check that new mirrored rates have been handled correctly
+    assert 'd_k_2*d_k_4' in transition_rates
+    assert 'd_k_2' in mc.rates
+    assert 'd_k_4' in mc.rates
+
+    # Check reversibility still holds for good measure
+    assert mc.is_reversible()
+
+
+@pytest.mark.parametrize("mc", parameterised_models)
+def test_latex_printing(mc):
+    """ Test that we can generate LaTeX expressions for a model
+
+    TODO: Add more cases
+    """
+
+    logging.debug(f"Printing latex for {mc.name}")
+    logging.debug(mc.as_latex())
+    logging.debug(mc.as_latex(label_order=mc.get_states()))
+    logging.debug(mc.as_latex(include_auxiliary_expression=True))
+
+    if 'O' in mc.graph:
+        logging.debug(mc.as_latex('O', True))
+        logging.debug(mc.as_latex(state_to_remove='O'))
+
+
+@pytest.mark.parametrize("mc", parameterised_models)
+def test_sample_trajectories(mc):
+    """Simulate the 4-state Beattie Model using the Gillespie method
+
+    Also, compare these results against the equilibrium distribution
 
-        # Check that new mirrored rates have been handled correctly
-        self.assertIn('d_k_2*d_k_4', transition_rates)
-        self.assertIn('d_k_2', mc.rates)
-        self.assertIn('d_k_4', mc.rates)
+    TODO add more models
+    """
+    n_samples = 500
 
-        # Check reversibility still holds for good measure
-        self.assertTrue(mc.is_reversible())
+    mc = example_models.construct_four_state_chain()
 
-    def test_latex_printing(self):
-        """ Test that we can generate LaTeX expressions for a model
+    param_dict = mc.default_values
+    param_dict['V'] = 0
 
-        TODO: Add more cases
-        """
+    with pytest.raises(ValueError):
+        mc.get_equilibrium_distribution()
 
-        models = [example_models.construct_four_state_chain(), example_models.construct_wang_chain()]
-        for mc in models:
-            logging.debug(f"Printing latex for {mc.name}")
-            logging.debug(mc.as_latex())
-            logging.debug(mc.as_latex(label_order=mc.get_states()))
-            logging.debug(mc.as_latex(state_to_remove='O'))
-            logging.debug(mc.as_latex(include_auxiliary_expression=True))
-            logging.debug(mc.as_latex('O', True))
+    labels, eqm_dist = mc.get_equilibrium_distribution(param_dict=param_dict)
+    starting_distribution = [int(val) for val in n_samples * eqm_dist]
 
-    def test_sample_trajectories(self):
-        """Simulate the 4-state Beattie Model using the Gillespie method
+    fig = plt.figure()
+    ax = fig.subplots()
 
-        Also, compare these results against the equilibrium distribution
+    df = mc.sample_trajectories(n_samples, (0, 10))
 
-        TODO add more models
-        """
-        n_samples = 5
+    df = mc.sample_trajectories(n_samples, (0, 250), param_dict=param_dict,
+                                starting_distribution=starting_distribution)
 
-        mc = example_models.construct_four_state_chain()
+    for state in [lab for lab in df.columns if lab != 'time']:
+        df[state] = df[state] * 100.0 / n_samples
 
-        param_dict = mc.default_values
-        param_dict['V'] = 0
+    df = df.set_index('time')
+    logging.debug(f"sample trajectories results: {df}")
 
-        self.assertRaises(TypeError, mc.get_equilibrium_distribution)
+    ax.set_ylabel("State occupancy (%)")
+    ax.set_xlabel('time (ms)')
+    df.plot(ax=ax)
 
-        labels, eqm_dist = mc.get_equilibrium_distribution(param_dict=param_dict)
-        starting_distribution = [int(val) for val in n_samples * eqm_dist]
+    ax.spines[['top', 'right']].set_visible(False)
 
-        df = mc.sample_trajectories(n_samples, (0, 10))
+    fig.savefig(os.path.join(output_dir, 'beattie_model_sample_trajectories'),
+                transparent=True)
+    fig.clf()
 
-        df = mc.sample_trajectories(n_samples, (0, 250), param_dict=param_dict,
-                                    starting_distribution=starting_distribution)
-        df = df.set_index('time')
-        logging.debug(f"sample trajectories results: {df}")
-        df.plot()
 
-        for label, val in zip(labels, eqm_dist):
-            plt.axhline(val * n_samples, ls='--', alpha=0.5, label=label)
-        plt.legend()
+model_sizes = [(2, 3), (4, 5), (3, 4), (2, 8)]
 
-        plt.savefig(os.path.join(self.output_dir, 'beattie_model_sample_trajectories'))
-        logging.debug
 
-    def test_HH_models(self):
-        for i, j in [(2, 3), (4, 5), (3, 4), (2, 8)]:
-            mc = example_models.construct_HH_model(i, j)
-            mc.draw_graph(os.path.join(self.output_dir, mc.name + ".html"))
-            # update one rate
-            mc.add_both_transitions('O', 'C1', 'new_rate', 'new_rate2', update=True)
+@pytest.mark.parametrize("model_size", model_sizes)
+def test_HH_models(model_size):
+    i, j = model_size
+    mc = example_models.construct_HH_model(i, j)
+    mc.draw_graph(os.path.join(output_dir, mc.name + ".html"))
+    # update one rate
+    mc.add_both_transitions('O', 'C1', 'new_rate', 'new_rate2', update=True)
 
 
 if __name__ == "__main__":