Repository: myokit/myokit (language: python)

Issue:
SBML parser csymbol does not work
Parsing an SBML file with csymbols fails, e.g.

```
<?xml version="1.0" encoding="UTF-8"?>
<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" level="3" version="2">
  <model id="antibiotic_linear_growth_hill_kill_resistance_model" timeUnits="second">

    <listOfParameters>
      <parameter id="bacterial_count" value="1" constant="false"/>
    </listOfParameters>

    <listOfRules>
      <assignmentRule variable="bacterial_count">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <apply>
            <times/>
            <ci> bacterial_count </ci>
            <csymbol definitionURL="http://www.sbml.org/sbml/symbols/time"/>
          </apply>
        </math>
      </assignmentRule>
    </listOfRules>

  </model>
</sbml>
```

raises the following error

```
raise SBMLParsingError('Unable to parse rule: ' + str(e), element)
myokit.formats.sbml._parser.SBMLParsingError: Error on line 10. Unable to parse rule: Unable to parse MathML: Error on line 15. Unable to create Name from csymbol: Unknown or inaccessible symbol in MathML: "http://www.sbml.org/sbml/symbols/time".
```


Discussion / hints:
I checked the tests we wrote at the time and it seems that we only tested the SBML API for csymbols, but forgot to test the parser. 

It seems that the parser expects the csymbol to already be registered in the `variable_factory`, which is currently not the case. Would the best way to deal with this then to add supported csymbols to a global dictionary, such as the `assignable` dict in the parser @MichaelClerx ?
Good question!

You know more about this than me, but I suspect you might need to add a special variable to every sbml.Model to represent time, and then change the way that your `_add_time` method works when converting to a myokit.Model ?

(Related: looking back at the code now I was a little confused at the method _add_time, which doesn't change the sbml model but is used by the sbml-to-myokit conversion.Perhaps these methods should live in a separate class?)
Did we ever fix this @DavAug ?
@MichaelClerx Nope, but it's definitely on my to do list. I might have a look this weekend :)
No rush, just curious!

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c952fe8d6..eea3fa71f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -27,6 +27,7 @@ This page lists the main changes made to Myokit in each release.
   - [#786](https://github.com/MichaelClerx/myokit/pull/786) Added a method `myokit.Variable.clamp()` to fix any variable to a constant value.
   - [#786](https://github.com/MichaelClerx/myokit/pull/786) Added a method `myokit.Variable.remove_child_variables()` to remove all child (nested) variables from a given variable.
   - [#788](https://github.com/MichaelClerx/myokit/pull/788) The IDE will now suggest filenames for model exports.
+  - [#792](https://github.com/MichaelClerx/myokit/pull/792) SBML import now handles equations that explicitly use time.
 - Changed
   - [#581](https://github.com/MichaelClerx/myokit/pull/581) Powers are now rendered without spaces in mmt code, e.g. `x^2` instead of `x ^ 2`.
   - [#595](https://github.com/MichaelClerx/myokit/pull/595) The `Simulation` class now uses CVODES instead of CVODE as backend, which may require changes to your installation.
diff --git a/docs/source/api_formats/sbml_api.rst b/docs/source/api_formats/sbml_api.rst
index ef8ce4475..9a03629ca 100644
--- a/docs/source/api_formats/sbml_api.rst
+++ b/docs/source/api_formats/sbml_api.rst
@@ -20,6 +20,9 @@ SBML Model API
 .. autoclass:: Compartment
     :inherited-members:
 
+.. autoclass:: CSymbolVariable
+    :inherited-members:
+
 .. autoclass:: Model
 
 .. autoclass:: ModifierSpeciesReference
diff --git a/docs/source/api_index/myokit.formats.rst b/docs/source/api_index/myokit.formats.rst
index 2f1feb314..bdf3c207f 100644
--- a/docs/source/api_index/myokit.formats.rst
+++ b/docs/source/api_index/myokit.formats.rst
@@ -168,6 +168,7 @@ myokit.formats.python
 myokit.formats.sbml
 -------------------
 - :class:`myokit.formats.sbml.Compartment`
+- :class:`myokit.formats.sbml.CSymbolVariable`
 - :class:`myokit.formats.sbml.Model`
 - :class:`myokit.formats.sbml.ModifierSpeciesReference`
 - :meth:`myokit.formats.sbml.importers`
diff --git a/myokit/formats/sbml/__init__.py b/myokit/formats/sbml/__init__.py
index 682430216..1507fb65a 100644
--- a/myokit/formats/sbml/__init__.py
+++ b/myokit/formats/sbml/__init__.py
@@ -9,6 +9,7 @@
 
 from ._api import (  # noqa
     Compartment,
+    CSymbolVariable,
     Model,
     ModifierSpeciesReference,
     Parameter,
diff --git a/myokit/formats/sbml/_api.py b/myokit/formats/sbml/_api.py
index 2dbf61131..4097b813f 100644
--- a/myokit/formats/sbml/_api.py
+++ b/myokit/formats/sbml/_api.py
@@ -9,6 +9,8 @@
 # Regex for id checking
 _re_id = re.compile(r'^[a-zA-Z_]+[a-zA-Z0-9_]*$')
 
+_SBML_TIME = 'http://www.sbml.org/sbml/symbols/time'
+
 
 class SBMLError(Exception):
     """Raised if something goes wrong when working with an SBML model."""
@@ -62,7 +64,10 @@ def set_value(self, value, is_rate=False):
         Arguments:
 
         ``value``
-            An expression
+            An expression. The values of any :class:`myokit.Name` objects in
+            the expression must be either subclasses of
+            :class:`myokit.formats.sbml.api.Quantity` or of
+            :class:`myokit.formats.sbml.api.CSymbolVariable`.
         ``rate``
             Set to ``True`` if the expression gives the rate of change of this
             variable.
@@ -84,6 +89,22 @@ def value(self):
         return self._value
 
 
+class CSymbolVariable(object):
+    """
+    Represents a CSymbol that can appear in SBML expressions, but which has a
+    predetermind value and/or meaning, e.g. "time".
+    """
+    def __init__(self, definition_url):
+        self._definition_url = str(definition_url)
+
+    def definition_url(self):
+        """ Returns the ``definitionUrl`` for this ``CSymbolVariable``. """
+        return self._definition_url
+
+    def __str__(self):
+        return '<CSymbolVariable ' + self._definition_url + '>'
+
+
 class Compartment(Quantity):
     """
     Represents a compartment in SBML; to create a compartment use
@@ -191,6 +212,7 @@ def __init__(self, name=None):
         self._sids = set()
 
         # Assignables: Compartments, species, species references, parameters
+        # Mapping from sID to Quantity subclass instance
         self._assignables = {}
 
         # Compartments, parameters, species (all maps from sids to objects)
@@ -216,6 +238,9 @@ def __init__(self, name=None):
         # Global conversion factor for species
         self._conversion_factor = None
 
+        # CSymbolVariable for time
+        self._time = CSymbolVariable(_SBML_TIME)
+
     def add_compartment(self, sid):
         """Adds a :class:`myokit.formats.sbml.Compartment` to this model."""
 
@@ -297,6 +322,17 @@ def assignable(self, sid):
         """
         return self._assignables[sid]
 
+    def assignable_or_csymbol(self, identifier):
+        """
+        Like :meth:`assignable`, but will also return a
+        :class:`CSymbolVariable` if the given identifier is a known CSymbol
+        ``definitionUrl``.
+        """
+        # Continas "://"? Then can't be a sid. So this is OK.
+        if identifier == _SBML_TIME:
+            return self._time
+        return self._assignables[identifier]
+
     def base_unit(self, unitsid):
         """
         Returns an SBML base unit, raises an :class:`SBMLError` if not
@@ -348,79 +384,7 @@ def myokit_model(self):
 
         Compartments defined by the SBML file are mapped to Myokit Components.
         """
-
-        # Create myokit model
-        myokit_model = myokit.Model(self.name())
-
-        # Add notes
-        notes = self.notes()
-        if notes:
-            myokit_model.meta['desc'] = notes
-
-        # Create reference container that links sid's to myokit objects
-        component_references = {}
-        variable_references = {}
-
-        # Create additional container for species variables in amount
-        # (Reactions define rates exclusively for amount of species)
-        species_amount_references = {}
-
-        # Create reference from Model expressions to myokit model expressions
-        expression_references = {}
-
-        # Instantiate component and variable objects first without assigning
-        # RHS expressions. Myokit objects may have to be renamed, so
-        # expressions are added in a second step.
-
-        # Add SBML compartments to myokit model
-        self._add_compartments(
-            myokit_model, component_references, variable_references,
-            expression_references)
-
-        # Add myokit component to model
-        # It stores time bound variable and global parameters
-        myokit_component_name = self._add_myokit_component(
-            myokit_model, component_references)
-
-        # Add species to components
-        self._add_species(
-            component_references, species_amount_references,
-            variable_references, expression_references)
-
-        # Add parameters to myokit component
-        self._add_parameters(
-            component_references[myokit_component_name], variable_references,
-            expression_references)
-
-        # Add stoichiometries from reactions
-        self._add_stoichiometries(
-            component_references, variable_references, expression_references)
-
-        # Add time variable to myokit component
-        self._add_time(
-            component_references[myokit_component_name], variable_references)
-
-        # Set RHS of compartment sizes
-        self._set_rhs_sizes(variable_references, expression_references)
-
-        # Set RHS of species (initalAssignemnt, assignmentRule, rateRule)
-        self._set_rhs_species(
-            species_amount_references, variable_references,
-            expression_references)
-
-        # Set RHS of parameters
-        self._set_rhs_parameters(variable_references, expression_references)
-
-        # Set RHS of stoichiometries
-        self._set_rhs_stoichiometries(
-            variable_references, expression_references)
-
-        # Set RHS of species changed by reactions
-        self._set_reactions(
-            variable_references, species_amount_references,
-            expression_references)
-
-        return myokit_model
+        return _MyokitConverter.convert(self)
 
     def name(self):
         """Returns this model's name."""
@@ -534,6 +498,12 @@ def substance_units(self):
         """
         return self._substance_units
 
+    def time(self):
+        """
+        Returns the :class:`CSymbolVariable` representing time in this model.
+        """
+        return self._time
+
     def time_units(self):
         """Returns the default units for time, or dimensionless if not set."""
         return self._time_units
@@ -557,1062 +527,1166 @@ def volume_units(self):
         """
         return self._volume_units
 
-    def _add_compartments(
-            self, myokit_model, component_references, variable_references,
-            expression_references):
-        """
-        Creates components for each compartment.
-        """
-        for sid, compartment in self._compartments.items():
-            # Create component for compartment
-            component = myokit_model.add_component_allow_renaming(
-                convert_name(sid))
+    # SBML base units (except Celsius, because it's not defined in myokit)
+    _base_units = {
+        'ampere': myokit.units.A,
+        'avogadro': myokit.parse_unit('1 (6.02214179e23)'),
+        'becquerel': myokit.units.Bq,
+        'candela': myokit.units.cd,
+        'coulomb': myokit.units.C,
+        'dimensionless': myokit.units.dimensionless,
+        'farad': myokit.units.F,
+        'gram': myokit.units.g,
+        'gray': myokit.units.Gy,
+        'henry': myokit.units.H,
+        'hertz': myokit.units.Hz,
+        'item': myokit.units.dimensionless,  # Myokit does not have item unit
+        'joule': myokit.units.J,
+        'katal': myokit.units.kat,
+        'kelvin': myokit.units.K,
+        'kilogram': myokit.units.kg,
+        'liter': myokit.units.L,
+        'litre': myokit.units.L,
+        'lumen': myokit.units.lm,
+        'lux': myokit.units.lux,
+        'meter': myokit.units.m,
+        'metre': myokit.units.m,
+        'mole': myokit.units.mol,
+        'newton': myokit.units.N,
+        'ohm': myokit.units.ohm,
+        'pascal': myokit.units.Pa,
+        'radian': myokit.units.rad,
+        'second': myokit.units.s,
+        'siemens': myokit.units.S,
+        'sievert': myokit.units.Sv,
+        'steradian': myokit.units.sr,
+        'tesla': myokit.units.T,
+        'volt': myokit.units.V,
+        'watt': myokit.units.W,
+        'weber': myokit.units.Wb,
+    }
 
-            # Add component to reference list
-            component_references[sid] = component
 
-            # Add compartment size to component
-            var = component.add_variable_allow_renaming('size')
+class Parameter(Quantity):
+    """
+    Represents a parameter in SBML; to create a parameter use
+    :meth:`model.add_parameter()`.
 
-            # Set unit of size variable
-            var.set_unit(compartment.size_units())
+    Arguments:
 
-            # Add size variable to reference list
-            variable_references[sid] = var
-            expression_references[myokit.Name(compartment)] = myokit.Name(var)
+    ``model``
+        The model this parameter is in.
+    ``sid``
+        This parameter's SId.
 
-    def _add_myokit_component(self, myokit_model, component_references):
-        """
-        Creates myokit component. It is used to store global parameters of the
-        model and the time bound variable.
-        """
+    """
+    def __init__(self, model, sid):
+        super(Parameter, self).__init__()
 
-        # Create default name of component
-        myokit_component_name = 'myokit'
+        if not isinstance(model, Model):
+            raise SBMLError(
+                '<' + str(model) + '> needs to be an instance of '
+                'myokit.formats.sbml.Model.')
 
-        # Add component
-        component = myokit_model.add_component_allow_renaming(
-            myokit_component_name)
+        self._model = model
+        self._sid = str(sid)
+        self._units = None
 
-        # Add myokit component to referece list
-        component_references[myokit_component_name] = component
+    def set_units(self, units):
+        """Sets this parameters units to the given ``units``."""
+        if not isinstance(units, myokit.Unit):
+            raise SBMLError(
+                '<' + str(units) + '> needs to be instance of myokit.Unit')
 
-        return myokit_component_name
+        self._units = units
 
-    def _add_parameters(
-            self, component, variable_references, expression_references):
-        """
-        Adds global parameters to the myokit component.
-        """
+    def sid(self):
+        """Returns this parameter's sid."""
+        return self._sid
 
-        for sid, parameter in self._parameters.items():
-            # Add parameter to component
-            var = component.add_variable_allow_renaming(convert_name(sid))
+    def __str__(self):
+        return '<Parameter ' + self._sid + '>'
 
-            # Set unit of parameter
-            var.set_unit(parameter.units())
+    def units(self):
+        """Returns the units this parameter is in, or ``None`` if not set."""
+        return self._units
 
-            # Add reference to parameter
-            variable_references[sid] = var
-            expression_references[myokit.Name(parameter)] = myokit.Name(var)
 
-    def _add_species(
-            self, component_references, species_amount_references,
-            variable_references, expression_references):
-        """
-        Adds amount (and potentially concentration) variables for each species
-        to the respective components.
+class Reaction(object):
+    """
+    Represents an SBML reaction; to create a reaction use
+    :meth:`Model.add_reaction()`.
 
-        RHS of species concentration is defined as species amount / size
-        compartment.
-        """
+    Arguments:
 
-        for sid, species in self._species.items():
-            # Get component from reference list
-            compartment = species.compartment()
-            compartment_sid = compartment.sid()
-            try:
-                component = component_references[compartment_sid]
-            except KeyError:
-                raise SBMLError(
-                    'The <' + str(compartment) + '> for <' + str(species) + '>'
-                    ' is not referenced in model. Please use the '
-                    '`add_compartment` method to reference a compartment.')
+    ``model``
+        The model this reaction is in.
+    ``sid``
+        This reaction's SId.
 
-            # Add species in amount to component
-            # (needed for reactions, even if species is defined in
-            # concentration)
-            var = component.add_variable_allow_renaming(
-                convert_name(sid + '_amount'))
+    """
+    def __init__(self, model, sid):
 
-            # Set unit of amount
-            var.set_unit(species.substance_units())
+        if not isinstance(model, Model):
+            raise SBMLError(
+                '<' + str(model) + '> needs to be an instance of '
+                'myokit.formats.sbml.Model.')
 
-            # Add reference to amount variable
-            species_amount_references[sid] = var
+        self._model = model
+        self._sid = str(sid)
 
-            # Add species in concentration if measured in concentration
-            if not species.is_amount():
-                # Add species in concentration
-                var = component.add_variable_allow_renaming(
-                    convert_name(sid + '_concentration'))
+        # Reactants, reaction products, and modifiers: all as SpeciesReference
+        self._reactants = []
+        self._products = []
+        self._modifiers = []
 
-                # Get myokit amount and size variable
-                amount = species_amount_references[sid]
-                size = variable_references[compartment_sid]
+        # All species involved in this reaction (sid to object)
+        self._species = collections.OrderedDict()
 
-                # Set unit of concentration
-                var.set_unit(amount.unit() / size.unit())
+        # The kinetic law specifying this reaction's rate (if set)
+        self._kinetic_law = None
 
-                # Define RHS of concentration as amount / size
-                rhs = myokit.Divide(
-                    myokit.Name(amount),
-                    myokit.Name(size))
-                var.set_rhs(rhs)
+    def add_modifier(self, species, sid=None):
+        """Adds a modifier to this reaction and returns the created object."""
+        if not isinstance(species, Species):
+            raise SBMLError(
+                '<' + str(species) + '> needs to be an instance of '
+                'myokit.formats.sbml.Species')
 
-            # Add reference to species (either in amount or concentration)
-            variable_references[sid] = var
-            expression_references[myokit.Name(species)] = myokit.Name(var)
+        sid = str(sid) if sid else sid
+        if sid is not None:
+            self._model._register_sid(sid)
+        ref = ModifierSpeciesReference(species, sid)
+        self._modifiers.append(ref)
+        self._species[species.sid()] = species
+        return ref
 
-    def _add_stoichiometries(
-            self, component_references, variable_references,
-            expression_references):
+    def add_product(self, species, sid=None):
         """
-        Adds stoichiometry parameters to the 'compartment' component to which
-        the species belongs.
+        Adds a reaction product to this reaction and returns the created
+        object.
         """
+        if not isinstance(species, Species):
+            raise SBMLError(
+                '<' + str(species) + '> needs to be an instance of '
+                'myokit.formats.sbml.Species')
 
-        for reaction in self._reactions.values():
-            # Get all reactants and products (modifier do not have
-            # stoichiometries)
-            species_references = reaction.reactants() + reaction.products()
+        sid = str(sid) if sid else sid
+        if sid is not None:
+            self._model._register_sid(sid)
+        ref = SpeciesReference(species, sid)
+        self._products.append(ref)
+        self._species[species.sid()] = species
+        if sid is not None:
+            self._model._assignables[sid] = ref
+        return ref
 
-            for species_reference in species_references:
-                sid = species_reference.sid()
-                if sid is not None:
-                    # Get component
-                    species = species_reference.species()
-                    compartment = species.compartment()
-                    compartment_sid = compartment.sid()
-                    try:
-                        component = component_references[compartment_sid]
-                    except KeyError:
-                        raise SBMLError(
-                            'The <' + str(compartment) + '> of <' +
-                            str(species) + '> in <' + str(reaction) + '> is '
-                            'not referenced in the model. Please use the '
-                            '`add_compartment` method to reference the '
-                            'compartment.')
+    def add_reactant(self, species, sid=None):
+        """Adds a reactant to this reaction and returns the created object."""
+        if not isinstance(species, Species):
+            raise SBMLError(
+                '<' + str(species) + '> needs to be an instance of '
+                'myokit.formats.sbml.Species')
 
-                    # Add variable to component
-                    var = component.add_variable_allow_renaming(
-                        convert_name(sid))
+        sid = str(sid) if sid else sid
+        if sid is not None:
+            self._model._register_sid(sid)
+        ref = SpeciesReference(species, sid)
+        self._reactants.append(ref)
+        self._species[species.sid()] = species
+        if sid is not None:
+            self._model._assignables[sid] = ref
+        return ref
 
-                    # Set unit of variable
-                    # (SBML defines stoichiometries as dimensionless)
-                    var.set_unit(myokit.units.dimensionless)
+    def kinetic_law(self):
+        """
+        Returns the kinetic law set for this reaction, or ``None`` if not set.
+        """
+        return self._kinetic_law
 
-                    # Add reference to variable
-                    variable_references[sid] = var
-                    expression_references[
-                        myokit.Name(species_reference)] = myokit.Name(var)
+    def modifiers(self):
+        """
+        Returns the list of modifiers used by this reaction (as
+        :class:`SpeciesReference` objects).
+        """
+        return list(self._modifiers)
 
-    def _add_time(self, component, variable_references):
+    def products(self):
         """
-        Adds time bound variable to the myokit compartment.
+        Returns the list of products created by this reaction (as
+        :class:`SpeciesReference` objects).
         """
+        return list(self._products)
 
-        # Add variable
-        var = component.add_variable_allow_renaming('time')
+    def reactants(self):
+        """
+        Returns the list of reactions used by this reaction (as
+        :class:`SpeciesReference` objects).
+        """
+        return list(self._reactants)
 
-        # Bind time variable to time in myokit model
-        var.set_binding('time')
+    def set_kinetic_law(self, expression):
+        """
+        Sets this reaction's kinetic law (as a :class:`myokit.Expression`).
+        """
+        if not isinstance(expression, myokit.Expression):
+            raise SBMLError(
+                '<' + str(expression) + '> needs to be an instance of '
+                'myokit.Expression.')
 
-        # Set time unit and initial value (SBML default: t=0)
-        var.set_unit(self.time_units())
-        var.set_rhs(0)
+        self._kinetic_law = expression
 
-        # Add reference to time variable
-        # (SBML referes to time by a csymbol:
-        # 'http://www.sbml.org/sbml/symbols/time/')
-        variable_references[
-            'http://www.sbml.org/sbml/symbols/time'] = var
+    def sid(self):
+        """Returns this reaction's sid."""
+        return self._sid
 
-    def _set_rhs_sizes(self, variable_references, expression_references):
+    def species(self, sid):
         """
-        Sets right hand side of compartments' size variables.
+        Finds and returns a :class:`Species` used in this reaction.
         """
+        return self._species[sid]
 
-        for sid, compartment in self._compartments.items():
-            # Get myokit variable
-            var = variable_references[sid]
+    def __str__(self):
+        return '<Reaction ' + self._sid + '>'
 
-            # Set initial value
-            expr = compartment.initial_value()
-            if expr is not None:
-                try:
-                    var.set_rhs(expr.clone(
-                        subst=expression_references))
-                except AttributeError:
-                    raise SBMLError(
-                        'Initial value for the size of <' + str(compartment) +
-                        '> contains unreferenced parameters/variables. Please'
-                        ' use e.g. the `add_parameter` method to add reference'
-                        ' to parameters in the model.')
 
-            if compartment.is_rate():
-                # Get initial state
-                try:
-                    state_value = var.eval()
-                except AttributeError:
-                    state_value = 1
-                    warnings.warn(
-                        'Size of compartment <' + str(compartment) + '> is '
-                        'promoted to state variable without being assigned '
-                        'with an initial value. Default is set to 1.')
+class Species(Quantity):
+    """
+    Represents an SBML species; to create a species use
+    :meth:`Compartment.add_species()`.
 
-                # Promote size to state variable
-                var.promote(state_value=state_value)
+    Arguments:
 
-            # Set RHS
-            # (assignmentRule overwrites initialAssignment)
-            expr = compartment.value()
-            if expr is not None:
-                try:
-                    var.set_rhs(expr.clone(
-                        subst=expression_references))
-                except AttributeError:
-                    raise SBMLError(
-                        'Value for the size of <' + str(compartment) +
-                        '> contains unreferenced parameters/variables. Please'
-                        ' use e.g. the `add_parameter` method to add reference'
-                        ' to parameters in the model.')
+    ``compartment``
+        The :class:`Compartment` that this species is in.
+    ``sid``
+        This species's SId.
+    ``is_amount``
+        Whether this species value is represented as an amount (if false, it is
+        represented as a concentration, which depends on the size of the
+        compartment it is in).
+    ``is_constant``
+        Whether or not this species is constant.
+    ``is_boundary``
+        Whether or not this species is at the boundary of a reaction.
 
-    def _set_reactions(
-            self, variable_references, species_amount_references,
-            expression_references):
+    """
+    def __init__(self, compartment, sid, is_amount, is_constant, is_boundary):
+        super(Species, self).__init__()
 
-        for reaction in self._reactions.values():
-            if reaction.kinetic_law() is None:
-                # Skip to the next reaction
-                continue
+        if not isinstance(compartment, Compartment):
+            raise SBMLError(
+                '<' + compartment + '> needs to be instance of'
+                'myokit.formats.sbml.Compartment')
+        if not isinstance(is_amount, bool):
+            raise SBMLError(
+                'Is_amount <' + str(is_amount) + '> needs to be a boolean.')
+        if not isinstance(is_constant, bool):
+            raise SBMLError(
+                'Is_constant <' + str(is_constant) + '> needs to be a boolean.'
+            )
+        if not isinstance(is_boundary, bool):
+            raise SBMLError(
+                'Is_boundary <' + str(is_boundary) + '> needs to be a boolean.'
+            )
+        self._compartment = compartment
 
-            # Set right hand side of reactants
-            self._set_rhs_reactants(
-                reaction, variable_references, species_amount_references,
-                expression_references)
+        self._sid = str(sid)
+        self._is_amount = bool(is_amount)
+        self._is_constant = bool(is_constant)
+        self._is_boundary = bool(is_boundary)
 
-            # Set right hand side of products
-            self._set_rhs_products(
-                reaction, variable_references, species_amount_references,
-                expression_references)
+        # Units for the amount, not the concentration
+        self._units = None
 
-    def _set_rhs_reactants(
-            self, reaction, variable_references, species_amount_references,
-            expression_references):
+        # Optional conversion factor from substance to extent units
+        self._conversion_factor = None
+
+        # Flag whether initial value is provided in amount or concentration
+        self._initial_value_in_amount = None
+
+    def is_amount(self):
         """
-        Sets right hand side of species acting as reactants in a reaction.
+        Returns ``True`` only if this species is defined as an amount (not a
+        concentration).
         """
+        return self._is_amount
 
-        for reactant in reaction.reactants():
-            # Get species object
-            species = reactant.species()
+    def is_boundary(self):
+        """Returns ``True`` only if this species is at a reaction boundary."""
+        return self._is_boundary
 
-            if species.is_constant() or species.is_boundary():
-                # Species is not altered by the reaction
-                # Skip to the next reactant
-                continue
+    def compartment(self):
+        """Returns the :class:`Compartment` that this species is in."""
+        return self._compartment
 
-            # Instantiate rate expression
-            expr = reaction.kinetic_law().clone()
+    def initial_value(self):
+        """
+        Returns a :class:`myokit.Expression` (or None) for this species'
+        initial value, and a boolean indicating whether initial value is
+        provided in amount (True) or in concentration (False).
+        """
+        return self._initial_value, self._initial_value_in_amount
 
-            # Get stoichiometry of reactant
-            try:
-                stoichiometry = myokit.Name(
-                    variable_references[reactant.sid()])
-            except KeyError:
-                stoichiometry = reactant.initial_value()
+    def is_constant(self):
+        """Returns ``True`` only if this species is constant."""
+        return self._is_constant
 
-            if stoichiometry is not None:
-                # Weight rate expression by stoichiometry
-                expr = myokit.Multiply(stoichiometry, expr)
+    def conversion_factor(self):
+        """
+        Returns the :class:`Parameter` acting as conversion factor for this
+        species, or the model default if not set (and ``None`` if that isn't
+        set either).
 
-            factor = species.conversion_factor()
-            if factor is not None:
-                # Get conversion factor variable
-                sid = factor.sid()
-                try:
-                    conversion_factor = variable_references[sid]
-                except KeyError:
-                    raise SBMLError(
-                        'Species <' + str(species) + '> has conversion factor '
-                        '<' + str(factor) + '> which is not referenced in the '
-                        'model.')
+        When calculating the rate of change of the species, multiplying by this
+        factor converts from the species units to "number of reaction events".
+        """
+        if self._conversion_factor is not None:
+            return self._conversion_factor
+        return self._compartment._model.conversion_factor()
 
-                # Convert rate expression from units of reaction extent
-                # to amount units
-                conversion_factor = myokit.Name(conversion_factor)
-                expr = myokit.Multiply(conversion_factor, expr)
+    def set_conversion_factor(self, factor):
+        """
+        Sets a :class:`Parameter` as conversion factor for this species,
+        see :meth:`conversion_factor()`.
+        """
+        if not isinstance(factor, Parameter):
+            raise SBMLError(
+                '<' + str(factor) + '> needs to be instance of'
+                'myokit.formats.sbml.Parameter.')
 
-            # Get myokit amount variable
-            try:
-                var = species_amount_references[species.sid()]
-            except KeyError:
-                raise SBMLError(
-                    'Kinetic law of <' + str(reaction) + '> contains '
-                    'unreferenced species <' + str(species) + '>. Please, '
-                    'use the `add_species` method to reference species.')
+        self._conversion_factor = factor
 
-            if not var.is_state():
-                # Get initial state
-                try:
-                    state_value = var.eval()
-                except AttributeError:
-                    state_value = 0
-                    warnings.warn(
-                        'Species <' + str(species) + '> is promoted to state '
-                        'variable without being assigned with an initial '
-                        'value. Default is set to 0.')
+    def set_initial_value(self, value, in_amount=None):
+        """
+        Sets a :class:`myokit.Expression` for this species' initial value.
 
-                # Promote size to state variable
-                var.promote(state_value=state_value)
-                var.set_rhs(myokit.Number(0))
+        `in_amount` is a boolean that indicates whether the initial value is
+        measured in amount (True) or concentration (False). If set to `None`
+        value is treated to have the same units as the species.
+        """
+        if not isinstance(value, myokit.Expression):
+            raise SBMLError(
+                '<' + str(value) + '> needs to be an instance of '
+                'myokit.Expression.')
+        if (in_amount is not None) and (not isinstance(in_amount, bool)):
+            raise SBMLError(
+                '<in_amount> needs to be an instance of bool or None.')
 
-            if var.eval():
-                # Subtract rate contributions
-                # (Reaction removes species from compartment)
-                expr = myokit.Minus(var.rhs(), expr)
-            else:
-                expr = myokit.PrefixMinus(expr)
+        self._initial_value = value
+        self._initial_value_in_amount = in_amount
 
-            # Set RHS
-            try:
-                var.set_rhs(expr.clone(subst=expression_references))
-            except AttributeError:
-                raise SBMLError(
-                    'Reaction rate expression of <' + str(reaction) + '> for <'
-                    + str(species) + '> contains unreferenced parameters/'
-                    'variables. Please use e.g. the `add_parameter` method to '
-                    'add reference to parameters in the model.')
+    def set_substance_units(self, units):
+        """Sets the units this species amount (not concentration) is in."""
+        if not isinstance(units, myokit.Unit):
+            raise SBMLError(
+                '<' + str(units) + '> needs to be instance of myokit.Unit')
 
-    def _set_rhs_products(
-            self, reaction, variable_references, species_amount_references,
-            expression_references):
+        self._units = units
+
+    def sid(self):
+        """Returns this species's sid."""
+        return self._sid
+
+    def __str__(self):
+        return '<Species ' + self._sid + '>'
+
+    def substance_units(self):
         """
-        Sets right hand side of species acting as products in a reaction.
+        Returns the units an amount of this species (not a concentration) is
+        in, or the model default if not set.
         """
+        if self._units is not None:
+            return self._units
+        return self._compartment._model.substance_units()
 
-        for product in reaction.products():
-            # Get species object
-            species = product.species()
 
-            if species.is_constant() or species.is_boundary():
-                # Species is not altered by the reaction
-                # Skip to the next product
-                continue
+class SpeciesReference(Quantity):
+    """
+    Represents a reference to a reactant or product species in an SBML
+    reaction.
 
-            # Instantiate rate expression
-            expr = reaction.kinetic_law().clone()
+    A species reference acts as a :class:`Quantity`, where the value represents
+    the species' stoichiometry.
+    """
+    REACTANT = 0
+    PRODUCT = 1
+    MODIFIER = 2
 
-            # Get stoichiometry of product
-            try:
-                stoichiometry = myokit.Name(
-                    variable_references[product.sid()])
-            except KeyError:
-                stoichiometry = product.initial_value()
+    def __init__(self, species, sid=None):
+        super(SpeciesReference, self).__init__()
 
-            if stoichiometry is not None:
-                # Weight rate expression by stoichiometry
-                expr = myokit.Multiply(stoichiometry, expr)
+        if not isinstance(species, Species):
+            raise SBMLError(
+                '<' + species + '> needs to be instance of'
+                'myokit.formats.sbml.Species')
 
-            factor = species.conversion_factor()
-            if factor is not None:
-                # Get conversion factor variable
-                sid = species.conversion_factor().sid()
-                try:
-                    conversion_factor = variable_references[sid]
-                except KeyError:
-                    raise SBMLError(
-                        'Species <' + str(species) + '> has conversion factor '
-                        '<' + str(factor) + '> which is not referenced in the '
-                        'model.')
+        self._species = species
+        self._sid = str(sid) if sid else sid
 
-                # Convert rate expression from units of reaction extent
-                # to amount units
-                conversion_factor = myokit.Name(conversion_factor)
-                expr = myokit.Multiply(conversion_factor, expr)
+    def species(self):
+        """Returns the species this object refers to."""
+        return self._species
 
-            # Get myokit amount variable
-            try:
-                var = species_amount_references[species.sid()]
-            except KeyError:
-                raise SBMLError(
-                    'Kinetic law of <' + str(reaction) + '> contains '
-                    'unreferenced species <' + str(species) + '>. Please, '
-                    'use the `add_species` method to reference species.')
+    def sid(self):
+        """Returns this species reference's SId, or ``None`` if not set."""
+        return self._sid
 
-            if not var.is_state():
-                # Get initial state
-                try:
-                    state_value = var.eval()
-                except AttributeError:
-                    state_value = 0
-                    warnings.warn(
-                        'Species <' + str(species) + '> is promoted to state '
-                        'variable without being assigned with an initial '
-                        'value. Default is set to 0.')
+    def __str__(self):
+        return '<SpeciesReference ' + self._sid + '>'
 
-                # Promote size to state variable
-                var.promote(state_value=state_value)
-                var.set_rhs(myokit.Number(0))
 
-            if var.rhs().eval():
-                # Add rate contributions
-                expr = myokit.Plus(var.rhs(), expr)
+class ModifierSpeciesReference(object):
+    """Represents a reference to a modifier species in an SBML reaction."""
 
-            # Set RHS
-            try:
-                var.set_rhs(expr.clone(subst=expression_references))
-            except AttributeError:
-                raise SBMLError(
-                    'Reaction rate expression of <' + str(reaction) + '> for <'
-                    + str(species) + '> contains unreferenced parameters/'
-                    'variables. Please use e.g. the `add_parameter` method to '
-                    'add reference to parameters in the model.')
+    def __init__(self, species, sid=None):
+        super(ModifierSpeciesReference, self).__init__()
 
-    def _set_rhs_species(
-            self, species_amount_references, variable_references,
-            expression_references):
-        """
-        Sets right hand side of species amount variables defined by
-        assignments.
+        if not isinstance(species, Species):
+            raise SBMLError(
+                '<' + species + '> needs to be instance of'
+                'myokit.formats.sbml.Species')
 
-        Rate expressions defined by reactions are dealt with in
-        :meth:`_set_reactions`.
-        """
+        self._species = species
+        self._sid = str(sid) if sid else sid
 
-        for sid, species in self._species.items():
-            # Get myokit variable
-            # We only adapt amount of species
-            var = species_amount_references[sid]
+    def species(self):
+        """Returns the species this object refers to."""
+        return self._species
 
-            # Set initial value
-            expr, expr_in_amount = species.initial_value()
-            if expr is not None:
-                # Need to convert initial value if initial value is provided
-                # in concentration
-                if expr_in_amount is None:
-                    expr_in_amount = species.is_amount()
+    def sid(self):
+        """Returns this species reference's SId, or ``None`` if not set."""
+        return self._sid
 
-                try:
-                    if expr_in_amount is False:
-                        # Get initial compartment size
-                        compartment = species.compartment()
-                        size = variable_references[compartment.sid()]
 
-                        # Convert initial value from concentration to amount
-                        expr = myokit.Multiply(expr, myokit.Name(size))
+def convert_name(name):
+    """
+    Converts SBML name to myokit compatible name.
 
-                    # Set initial value
-                    var.set_rhs(expr.clone(subst=expression_references))
-                except AttributeError:
-                    raise SBMLError(
-                        'Initial value of <' + str(species) + '> contains '
-                        'unreferenced parameters/variables. Please use e.g. '
-                        'the `add_parameter` method to reference expressions.')
+    Unlike the SBML SId type, Myokit names cannot start with an
+    underscore.
+    """
+    if name[:1] == '_':
+        name = 'underscore' + name
+    return name
 
-            if species.is_rate():
-                # Get initial state
-                try:
-                    state_value = var.eval()
-                except AttributeError:
-                    state_value = 0
-                    warnings.warn(
-                        'Species <' + str(species) + '> is promoted to state '
-                        'variable without being assigned with an initial '
-                        'value. Default is set to 0.')
 
-                # Promote size to state variable
-                var.promote(state_value=state_value)
+class _MyokitConverter(object):
+    """
+    Converts SBML Models to Myokit models.
+    """
+    @staticmethod
+    def convert(sbml_model):
+        """Converts the given SBML model to a Myokit model."""
 
-            # Set RHS (reactions are dealt with elsewhere)
-            expr = species.value()
-            if expr is not None:
-                # Need to convert initial value if species is measured in
-                # concentration (assignments match unit of measurement)
-                try:
-                    if not species.is_amount():
-                        # Get initial compartment size
-                        compartment = species.compartment()
-                        size = variable_references[compartment.sid()]
+        # Create myokit model
+        myokit_model = myokit.Model(sbml_model.name())
 
-                        # Convert initial value from concentration to amount
-                        expr = myokit.Multiply(expr, myokit.Name(size))
+        # Add notes
+        notes = sbml_model.notes()
+        if notes:
+            myokit_model.meta['desc'] = notes
 
-                    # Set initial value
-                    var.set_rhs(expr.clone(subst=expression_references))
-                except AttributeError:
-                    raise SBMLError(
-                        'Value of <' + str(species) + '> contains '
-                        'unreferenced parameters/variables. Please use e.g. '
-                        'the `add_parameter` method to reference expressions.')
+        # Create reference container that maps sid's to myokit objects
+        component_references = {}
+        variable_references = {}
 
-    def _set_rhs_parameters(self, variable_references, expression_references):
-        """
-        Sets right hand side of global parameters.
-        """
-
-        for sid, parameter in self._parameters.items():
-            # Get myokit variable
-            var = variable_references[sid]
+        # Create additional container for species variables in amount
+        # (Reactions define rates exclusively for amount of species)
+        species_amount_references = {}
 
-            # Set initial value
-            expr = parameter.initial_value()
-            if expr is not None:
-                try:
-                    var.set_rhs(expr.clone(
-                        subst=expression_references))
-                except AttributeError:
-                    raise SBMLError(
-                        'Initial value of <' + str(parameter) +
-                        '> contains unreferenced parameters/variables. Please'
-                        ' use e.g. the `add_parameter` method to add reference'
-                        ' to parameters in the model.')
+        # Create reference from Model expressions to myokit model expressions
+        expression_references = {}
 
-            if parameter.is_rate():
-                # Get initial state
-                try:
-                    state_value = var.eval()
-                except AttributeError:
-                    state_value = 0
-                    warnings.warn(
-                        'Parameter <' + str(parameter) + '> is promoted to '
-                        'state variable without being assigned with an initial'
-                        ' value. Default is set to 1.')
+        # Instantiate component and variable objects first without assigning
+        # RHS expressions. Myokit objects may have to be renamed, so
+        # expressions are added in a second step.
 
-                # Promote size to state variable
-                var.promote(state_value=state_value)
+        # Add SBML compartments to myokit model
+        _MyokitConverter.add_compartments(
+            sbml_model, myokit_model, component_references,
+            variable_references, expression_references)
 
-            # Set RHS
-            # (assignmentRule overwrites initialAssignment)
-            expr = parameter.value()
-            if expr is not None:
-                try:
-                    var.set_rhs(expr.clone(
-                        subst=expression_references))
-                except AttributeError:
-                    raise SBMLError(
-                        'Value of <' + str(parameter) + '> contains '
-                        'unreferenced parameters/variables. Please use e.g.'
-                        ' the `add_parameter` method to add reference '
-                        'to parameters in the model.')
+        # Add global component to model, to store time variable and any global
+        # parameters. Do this after adding the ordinary compartments, so that
+        # they will get precedence in any naming conflicts.
+        global_component = _MyokitConverter.add_global_component(
+            myokit_model, component_references)
 
-    def _set_rhs_stoichiometries(
-            self, variable_references, expression_references):
-        """
-        Sets right hand side of stoichiometry variables.
-        """
+        # Add time variable to myokit component
+        _MyokitConverter.add_time_variable(
+            sbml_model, global_component, expression_references)
 
-        for reaction in self._reactions.values():
-            # Get all reactants and products (modifier do not have
-            # stoichiometries)
-            species_references = reaction.reactants() + reaction.products()
+        # Add species to components
+        _MyokitConverter.add_species(
+            sbml_model, component_references, species_amount_references,
+            variable_references, expression_references)
 
-            for species_reference in species_references:
-                # Get sid
-                sid = species_reference.sid()
+        # Add parameters to myokit component
+        _MyokitConverter.add_parameters(
+            sbml_model, global_component,
+            variable_references, expression_references)
 
-                # If sid does not exist, there is no rhs to set
-                if sid is None:
-                    continue
+        # Add stoichiometries from reactions
+        _MyokitConverter.add_stoichiometries(
+            sbml_model, component_references, variable_references,
+            expression_references)
 
-                # Get stoichiometry variable
-                var = variable_references[sid]
+        # Set RHS of compartment sizes
+        _MyokitConverter.set_rhs_sizes(
+            sbml_model, variable_references, expression_references)
 
-                # Set initial value
-                expr = species_reference.initial_value()
-                if expr is not None:
-                    try:
-                        var.set_rhs(expr.clone(
-                            subst=expression_references))
-                    except AttributeError:
-                        raise SBMLError(
-                            'Initial value of <' + str(species_reference) +
-                            '> (initial stoichiometry) contains unreferenced '
-                            'parameters/variables. Please use e.g. the '
-                            '`add_parameter` method to add reference to '
-                            'parameters in the model.')
+        # Set RHS of species (initalAssignemnt, assignmentRule, rateRule)
+        _MyokitConverter.set_rhs_species(
+            sbml_model, species_amount_references, variable_references,
+            expression_references)
 
-                if species_reference.is_rate():
-                    # Get initial state
-                    try:
-                        state_value = var.eval()
-                    except AttributeError:
-                        state_value = 0
-                        warnings.warn(
-                            'Stoichiometry of <' + str(species_reference) + '>'
-                            ' is promoted to state variable without being '
-                            'assigned with an initial value. Default is set to'
-                            ' 1.')
+        # Set RHS of parameters
+        _MyokitConverter.set_rhs_parameters(
+            sbml_model, variable_references, expression_references)
 
-                    # Promote size to state variable
-                    var.promote(state_value=state_value)
+        # Set RHS of stoichiometries
+        _MyokitConverter.set_rhs_stoichiometries(
+            sbml_model, variable_references, expression_references)
 
-                # Set RHS
-                # (assignmentRule overwrites initialAssignment)
-                expr = species_reference.value()
-                if expr is not None:
-                    try:
-                        var.set_rhs(expr.clone(
-                            subst=expression_references))
-                    except AttributeError:
-                        raise SBMLError(
-                            'Value of <' + str(species_reference) + '> '
-                            '(stroichiometry) contains unreferenced parameters'
-                            '/variables. Please use e.g. the `add_parameter` '
-                            'method to add reference to parameters in the '
-                            'model.')
+        # Set RHS of species changed by reactions
+        _MyokitConverter.set_reactions(
+            sbml_model, variable_references, species_amount_references,
+            expression_references)
 
-    # SBML base units (except Celsius, because it's not defined in myokit)
-    _base_units = {
-        'ampere': myokit.units.A,
-        'avogadro': myokit.parse_unit('1 (6.02214179e23)'),
-        'becquerel': myokit.units.Bq,
-        'candela': myokit.units.cd,
-        'coulomb': myokit.units.C,
-        'dimensionless': myokit.units.dimensionless,
-        'farad': myokit.units.F,
-        'gram': myokit.units.g,
-        'gray': myokit.units.Gy,
-        'henry': myokit.units.H,
-        'hertz': myokit.units.Hz,
-        'item': myokit.units.dimensionless,  # Myokit does not have item unit
-        'joule': myokit.units.J,
-        'katal': myokit.units.kat,
-        'kelvin': myokit.units.K,
-        'kilogram': myokit.units.kg,
-        'liter': myokit.units.L,
-        'litre': myokit.units.L,
-        'lumen': myokit.units.lm,
-        'lux': myokit.units.lux,
-        'meter': myokit.units.m,
-        'metre': myokit.units.m,
-        'mole': myokit.units.mol,
-        'newton': myokit.units.N,
-        'ohm': myokit.units.ohm,
-        'pascal': myokit.units.Pa,
-        'radian': myokit.units.rad,
-        'second': myokit.units.s,
-        'siemens': myokit.units.S,
-        'sievert': myokit.units.Sv,
-        'steradian': myokit.units.sr,
-        'tesla': myokit.units.T,
-        'volt': myokit.units.V,
-        'watt': myokit.units.W,
-        'weber': myokit.units.Wb,
-    }
+        return myokit_model
 
+    @staticmethod
+    def add_global_component(
+            myokit_model, component_references, name='myokit'):
+        """
+        Creates a component used to store global parameters and the time
+        variable.
 
-class Parameter(Quantity):
-    """
-    Represents a parameter in SBML; to create a parameter use
-    :meth:`model.add_parameter()`.
+        An attempt will be made to use the name given by the parameter
+        ``name``, but if this is not available another similar name will be
+        found (using :meth:`myokit.add_component_allow_renaming`).
 
-    Arguments:
+        Returns the newly added component.
+        """
+        # Add component
+        component = myokit_model.add_component_allow_renaming(name)
 
-    ``model``
-        The model this parameter is in.
-    ``sid``
-        This parameter's SId.
+        # Add myokit component to reference list
+        component_references[component.name()] = component
 
-    """
-    def __init__(self, model, sid):
-        super(Parameter, self).__init__()
+        return component
 
-        if not isinstance(model, Model):
-            raise SBMLError(
-                '<' + str(model) + '> needs to be an instance of '
-                'myokit.formats.sbml.Model.')
+    @staticmethod
+    def add_time_variable(sbml_model, component, expression_references):
+        """
+        Adds time bound variable to the myokit compartment.
+        """
 
-        self._model = model
-        self._sid = str(sid)
-        self._units = None
+        # Add variable
+        var = component.add_variable_allow_renaming('time')
 
-    def set_units(self, units):
-        """Sets this parameters units to the given ``Units``."""
-        if not isinstance(units, myokit.Unit):
-            raise SBMLError(
-                '<' + str(units) + '> needs to be instance of myokit.Unit')
+        # Bind time variable to time in myokit model
+        var.set_binding('time')
 
-        self._units = units
+        # Set time unit and initial value (SBML default: t=0)
+        units = sbml_model.time_units()
+        var.set_unit(units)
+        var.set_rhs(myokit.Number(0, units))
 
-    def sid(self):
-        """Returns this parameter's sid."""
-        return self._sid
+        # Add reference to time variable
+        expression_references[myokit.Name(sbml_model.time())] = \
+            myokit.Name(var)
 
-    def __str__(self):
-        return '<Parameter ' + self._sid + '>'
+    @staticmethod
+    def add_compartments(
+            sbml_model, myokit_model, component_references,
+            variable_references, expression_references):
+        """
+        Creates components for each compartment.
+        """
+        for sid, compartment in sbml_model._compartments.items():
+            # Create component for compartment
+            component = myokit_model.add_component(convert_name(sid))
 
-    def units(self):
-        """Returns the units this parameter is in, or ``None`` if not set."""
-        return self._units
+            # Add component to reference list
+            component_references[sid] = component
 
+            # Add compartment size to component
+            var = component.add_variable_allow_renaming('size')
 
-class Reaction(object):
-    """
-    Represents an SBML reaction; to create a reaction use
-    :meth:`Model.add_reaction()`.
+            # Set unit of size variable
+            var.set_unit(compartment.size_units())
 
-    Arguments:
+            # Add size variable to reference list
+            variable_references[sid] = var
+            expression_references[myokit.Name(compartment)] = myokit.Name(var)
 
-    ``model``
-        The model this reaction is in.
-    ``sid``
-        This reaction's SId.
+    @staticmethod
+    def add_parameters(
+            sbml_model, component, variable_references, expression_references):
+        """
+        Adds global parameters to the myokit component.
+        """
 
-    """
-    def __init__(self, model, sid):
+        for sid, parameter in sbml_model._parameters.items():
+            # Add parameter to component
+            var = component.add_variable_allow_renaming(convert_name(sid))
 
-        if not isinstance(model, Model):
-            raise SBMLError(
-                '<' + str(model) + '> needs to be an instance of '
-                'myokit.formats.sbml.Model.')
+            # Set unit of parameter
+            var.set_unit(parameter.units())
 
-        self._model = model
-        self._sid = str(sid)
+            # Add reference to parameter
+            variable_references[sid] = var
+            expression_references[myokit.Name(parameter)] = myokit.Name(var)
 
-        # Reactants, reaction products, and modifiers: all as SpeciesReference
-        self._reactants = []
-        self._products = []
-        self._modifiers = []
+    @staticmethod
+    def add_species(
+            sbml_model, component_references, species_amount_references,
+            variable_references, expression_references):
+        """
+        Adds amount (and potentially concentration) variables for each species
+        to the respective components.
 
-        # All species involved in this reaction (sid to object)
-        self._species = collections.OrderedDict()
+        RHS of species concentration is defined as species amount / size
+        compartment.
+        """
 
-        # The kinetic law specifying this reaction's rate (if set)
-        self._kinetic_law = None
+        for sid, species in sbml_model._species.items():
+            # Get component from reference list
+            compartment = species.compartment()
+            compartment_sid = compartment.sid()
+            try:
+                component = component_references[compartment_sid]
+            except KeyError:
+                raise SBMLError(
+                    'The <' + str(compartment) + '> for <' + str(species) + '>'
+                    ' is not referenced in model. Please use the '
+                    '`add_compartment` method to reference a compartment.')
 
-    def add_modifier(self, species, sid=None):
-        """Adds a modifier to this reaction and returns the created object."""
-        if not isinstance(species, Species):
-            raise SBMLError(
-                '<' + str(species) + '> needs to be an instance of '
-                'myokit.formats.sbml.Species')
+            # Add species in amount to component
+            # (needed for reactions, even if species is defined in
+            # concentration)
+            var = component.add_variable_allow_renaming(
+                convert_name(sid + '_amount'))
 
-        sid = str(sid) if sid else sid
-        if sid is not None:
-            self._model._register_sid(sid)
-        ref = ModifierSpeciesReference(species, sid)
-        self._modifiers.append(ref)
-        self._species[species.sid()] = species
-        return ref
+            # Set unit of amount
+            var.set_unit(species.substance_units())
 
-    def add_product(self, species, sid=None):
-        """
-        Adds a reaction product to this reaction and returns the created
-        object.
-        """
-        if not isinstance(species, Species):
-            raise SBMLError(
-                '<' + str(species) + '> needs to be an instance of '
-                'myokit.formats.sbml.Species')
+            # Add reference to amount variable
+            species_amount_references[sid] = var
 
-        sid = str(sid) if sid else sid
-        if sid is not None:
-            self._model._register_sid(sid)
-        ref = SpeciesReference(species, sid)
-        self._products.append(ref)
-        self._species[species.sid()] = species
-        if sid is not None:
-            self._model._assignables[sid] = ref
-        return ref
+            # Add species in concentration if measured in concentration
+            if not species.is_amount():
+                # Add species in concentration
+                var = component.add_variable_allow_renaming(
+                    convert_name(sid + '_concentration'))
 
-    def add_reactant(self, species, sid=None):
-        """Adds a reactant to this reaction and returns the created object."""
-        if not isinstance(species, Species):
-            raise SBMLError(
-                '<' + str(species) + '> needs to be an instance of '
-                'myokit.formats.sbml.Species')
+                # Get myokit amount and size variable
+                amount = species_amount_references[sid]
+                size = variable_references[compartment_sid]
 
-        sid = str(sid) if sid else sid
-        if sid is not None:
-            self._model._register_sid(sid)
-        ref = SpeciesReference(species, sid)
-        self._reactants.append(ref)
-        self._species[species.sid()] = species
-        if sid is not None:
-            self._model._assignables[sid] = ref
-        return ref
+                # Set unit of concentration
+                var.set_unit(amount.unit() / size.unit())
 
-    def kinetic_law(self):
-        """
-        Returns the kinetic law set for this reaction, or ``None`` if not set.
-        """
-        return self._kinetic_law
+                # Define RHS of concentration as amount / size
+                rhs = myokit.Divide(
+                    myokit.Name(amount),
+                    myokit.Name(size))
+                var.set_rhs(rhs)
 
-    def modifiers(self):
-        """
-        Returns the list of modifiers used by this reaction (as
-        :class:`SpeciesReference` objects).
-        """
-        return list(self._modifiers)
+            # Add reference to species (either in amount or concentration)
+            variable_references[sid] = var
+            expression_references[myokit.Name(species)] = myokit.Name(var)
 
-    def products(self):
+    @staticmethod
+    def add_stoichiometries(
+            sbml_model, component_references, variable_references,
+            expression_references):
         """
-        Returns the list of products created by this reaction (as
-        :class:`SpeciesReference` objects).
+        Adds stoichiometry parameters to the 'compartment' component to which
+        the species belongs.
         """
-        return list(self._products)
 
-    def reactants(self):
-        """
-        Returns the list of reactions used by this reaction (as
-        :class:`SpeciesReference` objects).
-        """
-        return list(self._reactants)
+        for reaction in sbml_model._reactions.values():
+            # Get all reactants and products (modifier do not have
+            # stoichiometries)
+            species_references = reaction.reactants() + reaction.products()
 
-    def set_kinetic_law(self, expression):
+            for species_reference in species_references:
+                sid = species_reference.sid()
+                if sid is not None:
+                    # Get component
+                    species = species_reference.species()
+                    compartment = species.compartment()
+                    compartment_sid = compartment.sid()
+                    try:
+                        component = component_references[compartment_sid]
+                    except KeyError:
+                        raise SBMLError(
+                            'The <' + str(compartment) + '> of <' +
+                            str(species) + '> in <' + str(reaction) + '> is '
+                            'not referenced in the model. Please use the '
+                            '`add_compartment` method to reference the '
+                            'compartment.')
+
+                    # Add variable to component
+                    var = component.add_variable_allow_renaming(
+                        convert_name(sid))
+
+                    # Set unit of variable
+                    # (SBML defines stoichiometries as dimensionless)
+                    var.set_unit(myokit.units.dimensionless)
+
+                    # Add reference to variable
+                    variable_references[sid] = var
+                    expression_references[
+                        myokit.Name(species_reference)] = myokit.Name(var)
+
+    @staticmethod
+    def set_reactions(
+            sbml_model, variable_references, species_amount_references,
+            expression_references):
+
+        for reaction in sbml_model._reactions.values():
+            if reaction.kinetic_law() is None:
+                # Skip to the next reaction
+                continue
+
+            # Set right hand side of reactants
+            _MyokitConverter.set_rhs_reactants(
+                reaction, variable_references, species_amount_references,
+                expression_references)
+
+            # Set right hand side of products
+            _MyokitConverter.set_rhs_products(
+                reaction, variable_references, species_amount_references,
+                expression_references)
+
+    @staticmethod
+    def set_rhs_sizes(sbml_model, variable_references, expression_references):
         """
-        Sets this reaction's kinetic law (as a :class:`myokit.Expression`).
+        Sets right hand side of compartments' size variables.
         """
-        if not isinstance(expression, myokit.Expression):
-            raise SBMLError(
-                '<' + str(expression) + '> needs to be an instance of '
-                'myokit.Expression.')
 
-        self._kinetic_law = expression
+        for sid, compartment in sbml_model._compartments.items():
+            # Get myokit variable
+            var = variable_references[sid]
 
-    def sid(self):
-        """Returns this reaction's sid."""
-        return self._sid
+            # Set initial value
+            expr = compartment.initial_value()
+            if expr is not None:
+                expr = expr.clone(subst=expression_references)
+                try:
+                    var.set_rhs(expr)
+                except AttributeError:
+                    raise SBMLError(
+                        'Initial value for the size of <' + str(compartment) +
+                        '> contains unreferenced parameters/variables. Please'
+                        ' use e.g. the `add_parameter` method to add reference'
+                        ' to parameters in the model.')
 
-    def species(self, sid):
+            if compartment.is_rate():
+                # Get initial state
+                try:
+                    state_value = var.eval()
+                except AttributeError:
+                    state_value = 1
+                    warnings.warn(
+                        'Size of compartment <' + str(compartment) + '> is '
+                        'promoted to state variable without being assigned '
+                        'with an initial value. Default is set to 1.')
+
+                # Promote size to state variable
+                var.promote(state_value=state_value)
+
+            # Set RHS
+            # (assignmentRule overwrites initialAssignment)
+            expr = compartment.value()
+            if expr is not None:
+                expr = expr.clone(subst=expression_references)
+                try:
+                    var.set_rhs(expr)
+                except AttributeError:
+                    raise SBMLError(
+                        'Value for the size of <' + str(compartment) +
+                        '> contains unreferenced parameters/variables. Please'
+                        ' use e.g. the `add_parameter` method to add reference'
+                        ' to parameters in the model.')
+
+    @staticmethod
+    def set_rhs_reactants(
+            reaction, variable_references, species_amount_references,
+            expression_references):
         """
-        Finds and returns a :class:`Species` used in this reaction.
+        Sets right hand side of species acting as reactants in a reaction.
         """
-        return self._species[sid]
 
-    def __str__(self):
-        return '<Reaction ' + self._sid + '>'
+        for reactant in reaction.reactants():
+            # Get species object
+            species = reactant.species()
 
+            if species.is_constant() or species.is_boundary():
+                # Species is not altered by the reaction
+                # Skip to the next reactant
+                continue
 
-class Species(Quantity):
-    """
-    Represents an SBML species; to create a species use
-    :meth:`Compartment.add_species()`.
+            # Instantiate rate expression
+            expr = reaction.kinetic_law().clone()
 
-    Arguments:
+            # Get stoichiometry of reactant
+            try:
+                stoichiometry = myokit.Name(
+                    variable_references[reactant.sid()])
+            except KeyError:
+                stoichiometry = reactant.initial_value()
 
-    ``compartment``
-        The :class:`Compartment` that this species is in.
-    ``sid``
-        This species's SId.
-    ``is_amount``
-        Whether this species value is represented as an amount (if false, it is
-        represented as a concentration, which depends on the size of the
-        compartment it is in).
-    ``is_constant``
-        Whether or not this species is constant.
-    ``is_boundary``
-        Whether or not this species is at the boundary of a reaction.
+            if stoichiometry is not None:
+                # Weight rate expression by stoichiometry
+                expr = myokit.Multiply(stoichiometry, expr)
 
-    """
-    def __init__(self, compartment, sid, is_amount, is_constant, is_boundary):
-        super(Species, self).__init__()
+            factor = species.conversion_factor()
+            if factor is not None:
+                # Get conversion factor variable
+                sid = factor.sid()
+                try:
+                    conversion_factor = variable_references[sid]
+                except KeyError:
+                    raise SBMLError(
+                        'Species <' + str(species) + '> has conversion factor '
+                        '<' + str(factor) + '> which is not referenced in the '
+                        'model.')
 
-        if not isinstance(compartment, Compartment):
-            raise SBMLError(
-                '<' + compartment + '> needs to be instance of'
-                'myokit.formats.sbml.Compartment')
-        if not isinstance(is_amount, bool):
-            raise SBMLError(
-                'Is_amount <' + str(is_amount) + '> needs to be a boolean.')
-        if not isinstance(is_constant, bool):
-            raise SBMLError(
-                'Is_constant <' + str(is_constant) + '> needs to be a boolean.'
-            )
-        if not isinstance(is_boundary, bool):
-            raise SBMLError(
-                'Is_boundary <' + str(is_boundary) + '> needs to be a boolean.'
-            )
-        self._compartment = compartment
+                # Convert rate expression from units of reaction extent
+                # to amount units
+                conversion_factor = myokit.Name(conversion_factor)
+                expr = myokit.Multiply(conversion_factor, expr)
 
-        self._sid = str(sid)
-        self._is_amount = bool(is_amount)
-        self._is_constant = bool(is_constant)
-        self._is_boundary = bool(is_boundary)
+            # Get myokit amount variable
+            try:
+                var = species_amount_references[species.sid()]
+            except KeyError:
+                raise SBMLError(
+                    'Kinetic law of <' + str(reaction) + '> contains '
+                    'unreferenced species <' + str(species) + '>. Please, '
+                    'use the `add_species` method to reference species.')
 
-        # Units for the amount, not the concentration
-        self._units = None
+            if not var.is_state():
+                # Get initial state
+                try:
+                    state_value = var.eval()
+                except AttributeError:
+                    state_value = 0
+                    warnings.warn(
+                        'Species <' + str(species) + '> is promoted to state '
+                        'variable without being assigned with an initial '
+                        'value. Default is set to 0.')
 
-        # Optional conversion factor from substance to extent units
-        self._conversion_factor = None
+                # Promote size to state variable
+                var.promote(state_value=state_value)
+                var.set_rhs(myokit.Number(0))
 
-        # Flag whether initial value is provided in amount or concentration
-        self._initial_value_in_amount = None
+            if var.eval():
+                # Subtract rate contributions
+                # (Reaction removes species from compartment)
+                expr = myokit.Minus(var.rhs(), expr)
+            else:
+                expr = myokit.PrefixMinus(expr)
 
-    def is_amount(self):
+            # Set RHS
+            expr = expr.clone(subst=expression_references)
+            try:
+                var.set_rhs(expr)
+            except AttributeError:
+                raise SBMLError(
+                    'Reaction rate expression of <' + str(reaction) + '> for <'
+                    + str(species) + '> contains unreferenced parameters/'
+                    'variables. Please use e.g. the `add_parameter` method to '
+                    'add reference to parameters in the model.')
+
+    @staticmethod
+    def set_rhs_products(
+            reaction, variable_references, species_amount_references,
+            expression_references):
         """
-        Returns ``True`` only if this species is defined as an amount (not a
-        concentration).
+        Sets right hand side of species acting as products in a reaction.
         """
-        return self._is_amount
 
-    def is_boundary(self):
-        """Returns ``True`` only if this species is at a reaction boundary."""
-        return self._is_boundary
+        for product in reaction.products():
+            # Get species object
+            species = product.species()
 
-    def compartment(self):
-        """Returns the :class:`Compartment` that this species is in."""
-        return self._compartment
+            if species.is_constant() or species.is_boundary():
+                # Species is not altered by the reaction
+                # Skip to the next product
+                continue
 
-    def initial_value(self):
-        """
-        Returns a :class:`myokit.Expression` (or None) for this species'
-        initial value, and a boolean indicating whether initial value is
-        provided in amount (True) or in concentration (False).
-        """
-        return self._initial_value, self._initial_value_in_amount
+            # Instantiate rate expression
+            expr = reaction.kinetic_law().clone()
 
-    def is_constant(self):
-        """Returns ``True`` only if this species is constant."""
-        return self._is_constant
+            # Get stoichiometry of product
+            try:
+                stoichiometry = myokit.Name(
+                    variable_references[product.sid()])
+            except KeyError:
+                stoichiometry = product.initial_value()
+
+            if stoichiometry is not None:
+                # Weight rate expression by stoichiometry
+                expr = myokit.Multiply(stoichiometry, expr)
+
+            factor = species.conversion_factor()
+            if factor is not None:
+                # Get conversion factor variable
+                sid = species.conversion_factor().sid()
+                try:
+                    conversion_factor = variable_references[sid]
+                except KeyError:
+                    raise SBMLError(
+                        'Species <' + str(species) + '> has conversion factor '
+                        '<' + str(factor) + '> which is not referenced in the '
+                        'model.')
+
+                # Convert rate expression from units of reaction extent
+                # to amount units
+                conversion_factor = myokit.Name(conversion_factor)
+                expr = myokit.Multiply(conversion_factor, expr)
+
+            # Get myokit amount variable
+            try:
+                var = species_amount_references[species.sid()]
+            except KeyError:
+                raise SBMLError(
+                    'Kinetic law of <' + str(reaction) + '> contains '
+                    'unreferenced species <' + str(species) + '>. Please, '
+                    'use the `add_species` method to reference species.')
 
-    def conversion_factor(self):
-        """
-        Returns the :class:`Parameter` acting as conversion factor for this
-        species, or the model default if not set (and ``None`` if that isn't
-        set either).
+            if not var.is_state():
+                # Get initial state
+                try:
+                    state_value = var.eval()
+                except AttributeError:
+                    state_value = 0
+                    warnings.warn(
+                        'Species <' + str(species) + '> is promoted to state '
+                        'variable without being assigned with an initial '
+                        'value. Default is set to 0.')
 
-        When calculating the rate of change of the species, multiplying by this
-        factor converts from the species units to "number of reaction events".
-        """
-        if self._conversion_factor is not None:
-            return self._conversion_factor
-        return self._compartment._model.conversion_factor()
+                # Promote size to state variable
+                var.promote(state_value=state_value)
+                var.set_rhs(myokit.Number(0))
 
-    def set_conversion_factor(self, factor):
-        """
-        Sets a :class:`Parameter` as conversion factor for this species,
-        see :meth:`conversion_factor()`.
-        """
-        if not isinstance(factor, Parameter):
-            raise SBMLError(
-                '<' + str(factor) + '> needs to be instance of'
-                'myokit.formats.sbml.Parameter.')
+            if var.rhs().eval():
+                # Add rate contributions
+                expr = myokit.Plus(var.rhs(), expr)
 
-        self._conversion_factor = factor
+            # Set RHS
+            expr = expr.clone(subst=expression_references)
+            try:
+                var.set_rhs(expr)
+            except AttributeError:
+                raise SBMLError(
+                    'Reaction rate expression of <' + str(reaction) + '> for <'
+                    + str(species) + '> contains unreferenced parameters/'
+                    'variables. Please use e.g. the `add_parameter` method to '
+                    'add reference to parameters in the model.')
 
-    def set_initial_value(self, value, in_amount=None):
+    @staticmethod
+    def set_rhs_species(
+            sbml_model, species_amount_references, variable_references,
+            expression_references):
         """
-        Sets a :class:`myokit.Expression` for this species' initial value.
+        Sets right hand side of species amount variables defined by
+        assignments.
 
-        `in_amount` is a boolean that indicates whether the initial value is
-        measured in amount (True) or concentration (False). If set to `None`
-        value is treated to have the same units as the species.
+        Rate expressions defined by reactions are dealt with in
+        :meth:`_set_reactions`.
         """
-        if not isinstance(value, myokit.Expression):
-            raise SBMLError(
-                '<' + str(value) + '> needs to be an instance of '
-                'myokit.Expression.')
-        if (in_amount is not None) and (not isinstance(in_amount, bool)):
-            raise SBMLError(
-                '<in_amount> needs to be an instance of bool or None.')
 
-        self._initial_value = value
-        self._initial_value_in_amount = in_amount
+        for sid, species in sbml_model._species.items():
+            # Get myokit variable
+            # We only adapt amount of species
+            var = species_amount_references[sid]
 
-    def set_substance_units(self, units):
-        """Sets the units this species amount (not concentration) is in."""
-        if not isinstance(units, myokit.Unit):
-            raise SBMLError(
-                '<' + str(units) + '> needs to be instance of myokit.Unit')
+            # Set initial value
+            expr, expr_in_amount = species.initial_value()
+            if expr is not None:
+                # Need to convert initial value if initial value is provided
+                # in concentration
+                if expr_in_amount is None:
+                    expr_in_amount = species.is_amount()
 
-        self._units = units
+                if expr_in_amount is False:
+                    # Get initial compartment size
+                    compartment = species.compartment()
+                    size = variable_references[compartment.sid()]
 
-    def sid(self):
-        """Returns this species's sid."""
-        return self._sid
+                    # Convert initial value from concentration to amount
+                    expr = myokit.Multiply(expr, myokit.Name(size))
 
-    def __str__(self):
-        return '<Species ' + self._sid + '>'
+                # Set initial value
+                expr = expr.clone(subst=expression_references)
+                try:
+                    var.set_rhs(expr)
+                except AttributeError:
+                    raise SBMLError(
+                        'Initial value of <' + str(species) + '> contains '
+                        'unreferenced parameters/variables. Please use e.g. '
+                        'the `add_parameter` method to reference expressions.')
 
-    def substance_units(self):
-        """
-        Returns the units an amount of this species (not a concentration) is
-        in, or the model default if not set.
-        """
-        if self._units is not None:
-            return self._units
-        return self._compartment._model.substance_units()
+            if species.is_rate():
+                # Get initial state
+                try:
+                    state_value = var.eval()
+                except AttributeError:
+                    state_value = 0
+                    warnings.warn(
+                        'Species <' + str(species) + '> is promoted to state '
+                        'variable without being assigned with an initial '
+                        'value. Default is set to 0.')
 
+                # Promote size to state variable
+                var.promote(state_value=state_value)
 
-class SpeciesReference(Quantity):
-    """
-    Represents a reference to a reactant or product species in an SBML
-    reaction.
+            # Set RHS (reactions are dealt with elsewhere)
+            expr = species.value()
+            if expr is not None:
+                # Need to convert initial value if species is measured in
+                # concentration (assignments match unit of measurement)
+                if not species.is_amount():
+                    # Get initial compartment size
+                    compartment = species.compartment()
+                    size = variable_references[compartment.sid()]
 
-    A species reference acts as a :class:`Quantity`, where the value represents
-    the species' stoichiometry.
-    """
-    REACTANT = 0
-    PRODUCT = 1
-    MODIFIER = 2
+                    # Convert initial value from concentration to amount
+                    expr = myokit.Multiply(expr, myokit.Name(size))
 
-    def __init__(self, species, sid=None):
-        super(SpeciesReference, self).__init__()
+                # Set initial value
+                expr = expr.clone(subst=expression_references)
+                try:
+                    var.set_rhs(expr)
+                except AttributeError:
+                    raise SBMLError(
+                        'Value of <' + str(species) + '> contains '
+                        'unreferenced parameters/variables. Please use e.g. '
+                        'the `add_parameter` method to reference expressions.')
 
-        if not isinstance(species, Species):
-            raise SBMLError(
-                '<' + species + '> needs to be instance of'
-                'myokit.formats.sbml.Species')
+    @staticmethod
+    def set_rhs_parameters(
+            sbml_model, variable_references, expression_references):
+        """
+        Sets right hand side of global parameters.
+        """
 
-        self._species = species
-        self._sid = str(sid) if sid else sid
+        for sid, parameter in sbml_model._parameters.items():
+            # Get myokit variable
+            var = variable_references[sid]
 
-    def species(self):
-        """Returns the species this object refers to."""
-        return self._species
+            # Set initial value
+            expr = parameter.initial_value()
+            if expr is not None:
+                expr = expr.clone(subst=expression_references)
+                try:
+                    var.set_rhs(expr)
+                except AttributeError:
+                    raise SBMLError(
+                        'Initial value of ' + str(parameter) +
+                        ' contains unreferenced parameters/variables. Please'
+                        ' use e.g. the `add_parameter` method to add reference'
+                        ' to parameters in the model.')
 
-    def sid(self):
-        """Returns this species reference's SId, or ``None`` if not set."""
-        return self._sid
+            if parameter.is_rate():
+                # Get initial state
+                try:
+                    state_value = var.eval()
+                except AttributeError:
+                    state_value = 0
+                    warnings.warn(
+                        'Parameter ' + str(parameter) + ' is promoted to'
+                        ' state variable without being assigned with an'
+                        ' initial value. Default is set to 1.')
 
-    def __str__(self):
-        return '<SpeciesReference ' + self._sid + '>'
+                # Promote size to state variable
+                var.promote(state_value=state_value)
 
+            # Set RHS
+            # (assignmentRule overwrites initialAssignment)
+            expr = parameter.value()
+            if expr is not None:
+                expr = expr.clone(subst=expression_references)
+                try:
+                    var.set_rhs(expr)
+                except AttributeError:
+                    raise SBMLError(
+                        'Value of ' + str(parameter) + ' contains unreferenced'
+                        ' parameters/variables. Please use e.g. the'
+                        ' `add_parameter` method to add reference to'
+                        ' parameters in the model.')
 
-class ModifierSpeciesReference(object):
-    """Represents a reference to a modifier species in an SBML reaction."""
+    @staticmethod
+    def set_rhs_stoichiometries(
+            sbml_model, variable_references, expression_references):
+        """
+        Sets right hand side of stoichiometry variables.
+        """
 
-    def __init__(self, species, sid=None):
-        super(ModifierSpeciesReference, self).__init__()
+        for reaction in sbml_model._reactions.values():
+            # Get all reactants and products (modifier do not have
+            # stoichiometries)
+            species_references = reaction.reactants() + reaction.products()
 
-        if not isinstance(species, Species):
-            raise SBMLError(
-                '<' + species + '> needs to be instance of'
-                'myokit.formats.sbml.Species')
+            for species_reference in species_references:
+                # Get sid
+                sid = species_reference.sid()
 
-        self._species = species
-        self._sid = str(sid) if sid else sid
+                # If sid does not exist, there is no rhs to set
+                if sid is None:
+                    continue
 
-    def species(self):
-        """Returns the species this object refers to."""
-        return self._species
+                # Get stoichiometry variable
+                var = variable_references[sid]
 
-    def sid(self):
-        """Returns this species reference's SId, or ``None`` if not set."""
-        return self._sid
+                # Set initial value
+                expr = species_reference.initial_value()
+                if expr is not None:
+                    expr = expr.clone(subst=expression_references)
+                    try:
+                        var.set_rhs(expr)
+                    except AttributeError:
+                        raise SBMLError(
+                            'Initial value of <' + str(species_reference) +
+                            '> (initial stoichiometry) contains unreferenced '
+                            'parameters/variables. Please use e.g. the '
+                            '`add_parameter` method to add reference to '
+                            'parameters in the model.')
 
+                if species_reference.is_rate():
+                    # Get initial state
+                    try:
+                        state_value = var.eval()
+                    except AttributeError:
+                        state_value = 0
+                        warnings.warn(
+                            'Stoichiometry of <' + str(species_reference) + '>'
+                            ' is promoted to state variable without being '
+                            'assigned with an initial value. Default is set to'
+                            ' 1.')
 
-def convert_name(name):
-    """
-    Converts SBML name to myokit compatible name.
+                    # Promote size to state variable
+                    var.promote(state_value=state_value)
 
-    Unlike the SBML SId type, Myokit names cannot start with an
-    underscore.
-    """
-    if name[:1] == '_':
-        name = 'underscore' + name
-    return name
+                # Set RHS
+                # (assignmentRule overwrites initialAssignment)
+                expr = species_reference.value()
+                if expr is not None:
+                    expr = expr.clone(subst=expression_references)
+                    try:
+                        var.set_rhs(expr)
+                    except AttributeError:
+                        raise SBMLError(
+                            'Value of <' + str(species_reference) + '> '
+                            '(stroichiometry) contains unreferenced parameters'
+                            '/variables. Please use e.g. the `add_parameter` '
+                            'method to add reference to parameters in the '
+                            'model.')
 
diff --git a/myokit/formats/sbml/_parser.py b/myokit/formats/sbml/_parser.py
index 2378ab67c..e6f93796e 100644
--- a/myokit/formats/sbml/_parser.py
+++ b/myokit/formats/sbml/_parser.py
@@ -251,7 +251,7 @@ def _parse_initial_assignment(self, element, model):
             # Parse rule and apply
             var.set_initial_value(_parse_mathml(
                 expr,
-                lambda x, y: myokit.Name(model.assignable(x)),
+                lambda x, y: myokit.Name(model.assignable_or_csymbol(x)),
             ))
 
         except myokit.formats.sbml.SBMLError as e:
@@ -274,7 +274,7 @@ def _parse_kinetic_law(self, element, model, reaction):
         # Kinetic law is not allowed use species unless they're defined in the
         # reaction.
         def name(x, y):
-            obj = model.assignable(x)
+            obj = model.assignable_or_csymbol(x)
             if isinstance(obj, myokit.formats.sbml.Species):
                 obj = reaction.species(x)
             return myokit.Name(obj)
@@ -561,7 +561,7 @@ def _parse_rule(self, element, model, rate=False):
             # Parse rule and apply
             var.set_value(_parse_mathml(
                 expr,
-                lambda x, y: myokit.Name(model.assignable(x)),
+                lambda x, y: myokit.Name(model.assignable_or_csymbol(x)),
             ), rate)
 
         except myokit.formats.sbml.SBMLError as e:


Test patch:
diff --git a/myokit/tests/test_sbml_api.py b/myokit/tests/test_sbml_api.py
index 1c0f9dc9c..e4a1d386f 100755
--- a/myokit/tests/test_sbml_api.py
+++ b/myokit/tests/test_sbml_api.py
@@ -14,6 +14,8 @@
 import myokit.formats
 import myokit.formats.sbml as sbml
 
+from myokit.formats.sbml._api import _MyokitConverter as X
+
 from shared import WarningCollector
 
 # Unit testing in Python 2 and 3
@@ -156,24 +158,41 @@ def test_area_units(self):
         self.assertEqual(model.area_units(), area_units)
 
     def test_assignable(self):
+        # Tests assignable() and assignable_or_csymbol()
 
         model = sbml.Model(name='model')
 
         # Add assignables to model
         c_sid = 'compartment'
-        model.add_compartment(c_sid)
+        c = model.add_compartment(c_sid)
 
         p_sid = 'parameters'
-        model.add_parameter(p_sid)
+        p = model.add_parameter(p_sid)
 
         s_sid = 'species'
         comp = sbml.Compartment(model=model, sid=c_sid)
-        model.add_species(compartment=comp, sid=s_sid)
+        s = model.add_species(compartment=comp, sid=s_sid)
 
         # Check that all assignables are accessible
-        self.assertIsNotNone(model.assignable(c_sid))
-        self.assertIsNotNone(model.assignable(p_sid))
-        self.assertIsNotNone(model.assignable(s_sid))
+        t = 'http://www.sbml.org/sbml/symbols/time'
+        self.assertEqual(model.assignable(c_sid), c)
+        self.assertEqual(model.assignable(p_sid), p)
+        self.assertEqual(model.assignable(s_sid), s)
+        self.assertRaises(KeyError, model.assignable, t)
+
+        # Check that assignables and csymbols can be obtained with
+        # assignable_or_csymbol:
+        self.assertEqual(model.assignable_or_csymbol(c_sid), c)
+        self.assertEqual(model.assignable_or_csymbol(p_sid), p)
+        self.assertEqual(model.assignable_or_csymbol(s_sid), s)
+        self.assertEqual(model.assignable_or_csymbol(t), model.time())
+        self.assertRaises(
+            KeyError, model.assignable,
+            'http://www.sbml.org/sbml/symbols/beer')
+        self.assertRaises(
+            KeyError, model.assignable, 'http://google.com')
+        self.assertRaises(
+            KeyError, model.assignable, 'https://google.com')
 
     def test_base_unit(self):
 
@@ -417,6 +436,15 @@ def test_substance_units(self):
 
         self.assertEqual(model.substance_units(), myokit.units.g)
 
+    def test_time_symbol(self):
+
+        model = sbml.Model(name='model')
+        time = model.time()
+        turl = 'http://www.sbml.org/sbml/symbols/time'
+        self.assertIsInstance(time, sbml.CSymbolVariable)
+        self.assertEqual(str(time), '<CSymbolVariable ' + turl + '>')
+        self.assertEqual(time.definition_url(), turl)
+
     def test_time_units(self):
 
         model = sbml.Model(name='model')
@@ -1260,7 +1288,7 @@ def test_parameters(self):
         m.add_parameter('c')
         m = m.myokit_model()
 
-        # Check that model created parameters in 'myokit' component
+        # Check that model created parameters in 'global' component
         self.assertTrue(m.has_variable('myokit.z'))
         self.assertTrue(m.has_variable('myokit.boat'))
         self.assertTrue(m.has_variable('myokit.c'))
@@ -1450,8 +1478,7 @@ def test_species_initial_value(self):
         mm = m.myokit_model()
         ms = mm.get('comp.spec_1_amount')
         self.assertFalse(ms.is_state())
-        self.assertEqual(
-            ms.rhs().code(), 'myokit.p1')
+        self.assertEqual(ms.rhs().code(), 'myokit.p1')
 
         # Species in concentration: unreferenced parameter
         p2 = sbml.Parameter(m, 'p2')
@@ -2064,15 +2091,15 @@ def test_reaction_stoichiometries(self):
         c_sid = 'compartment'
         c = m.add_compartment(c_sid)
 
-        m._add_compartments(
-            myokit_model, component_references, variable_references,
+        X.add_compartments(
+            m, myokit_model, component_references, variable_references,
             expression_references)
 
         sid = 'species'
         s = m.add_species(c, sid)
 
-        m._add_species(
-            component_references, species_amount_references,
+        X.add_species(
+            m, component_references, species_amount_references,
             variable_references, expression_references)
 
         # Case I: No stoichiometry reference
@@ -2085,8 +2112,9 @@ def test_reaction_stoichiometries(self):
         r.add_product(species=s, sid=p_sid)
 
         # Add stoichiometries to myokit model
-        m._add_stoichiometries(
-            component_references, variable_references, expression_references)
+        X.add_stoichiometries(
+            m, component_references, variable_references,
+            expression_references)
 
         # Check that stoichiometries do not exist in myokit model and are not
         # referenced in variable_references
@@ -2106,8 +2134,9 @@ def test_reaction_stoichiometries(self):
         ps = r.add_product(species=s, sid=p_sid)
 
         # Add stoichiometries to myokit model
-        m._add_stoichiometries(
-            component_references, variable_references, expression_references)
+        X.add_stoichiometries(
+            m, component_references, variable_references,
+            expression_references)
 
         # Check that stoichiometries exist in myokit model and are
         # referenced in variable_references
@@ -2131,8 +2160,8 @@ def test_reaction_stoichiometries(self):
 
         # Check that stoichiometries are referenced in
 
-    def test_time(self):
-        # Tests whether time variable is created properly.
+    def test_time_variable_is_created(self):
+        # Tests whether time variable is created.
 
         m = sbml.Model()
         m.set_time_units(myokit.units.ampere)
@@ -2146,41 +2175,25 @@ def test_time(self):
         self.assertEqual(var.unit(), myokit.units.ampere)
 
         # Check that initial value is set
-        self.assertEqual(var.eval(), 0)
+        self.assertEqual(var.rhs(), myokit.Number(0, myokit.units.ampere))
 
         # Chet that variable is time bound
         self.assertTrue(var.binding(), 'time')
 
-    def test_time_2(self):
-        # Tests whether time bound variable is added properly to myokit
-        # component.
-
-        # Create myokit component
-        myokit_model = myokit.Model()
-        component_references = {}
-        m = sbml.Model(name='model')
-        m._add_myokit_component(myokit_model, component_references)
-        c = component_references['myokit']
-
-        # Create remaining inputs
-        variable_references = {}
-
-        # Add time bound variable to myokit model
-        unit = myokit.units.s
-        m.set_time_units(unit)
-        m._add_time(c, variable_references)
-
-        # Check that time bound variable exists in myokit component
-        self.assertTrue(c.has_variable('time'))
-
-        # Check that variable is referenced under csymbol
-        sid = 'http://www.sbml.org/sbml/symbols/time'
-        self.assertIn(sid, variable_references.keys())
-
-        # Check that variable has correct units and initial value
-        var = variable_references[sid]
-        self.assertEqual(var.rhs(), myokit.Number(0))
-        self.assertEqual(var.unit(), unit)
+    def test_time_variable_is_usable(self):
+        # Tests that time variable can be used in equations
+
+        # Create SBML model
+        s = sbml.Model()
+        s.set_time_units(myokit.units.ms)
+        p = s.add_parameter('param')
+        p.set_value(myokit.Plus(myokit.Number(1), myokit.Name(s.time())))
+        m = s.myokit_model()
+        t = m.get('myokit.time')
+        self.assertEqual(t.unit(), myokit.units.ms)
+        self.assertEqual(
+            m.get('myokit.param').rhs(),
+            myokit.Plus(myokit.Number(1), myokit.Name(m.time())))
 
 
 if __name__ == '__main__':
diff --git a/myokit/tests/test_sbml_parser.py b/myokit/tests/test_sbml_parser.py
index 511546281..5ad82a24c 100755
--- a/myokit/tests/test_sbml_parser.py
+++ b/myokit/tests/test_sbml_parser.py
@@ -275,6 +275,28 @@ def test_parse_constraint(self):
             self.p.parse_string(self.wrap(xml))
         self.assertIn('Ignoring SBML constraints', w.text())
 
+    def test_parse_csymbol_time(self):
+        # Test parsing of the csymbol "time"
+        xml = (
+            '<model timeUnits="second">'
+            ' <listOfParameters>'
+            '  <parameter id="bacterial_count" value="1" constant="false"/>'
+            ' </listOfParameters>'
+            ' <listOfRules>'
+            '  <assignmentRule variable="bacterial_count">'
+            '   <math xmlns="http://www.w3.org/1998/Math/MathML">'
+            '    <apply>'
+            '     <times/>'
+            '     <ci>bacterial_count</ci>'
+            ' <csymbol definitionURL="http://www.sbml.org/sbml/symbols/time"/>'
+            '    </apply>'
+            '   </math>'
+            '  </assignmentRule>'
+            ' </listOfRules>'
+            '</model>')
+
+        self.p.parse_string(self.wrap(xml))
+
     def test_parse_events(self):
         # Test parsing events (these are ignored)