Repository: myokit/myokit (language: python)

Issue:
state variable with an initial condition dependent on another (constant) variable
I'm trying to define a state variable in a myokit model, that has a default state that is an expression containing other constant variables. Is this possible? From the docs I suspect it is not, and the default state must be a concrete value?


Discussion / hints:
if it does have to be a value, would it be feasible to change this so it could be an expression?
At the moment it's not possible.

I'm not sure how hard it would be to implement this. Some things to think about:

- The model API would have to allow it
  - Any existing checks for cyclic dependencies should be triggered when an initial value is evaluated
  - Can one initial value depend on another? If so, need to extend cyclical checks
- Probably want to restrict constants to non-nested ones? E.g. would it be problematic to have x with initial value x.y ?
- Parser would need updating
- Exports would all need checking:
  - code exports would need to ensure state initialisation is after constant initialisation
  - language exports would need to check the target language supports it
- Sims would need checking (same a runnable code exports above)
- CVODES sim would need checking: does this create weird situations? Not sure! For one thing, if an initial value is no longer a literal, then maybe dy/dq would not allow non-literal initial values as q ?

None of those seem insurmountable, but plenty to think about, and creates lots of new situations that are not covered by the existing unit tests.
Would make this ticket easier: https://github.com/myokit/myokit/issues/500

all sounds doable. If you don't mind I'll start having a go at this. I need to dive into how myokit does everything a bit more to understand what to do, so I'll start with that :)
I think for simplicitly we would restrict initial value expressions to only consist of constants and constant variables
Cool! But expect days rather than hours (but not weeks!)
I think you should be able to create branches on this repo? If you don't want to work in a fork. I'm happy to look at draft PRs and such!
Thanks, working on the repo directly would be preferred. I'll start work on this next week and will push a new branch then

Fix patch:
diff --git a/.gitignore b/.gitignore
index 5beb446b2..523b12cb1 100644
--- a/.gitignore
+++ b/.gitignore
@@ -13,6 +13,7 @@ __pycache__
 
 # Virtual environments
 venv
+env
 
 # DS Store
 .DS_Store
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 57c2813d4..1e639d067 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,15 +4,33 @@ This page lists the main changes made to Myokit in each release.
 
 ## Unreleased
 - Added
-  - [#909](https://github.com/MichaelClerx/myokit/pull/909) Myokit is now tested on Python 3.11.
+  - [#898](https://github.com/myokit/myokit/pull/899) Initial conditions for state variables can now be (constant-valued) expressions. Previously, only literals were allowed.
+  - [#918](https://github.com/myokit/myokit/pull/918) Added `Model.initial_values`, which returns the model's initial values as a list of expressions, floats, or equations.
+  - [#918](https://github.com/myokit/myokit/pull/918) Added `Model.set_initial_values`, which sets the model's initial values using a list of expressions or anything accepted by `map_to_state`.
+  - [#918](https://github.com/myokit/myokit/pull/918) Added `Variable.initial_value`, which returns a state variable's initial value as an expression or float.
+  - [#918](https://github.com/myokit/myokit/pull/918) Added `Variable.set_initial_value`, which sets a state variable's initial value using an expression, string expression, or float.
+  - [#918](https://github.com/myokit/myokit/pull/918) Added `IllegalReferenceInInitialValueError` for illegal references made in initial value expressions.
 - Changed
+  - [#898](https://github.com/myokit/myokit/pull/899) `Variable.promote` can now take an expression, a number, or a string (which will be parsed as an expression) for the argument `initial_value`.
+  - [#909](https://github.com/MichaelClerx/myokit/pull/909) Myokit is now tested on Python 3.11.
   - [#909](https://github.com/MichaelClerx/myokit/pull/909) Myokit is no longer tested on Python 3.6.
+  - [#918](https://github.com/myokit/myokit/pull/918) The only argument to `Variable.promote` has been renamed from `state_value` to `initial_value`.
+  - [#918](https://github.com/myokit/myokit/pull/918) `Model.import_component` now preserves the ordering of the imported states.
 - Deprecated
+  - [#918](https://github.com/myokit/myokit/pull/918) `Model.state` is deprecated in favour of `Model.initial_values`.
+  - [#918](https://github.com/myokit/myokit/pull/918) `Model.set_state` is deprecated in favour of `Model.set_initial_values`.
+  - [#918](https://github.com/myokit/myokit/pull/918) `Model.inits` is deprecated.
+  - [#918](https://github.com/myokit/myokit/pull/918) `Variable.state_value` is deprecated in favour of `Variable.initial_value`.
+  - [#918](https://github.com/myokit/myokit/pull/918) `Variable.set_state_value` is deprecated in favour of `Variable.set_initial_value`.
+  - [#918](https://github.com/myokit/myokit/pull/918) `Model.save_state` and `Model.load_state` are deprecated.
+  - [#918](https://github.com/myokit/myokit/pull/918) The `Variable.promote` keyword argument `state_value` is deprecated in favour of `initial_value`.
 - Removed
+  - [#898](https://github.com/myokit/myokit/pull/899) `NonLiteralValueError` is now no longer raised by Myokit and so has been removed.
   - [#921](https://github.com/myokit/myokit/pull/921) Removed deprecated `LegacySimulation` class.
   - [#921](https://github.com/myokit/myokit/pull/921) Removed deprecated `PSimulation` class.
   - [#921](https://github.com/myokit/myokit/pull/921) Removed deprecated `ICSimulation` class.
 - Fixed
+  - [#918](https://github.com/myokit/myokit/pull/918) Fixed error when `var()` was called on a `PartialDerivative` with a `Derivative` as first argument.
   - [#933](https://github.com/myokit/myokit/pull/933) Made `myokit.step` base its error classification (large/small/none) on numerics, not representation.
 
 ## [1.33.9] - 2022-11-08
diff --git a/docs/source/api_core/exceptions.rst b/docs/source/api_core/exceptions.rst
index 77f601888..c76f47ce0 100644
--- a/docs/source/api_core/exceptions.rst
+++ b/docs/source/api_core/exceptions.rst
@@ -43,6 +43,8 @@ Inheriting classes
 
 .. autoclass:: IllegalReferenceError
 
+.. autoclass:: IllegalReferenceInInitialValueError
+
 .. autoclass:: ImportError
 
 .. autoclass:: IncompatibleModelError
@@ -65,8 +67,6 @@ Inheriting classes
 
 .. autoclass:: MissingTimeVariableError
 
-.. autoclass:: NonLiteralValueError
-
 .. autoclass:: NumericalError
 
 .. autoclass:: ParseError
diff --git a/docs/source/api_index/myokit.rst b/docs/source/api_index/myokit.rst
index 4054e8216..c05c342e3 100644
--- a/docs/source/api_index/myokit.rst
+++ b/docs/source/api_index/myokit.rst
@@ -53,6 +53,7 @@ myokit
 - :class:`myokit.If`
 - :class:`myokit.IllegalAliasError`
 - :class:`myokit.IllegalReferenceError`
+- :class:`myokit.IllegalReferenceInInitialValueError`
 - :class:`myokit.ImportError`
 - :class:`myokit.IncompatibleModelError`
 - :class:`myokit.IncompatibleUnitError`
@@ -89,7 +90,6 @@ myokit
 - :class:`myokit.Multiply`
 - :class:`myokit.MyokitError`
 - :class:`myokit.Name`
-- :class:`myokit.NonLiteralValueError`
 - :class:`myokit.Not`
 - :class:`myokit.NotEqual`
 - :class:`myokit.Number`
diff --git a/docs/source/guide/simulations.rst b/docs/source/guide/simulations.rst
index 68b3b676e..42ccdc3a6 100644
--- a/docs/source/guide/simulations.rst
+++ b/docs/source/guide/simulations.rst
@@ -106,7 +106,7 @@ Example::
     s = myokit.Simulation(m, p)
 
     # The simulation time is now zero, the simulation state is the one returned
-    # by m.state()
+    # by m.initial_values(as_floats=True)
 
     log = s.run(1000)
 
diff --git a/docs/source/syntax/model.rst b/docs/source/syntax/model.rst
index 48681d92f..9abb62d78 100644
--- a/docs/source/syntax/model.rst
+++ b/docs/source/syntax/model.rst
@@ -44,8 +44,9 @@ Meta-data can be added to a model using the syntax ``field: value``
 All state variables require an initial value to be specified in the model
 header using the syntax ``component.variable = value``
 
-Initial values can be numbers or expressions, as long as they don't reference
-other variables.
+Initial values can be numbers or expressions. Expressions can make reference 
+to variables as long as they not nested (see below) and are constant in time.
+References must be made using the syntax ``component.variable``.
 
 Example::
 
@@ -55,9 +56,9 @@ Example::
           backslash notation.
     author: Identifies the author of the model implementation
     membrane.V = -84
-    na_fast.m  = 0
-    na_fast.h  = 1.0
-    na_fast.j  = 1.0
+    ina.m  = 0
+    ina.h  = 0.9
+    ina.j  = 1 / ina.parameter_1
 
 Component syntax
 ================
@@ -416,9 +417,9 @@ The following functions are defined:
 In addition, the expression ``dot(x)`` can be used to reference the time
 derivative of state variable ``x``.
 
-Conditional statements
-======================
-Conditional statements can be made using the ``if`` function::
+Conditional statements (if)
+===========================
+Simple conditional statements can be made using the ``if`` function::
 
     x = if(V < -50,
         0.2 * exp((V - 12) / 4.7),
@@ -432,8 +433,8 @@ Which should be read as::
         x = 0.5 * exp((V + 19) / 1.2)
 
 
-Advanced conditional statements
-===============================
+Piecewise conditional statements
+================================
 Conditional statements with more than 1 branch can be made using the
 ``piecewise`` construct::
 
@@ -586,7 +587,7 @@ myocyte::
           """
     # Template functions
     sig(V, Vstar, a, b) = exp(a * (Vstar - V)) / (1 + exp(b * (Vstar - V)))
-    # Initial conditions
+    # Initial values
     membrane.V         = -84.4
     na_fast.m          = 0.0017
     na_fast.h          = 0.98
diff --git a/myokit/__init__.py b/myokit/__init__.py
index e019bdde2..c12f90b05 100644
--- a/myokit/__init__.py
+++ b/myokit/__init__.py
@@ -341,6 +341,7 @@
     GenerationError,
     IllegalAliasError,
     IllegalReferenceError,
+    IllegalReferenceInInitialValueError,
     ImportError,
     IncompatibleModelError,
     IncompatibleUnitError,
@@ -354,7 +355,6 @@
     MissingRhsError,
     MissingTimeVariableError,
     MyokitError,
-    NonLiteralValueError,
     NumericalError,
     ParseError,
     ProtocolEventError,
diff --git a/myokit/__main__.py b/myokit/__main__.py
index 1ab0543d7..c60f942d9 100644
--- a/myokit/__main__.py
+++ b/myokit/__main__.py
@@ -240,7 +240,7 @@ def add_compiler_parser(subparsers):
 
 def debug(source, variable, deps=False):
     """
-    Shows how a single variable is calculated from the initial conditions.
+    Shows how a single variable is calculated from the initial values.
     """
     import myokit
 
@@ -261,8 +261,8 @@ def add_debug_parser(subparsers):
     parser = subparsers.add_parser(
         'debug',
         description='Shows how a single variable is calculated from the '
-                    'initial conditions. The variable\'s equation and value'
-                    ' are displayed, along with the value and formula of any'
+                    ' initial values. The variable\'s equation and value are'
+                    ' displayed, along with the value and formula of any'
                     ' nested variables and the values of all dependencies.',
         help='Shows how a single variable is calculated.',
     )
diff --git a/myokit/_aux.py b/myokit/_aux.py
index 3338f035a..6a0732d7b 100644
--- a/myokit/_aux.py
+++ b/myokit/_aux.py
@@ -281,12 +281,12 @@ def _meta(self, x1, x2):
 
     def _state(self, m1, m2):
         """
-        Compares two models' states.
+        Compares two models' states and initial values.
         """
         s1 = iter([v.qname() for v in m1.states()])
         s2 = iter([v.qname() for v in m2.states()])
-        c1 = m1.state()
-        c2 = m2.state()
+        c1 = m1.initial_values()
+        c2 = m2.initial_values()
         for k, v1 in enumerate(s1):
             try:
                 v2 = next(s2)
@@ -298,7 +298,7 @@ def _state(self, m1, m2):
                     '[x] Mismatched State at position ' + str(k) + ': [1]<'
                     + v1 + '> [2]<' + v2 + '>')
                 continue
-            if c1[k] != c2[k]:
+            if c1[k].code() != c2[k].code():
                 self._write('[x] Mismatched Initial value for <' + v1 + '>')
         n = m2.count_states()
         for k, v2 in enumerate(s2):
@@ -549,7 +549,7 @@ def step(model, initial=None, reference=None, ignore_errors=False):
     """
     # Get initial state
     if initial is None:
-        initial = model.state()
+        initial = model.initial_values(as_floats=True)
 
     # Get evaluation at initial state
     values = model.evaluate_derivatives(
diff --git a/myokit/_bin/example.mmt b/myokit/_bin/example.mmt
index 5190c568c..901e66599 100644
--- a/myokit/_bin/example.mmt
+++ b/myokit/_bin/example.mmt
@@ -3,7 +3,6 @@ name: Luo-Rudy 1991
 desc: """
 Implementation of the Luo-Rudy model I for the ventricular myocyte.
 Based on an updated formulation downloaded from http://rudylab.wustl.edu
-and licensed under a GNU GPL license version 2.
 """
 # Initial conditions
 membrane.V = -84.5286
@@ -19,7 +18,7 @@ ica.Ca_i   = 0.0002
 desc: """
 Membrane potential and stimulus
 """
-C = 1.0 [uF/cm^2]
+C = 1 [uF/cm^2]
     in [uF/cm^2]
 dot(V) = - (1/C) * (i_ion + i_diff + i_stim) in [mV]
     label membrane_potential
diff --git a/myokit/_err.py b/myokit/_err.py
index a58fcdc12..b0b8ba26d 100644
--- a/myokit/_err.py
+++ b/myokit/_err.py
@@ -169,8 +169,19 @@ class IllegalReferenceError(IntegrityError):
     def __init__(self, reference, owner):
         super(IllegalReferenceError, self).__init__(
             'Illegal reference: The referenced variable <' + reference.qname()
-            + '> is outside the scope of <' + owner.qname() + '>.'
-        )
+            + '> is outside the scope of <' + owner.qname() + '>.')
+
+
+class IllegalReferenceInInitialValueError(IllegalReferenceError):
+    """
+    Raised when an illegal reference is made in an initial value.
+
+    The only way this can occur is if the reference is to a nested variable.
+    """
+    def __init__(self, reference, owner):
+        super(IllegalReferenceError, self).__init__(
+            'Illegal reference made in initial value: The referenced variable'
+            ' <' + reference.qname() + '> is nested.')
 
 
 class ImportError(MyokitError):
@@ -292,14 +303,6 @@ def __init__(self):
         super(MissingTimeVariableError, self).__init__(msg)
 
 
-class NonLiteralValueError(IntegrityError):
-    """
-    Raised when a literal value is required but not given.
-
-    *Extends:* :class:`myokit.IntegrityError`
-    """
-
-
 class NumericalError(MyokitError):
     """
     Raised when a numerical error occurs during the evaluation of a myokit
diff --git a/myokit/_expressions.py b/myokit/_expressions.py
index 65a882287..4502b8541 100644
--- a/myokit/_expressions.py
+++ b/myokit/_expressions.py
@@ -1142,7 +1142,7 @@ def rhs(self):
         """See :meth:`LhsExpression.rhs()`."""
         if self._proper:
             if self._value.is_state():
-                return Number(self._value.state_value())
+                return self._value.initial_value()
             elif self._value.lhs() == self:
                 return self._value.rhs()
         return None
@@ -1390,7 +1390,7 @@ def var(self):
         As with time-derivatives, this returns the variable of which a
         derivative is taken (i.e. the dependent variable "y" in "dy/dx").
         """
-        return self._var1._value
+        return self._var1.var()
 
 
 class InitialValue(LhsExpression):
@@ -1410,7 +1410,7 @@ def __init__(self, var):
         super(InitialValue, self).__init__((var, ))
         if not isinstance(var, Name):
             raise IntegrityError(
-                'The first argument to an initial condition must be a variable'
+                'The first argument to an initial value must be a variable'
                 ' name.', self._token)
 
         self._var = var
@@ -1436,7 +1436,7 @@ def _code(self, b, c):
 
     def _diff(self, lhs, idstates):
         raise NotImplementedError(
-            'Partial derivatives of initial conditions are not supported.')
+            'Partial derivatives of initial values are not supported.')
 
     def _eval_unit(self, mode):
         return self._var._eval_unit(mode)
@@ -1453,7 +1453,7 @@ def rhs(self):
         See :meth:`LhsExpression.rhs()`.
 
         The RHS returned in this case will be ``None``, as there is no RHS
-        associated with initial conditions in the model.
+        associated with initial values in the model.
         """
         # Note: This _could_ return a Number(init, var unit) instead...
         return None
@@ -1471,7 +1471,7 @@ def _validate(self, trail):
         var = self._var._value
         if not (isinstance(var, myokit.Variable) and var.is_state()):
             raise IntegrityError(
-                'Initial conditions can only be defined for state variables.',
+                'Initial values can only be defined for state variables.',
                 self._token)
 
 
@@ -2866,7 +2866,7 @@ def _diff(self, lhs, idstates):
             'Conditions do not have partial derivatives.')
 
 
-class PrefixCondition(Condition, PrefixExpression):
+class PrefixCondition(PrefixExpression, Condition):
     """
     Interface for prefix conditions.
 
@@ -2920,7 +2920,7 @@ def _polishb(self, b):
         self._op._polishb(b)
 
 
-class InfixCondition(Condition, InfixExpression):
+class InfixCondition(InfixExpression, Condition):
     """
     Base class for infix expressions.
 
diff --git a/myokit/_io.py b/myokit/_io.py
index 36d7b2312..c268d7e2f 100644
--- a/myokit/_io.py
+++ b/myokit/_io.py
@@ -93,7 +93,7 @@ def load_script(filename):
 
 def load_state(filename, model=None):
     """
-    Loads an initial state from a file in one of the formats specified by
+    Loads a model state from a file in one of the formats specified by
     :func:`myokit.parse_state()`.
 
     If a :class:`Model` is provided the state will be run through
@@ -110,7 +110,8 @@ def load_state(filename, model=None):
 
 def load_state_bin(filename):
     """
-    Loads an initial state from a file in the binary format used by myokit.
+    Loads a model state from a file in the binary format used by Myokit.
+
     See :meth:`save_state_bin` for details.
     """
     filename = os.path.expanduser(filename)
@@ -244,7 +245,7 @@ def save_script(filename, script):
 
 def save_state(filename, state, model=None):
     """
-    Stores the given state in the file at ``filename``.
+    Stores a model state to the path ``filename``.
 
     If no ``model`` is specified ``state`` should be given as a list of
     floating point numbers and will be stored by simply placing each number on
@@ -272,8 +273,8 @@ def save_state(filename, state, model=None):
 
 def save_state_bin(filename, state, precision=myokit.DOUBLE_PRECISION):
     """
-    Stores the given state (or any list of floating point numbers) in the file
-    at ``filename``, using a binary format.
+    Stores a model state (or any given list of floating point numbers) to the
+    path ``filename``, using a binary format.
 
     The used format is a zip file, containing a single entry: ``state_x_y``,
     where ``x`` is the used data type (``d`` or ``f``) and ``y`` is the number
diff --git a/myokit/_model_api.py b/myokit/_model_api.py
index 3100b76ae..7cf9b6220 100644
--- a/myokit/_model_api.py
+++ b/myokit/_model_api.py
@@ -720,9 +720,7 @@ def _resolve(self, name):
         """ See :meth:`VarProvider._resolve(). """
         def sa(name):
             # Suggest alternative
-            m = self.model()
-            (var, sug, msg) = m.suggest_variable(name)
-            return msg
+            return self.model().suggest_variable(name)[2]
 
         # Try resolving as an alias
         try:
@@ -772,7 +770,7 @@ class Model(ObjectWithMeta, VarProvider):
     Variables stored inside components can be accessed using :meth:`get()` or
     :meth:`values()`. Values defined through their derivative make up the
     model state and can be accessed using :meth:`states()`. States have
-    initial values accessible through :meth:`inits()`.
+    initial values accessible through :meth:`initial_values()`.
 
     A model's validity can be checked using :meth:`is_valid()`, which returns
     the latest validation status and :meth:`validate()`, which (re)validates
@@ -797,10 +795,10 @@ def __init__(self, name=None):
         self._components = {}
 
         # The model's state variables
-        self._state = []
+        self._state_vars = []
 
-        # The model's current state (list of floats)
-        self._current_state = []
+        # The model's initial state, as a list of Expressions
+        self._state_init = []
 
         # A dict mapping binding names to variables
         self._bindings = {}
@@ -1087,9 +1085,11 @@ def clone(self):
         for c in self._components.values():
             c._clone1(clone)
 
-        # Clone state
-        for k, v in enumerate(self._state):
-            clone.get(v.qname()).promote(self._current_state[k])
+        # Create states
+        # Note that the order in which promote() is called determines the
+        # state ordering, so this happens here and not in the Variable class.
+        for k, v in enumerate(self._state_vars):
+            clone.get(v.qname()).promote()
 
         # Create mapping of old var references to new references
         var_map = {}
@@ -1110,6 +1110,11 @@ def clone(self):
         for prefix, prepend in self._reserved_uname_prefixes.items():
             clone.reserve_unique_name_prefix(prefix, prepend)
 
+        # Copy initial state expressions
+        for k, v in enumerate(self._state_vars):
+            clone.get(v.qname()).set_initial_value(
+                self._state_init[k].clone(subst=lhs_map))
+
         # Return
         return clone
 
@@ -1141,17 +1146,15 @@ def _code(self, b, t):
         self._code_meta(b, 0)
 
         # Initial state
-        if self._state:
+        if self._state_vars:
             pre = t * TAB
             b.write(pre + '# Initial values\n')
-            names = [eq.lhs.code() for eq in self.inits()]
+            names = [v.qname() for v in self._state_vars]
+            values = [e.code() for e in self._state_init]
             n = max([len(name) for name in names])
-            names = iter(names)
-            for eq in self.inits():
-                name = next(names)
-                b.write(
-                    pre + name + ' ' * (n - len(name)) + ' = ' + eq.rhs.code()
-                    + '\n')
+            for name, value in zip(names, values):
+                b.write(pre + name + ' ' * (n - len(name)) + ' = '
+                        + value + '\n')
             b.write(pre + '\n')
         else:
             # No initial state? Then add newline
@@ -1243,7 +1246,7 @@ def count_states(self):
         """
         Returns the number of state variables in this model.
         """
-        return len(self._state)
+        return len(self._state_vars)
 
     def create_unique_names(self):
         """
@@ -1357,7 +1360,7 @@ def evaluate_derivatives(
         # Insert new state (if required)
         if state is not None:
             new_state = self.map_to_state(state)
-            for state, value in zip(self._state, new_state):
+            for state, value in zip(self._state_vars, new_state):
                 values[myokit.Name(state)] = value
             state = None
 
@@ -1390,7 +1393,7 @@ def evaluate_derivatives(
                     values[eq.lhs] = eq.rhs.eval(values, precision=precision)
 
         # Return calculated state
-        return [values[state.lhs()] for state in self._state]
+        return [values[state.lhs()] for state in self._state_vars]
 
     def eval_state_derivatives(
             self, state=None, inputs=None, precision=myokit.DOUBLE_PRECISION,
@@ -1489,14 +1492,14 @@ def add_dep(lhs):
     def format_state(self, state=None, state2=None,
                      precision=myokit.DOUBLE_PRECISION):
         """
-        Converts the given list of floating point numbers to a string where
-        each line has the format ``<full_qualified_name> = <float_value>``.
+        Converts a sequence of floating point numbers to a string where each
+        line has the format ``<full_qualified_name> = <float_value>``.
 
         Arguments:
 
         ``state=None``
-            The state to show derivatives for. If no state is given the state
-            returned by :meth:`state` is used.
+            The state to display. If no state is given this model's (evaluated)
+            :meth:`<initial_values()>initial values` are used.
         ``state2=None``
             An optional second state, to be shown next to ``state`` for
             comparison.
@@ -1505,19 +1508,19 @@ def format_state(self, state=None, state2=None,
             :meth:`myokit.float.str` when formatting the state values.
 
         """
-        n = len(self._state)
-        if state is not None:
+        n = len(self._state_vars)
+        if state is None:
+            state = self.initial_values(as_floats=True)
+        else:
             if len(state) != n:
                 raise ValueError(
-                    'Argument `state` must be a list of (' + str(n)
+                    'Argument `state` must be a sequence of (' + str(n)
                     + ') floating point numbers.')
-        else:
-            state = self.state()
 
         if state2 is not None:
             if len(state2) != n:
                 raise ValueError(
-                    'Argument `state2` must be a list of (' + str(n)
+                    'Argument `state2` must be a sequence of (' + str(n)
                     + ') floating point numbers.')
 
         out = []
@@ -1541,35 +1544,35 @@ def format_state_derivatives(self, state=None, derivatives=None,
         Like :meth:`format_state` but displays the derivatives along with
         each state's value.
 
-
         Arguments:
 
         ``state=None``
-            The state to display. If no state is given the state returned by
-            :meth:`state` is used.
+            The state to show derivatives for. If no state is given this
+            model's (evaluated) :meth:`<initial_values()>initial values` are
+            used.
         ``derivatives=None``
-            An optional list of evaluated derivatives. If not given, the values
-            will be calculed from ``state`` using :meth:`eval_derivatives()`.
+            An optional list or other sequence of evaluated derivatives. If not
+            given, the values will be calculed from ``state`` using
+            :meth:`eval_derivatives()`.
         ``precision=myokit.DOUBLE_PRECISION``
             An optional precision argument to use when evaluating the state
             derivatives, and to pass into :meth:`myokit.float.str` when
             formatting the state values and derivatives.
 
         """
-        n = len(self._state)
+        n = len(self._state_vars)
         if state is None:
-            state = self.state()
+            state = self.initial_values(as_floats=True)
         elif len(state) != n:
             raise ValueError(
-                'Argument `state` must be a list of (' + str(n)
+                'Argument `state` must be a sequence of (' + str(n)
                 + ') floating point numbers.')
 
         if derivatives is None:
-            derivatives = self.evaluate_derivatives(
-                state, precision=precision)
+            derivatives = self.evaluate_derivatives(state, precision=precision)
         elif len(derivatives) != n:
             raise ValueError(
-                'Argument `deriv` must be a list of (' + str(n)
+                'Argument `derivatives` must be a sequence of (' + str(n)
                 + ') floating point numbers.')
 
         out = []
@@ -1802,6 +1805,9 @@ def import_component(self, external_component, new_name=None, var_map=None,
             for var in comp.variables():
                 vars_ref.update(var.refs_to(state_refs=False))
                 vars_ref.update(var.refs_to(state_refs=True))
+                if var.is_state():
+                    vars_ref.update(
+                        [e.var() for e in var.initial_value().references()])
             vars_ref.update(comp._alias_map.values())
             vars_ref -= set(comp.variables())
             vars_ref = [x for x in vars_ref if not x.is_nested()]
@@ -1877,8 +1883,8 @@ def import_component(self, external_component, new_name=None, var_map=None,
                 if ext_var in vars_to_map:
                     var_map[ext_var] = self_var
 
-        # add variables to var_map that map to other imported components
-        # but will be reassigned to the clone later
+        # Add variables to var_map that map to other imported components but
+        # will be reassigned to the clone later
         for l in map_to_clone:
             for ext_var in l:
                 if ext_var not in var_map:
@@ -1963,19 +1969,29 @@ def import_component(self, external_component, new_name=None, var_map=None,
                             ' mismatch: ' + str(e))
 
         # Clone component pt 1: create, meta data, empty variables
-        new_component = []
+        new_component = []  # List of components
         for i, comp in enumerate(external_component):
             new_component.append(comp._clone1(self, new_name[i]))
 
-            for var in comp.variables(state=True):
-                # Clone states
-                # TODO: Not sure why clone() code doesn't do this?
-                new_component[i].get(var.qname(comp)).promote(
-                    var.state_value())
             # Now we can add variable to var_map if needed
             for var in map_to_clone[i]:
                 var_map[var] = new_component[i].get(var.qname(comp))
 
+        # Clone states, preserving the state order
+        # Note: The order in which promote() is called determines the order
+        # of the states in the new component. (This is one of the reasons that
+        # the component._clone1 method called above doesn't call promote.)
+        new_states = []     # New states, from all components
+        state_map = {}      # New-state to old state
+        for old_comp, comp in zip(external_component, new_component):
+            for ext_var in old_comp.variables(state=True):
+                var = comp.get(ext_var.name())
+                state_map[var] = ext_var
+                new_states.append(var)
+        new_states.sort(key=lambda var: state_map[var].indice())
+        for var in new_states:
+            var.promote()   # Initial value is set later
+
         # Create mapping of old var references to new references
         # This is a mapping from Name(var) and Derivative(Name(var)) objects
         # to new myokit.Expressions referencing the target model's variables
@@ -1991,9 +2007,10 @@ def import_component(self, external_component, new_name=None, var_map=None,
                     lhs_map[myokit.Derivative(myokit.Name(ext_var))] = \
                         myokit.Derivative(myokit.Name(self_var))
 
-        # Next, add all entries in the var_map. If unit conversion is enabled,
-        # this may include the addition of unit conversion factors
-
+        # Next, add all entries in the var_map to the lhs_map. If unit
+        # conversion is enabled, this may include the addition of unit
+        # conversion factors (so some Names in lhs_map will be mapped onto
+        # Multiply expressions).
         for ext_var, self_var in var_map.items():
             # Substitute in either a reference to self_var, or an expression
             # that converts self_var to the units ext_var's equation expects.
@@ -2022,6 +2039,11 @@ def import_component(self, external_component, new_name=None, var_map=None,
         for i, comp in enumerate(external_component):
             comp._clone2(new_component[i], lhs_map, var_map)
 
+        # Clone initial values
+        for var in new_states:
+            var.set_initial_value(
+                state_map[var].initial_value().clone(subst=lhs_map))
+
         # Time unit conversion? Then update all derivatives.
         if time_factor is not None:
             for comp in new_component:
@@ -2031,15 +2053,32 @@ def import_component(self, external_component, new_name=None, var_map=None,
                         var.set_rhs(
                             myokit.Multiply(rhs, myokit.Number(time_factor)))
 
+    def initial_values(self, as_floats=False):
+        """
+        Returns a list of the model's initial values.
+
+        By default, expressions are returned, but this can be changed to
+        a list of floats by setting ``as_floats=True``.
+        """
+        if as_floats:
+            if any(not e.is_literal() for e in self._state_init):
+                self.validate()  # Check for cycles before evaluating
+            return [float(y) for y in self._state_init]
+        return list(self._state_init)
+
     def inits(self):
         """
-        Returns an iterator over the ``Equation`` objects defining this model's
-        current state.
+        Deprecated method: Returns an iterator over the ``Equation`` objects
+        defining this model's initial values.
         """
+        # Deprecated since 2023-02-22
+        import warnings
+        warnings.warn('The method `inits` is deprecated.')
+
         def StateDefIterator(model):
-            for k, var in enumerate(model._state):
-                yield Equation(
-                    myokit.Name(var), myokit.Number(self._current_state[k]))
+            for var, value in zip(model._state_vars, model._state_init):
+                yield Equation(myokit.Name(var), value)
+
         return StateDefIterator(self)
 
     def is_similar(self, other, check_unames=False):
@@ -2134,10 +2173,13 @@ def __len__(self):
 
     def load_state(self, filename):
         """
-        Sets the model state using data from a file formatted in any style
-        accepted by :func:`myokit.parse_state`.
+        Deprecated method: Sets the model's initial values using data from a
+        file formatted in any style accepted by :func:`myokit.map_to_state`.
         """
-        self.set_state(myokit.load_state(filename, self))
+        # Deprecated since 2023-02-22
+        import warnings
+        warnings.warn('The method `load_state` is deprecated.')
+        self.set_initial_values(myokit.load_state(filename, self))
 
     def map_component_dependencies(
             self, omit_states=True, omit_constants=False):
@@ -2701,7 +2743,7 @@ def reorder_state(self, order):
 
         This method does not affect the model's validation status.
         """
-        n = len(self._state)
+        n = len(self._state_vars)
         if len(order) != n:
             raise ValueError(
                 'The given list must contain the same number of entries as'
@@ -2719,9 +2761,9 @@ def reorder_state(self, order):
                     'Duplicate entry in order specification: "'
                     + str(v.qname()) + '".')
             state.append(v)
-            current.append(self._current_state[v._indice])
-        self._state = state
-        self._current_state = current
+            current.append(self._state_init[v._indice])
+        self._state_vars = state
+        self._state_init = current
         for k, v in enumerate(state):
             v._indice = k
 
@@ -2771,7 +2813,7 @@ def remove_derivative_references(self):
             time_unit = time.unit()
 
         # Scan all states
-        for state in self._state:
+        for state in self._state_vars:
 
             # Search for references to dot(state)
             refs = list(state.refs_by())
@@ -2842,7 +2884,7 @@ def _reset_indices(self):
         """
         Resets the indices of this model's state variables.
         """
-        for k, v in enumerate(self._state):
+        for k, v in enumerate(self._state_vars):
             v._indice = k
 
     def _reset_validation(self):
@@ -2853,7 +2895,11 @@ def _reset_validation(self):
 
     def _resolve(self, name):
         """ See :meth:`VarProvider._resolve(). """
-        return self.get(name)
+        try:
+            return self.get(name)
+        except KeyError:
+            raise myokit.UnresolvedReferenceError(
+                name, self.suggest_variable(name)[2])
 
     def resolve_interdependent_components(self):
         """
@@ -2887,9 +2933,39 @@ def resolve_interdependent_components(self):
 
     def save_state(self, filename):
         """
-        Saves the model state to a file.
+        Deprecated method: Saves the model state to a file (as floats).
+        """
+        # Deprecated since 2023-02-22
+        import warnings
+        warnings.warn('The method `save_state` is deprecated.')
+
+        return myokit.save_state(
+            filename, self.initial_values(as_floats=True), self)
+
+    def set_initial_values(self, values):
+        """
+        Sets this model's initial values.
+
+        The ``values`` must be specified as either a list of floats,
+        expressions, and/or strings; or as a dict or string in a format
+        accepted by :meth:`map_to_state`.
         """
-        return myokit.save_state(filename, self.state(), self)
+        # Use map to state?
+        if isinstance(values, basestring) or isinstance(values, dict):
+            self._state_init = [
+                myokit.Number(x) for x in self.map_to_state(values)]
+        elif len(values) != len(self._state_vars):
+            raise ValueError('Wrong number of initial values, expecting '
+                             + str(len(self._state_vars)) + '.')
+        else:
+            # Parsing of arguments without making changes, in case it fails.
+            expr = []
+            for var, value in zip(self._state_vars, values):
+                expr.append(var._set_initial_value(value, False))
+
+            # Set all at once, and reset validation status
+            self._state_init = expr
+            self._valid = None
 
     def set_name(self, name=None):
         """
@@ -2914,10 +2990,14 @@ def __setstate__(self, state):
 
     def set_state(self, state):
         """
-        Changes this model's state. Accepts any type of input handled by
-        :meth:`map_to_state`.
+        Deprecated method: use :meth:`set_initial_values` instead.
         """
-        self._current_state = self.map_to_state(state)
+        # Deprecated since 2023-02-22
+        import warnings
+        warnings.warn(
+            'The method `set_state` is deprecated. Please use'
+            ' `set_initial_values` instead.')
+        self.set_initial_values(state)
 
     def set_value(self, qname, value):
         """
@@ -2938,6 +3018,9 @@ def show_evaluation_of(self, var):
         The variable's equation and value are displayed, along with the value
         and formula of any nested variables and the values of all dependencies.
         """
+        # Model must be valid, or cycles can occur
+        self.validate()
+
         def format_float(number):
             s = str(number)
             if len(s) < 10:
@@ -2951,7 +3034,10 @@ def format_float(number):
         # Add initial value
         rhs = var.rhs()
         if var.is_state():
-            out.append('Initial value = ' + str(var.state_value()))
+            value = var.initial_value()
+            out.append('Initial value = ' + value.code())
+            if not isinstance(value, myokit.Number):
+                out.append('              = ' + format_float(value))
             out.append(spacer)
         varname = var.lhs().code()
 
@@ -3237,17 +3323,22 @@ def add_equation(eq_list, done, eq):
 
     def state(self):
         """
-        Returns the current state of the model as a list of floating point
-        numbers.
+        Deprecated method, use
+        :meth:`initial_values(as_floats=True)<initial_values>` instead.
         """
-        return list(self._current_state)
+        # Deprecated since 2023-02-22
+        import warnings
+        warnings.warn(
+            'The method `state` is deprecated. Please use'
+            ' `initial_values(as_floats=True)` instead.')
+        return self.initial_values(as_floats=True)
 
     def states(self):
         """
         Returns an iterator over this model's state :class:`variable
         <myokit.Variable>` objects.
         """
-        return iter(self._state)
+        return iter(self._state_vars)
 
     def suggest_variable(self, name):
         """
@@ -3328,10 +3419,7 @@ def time(self):
         "time". For a valid model, this method always returns a unique
         variable. If no time variable has been declared ``None`` is returned.
         """
-        try:
-            return self._bindings['time']
-        except KeyError:
-            return None
+        return self._bindings.get('time')
 
     def timex(self):
         """
@@ -3399,6 +3487,15 @@ def validate(self, remove_unused_variables=False):
                 'Invalid time variable set. Time variable must be bound to'
                 ' external value "time".')
 
+        # Test initial value expresions
+        n = len(self._state_vars)
+        if n != len(self._state_init):   # pragma: no cover
+            # Cover pragma: This can only happen if there's an API bug
+            self._valid = False
+            raise myokit.IntegrityError(
+                'Initial values list must have same size as state variables'
+                ' list.')
+
         # Validation of components, variables
         for c in self.components():
             if c._parent != self:   # pragma: no cover
@@ -3419,15 +3516,6 @@ def validate(self, remove_unused_variables=False):
                     'Component called <' + c.qname() + '> found at index <'
                     + n + '>.')
 
-        # Test current state values
-        n = len(self._state)
-        if n != len(self._current_state):   # pragma: no cover
-            # Cover pragma: This can only happen if there's an API bug
-            self._valid = False
-            raise myokit.IntegrityError(
-                'Current state values list must have same size as state'
-                ' variables list.')
-
         # Find cycles, warn of unused variables
         self._validate_solvability(remove_unused_variables)
 
@@ -3481,7 +3569,7 @@ def follow(lhs, var, trail, mark_used=True):
 
         # Follow all state variables (unless already visited), all bound
         # variables and all used variables.
-        used = [x for x in self._state]
+        used = [x for x in self._state_vars]
         used += [x for x in self._bindings.values()]
         used += [x for x in self._labels.values()]
 
@@ -3547,7 +3635,7 @@ def _var_info(self, var, spacer=None):
         if not isinstance(var, ModelPart):
             var = self.suggest_variable(var)
             if var[0] is None:
-                if var[1] is None:
+                if var[1] is None:  # pragma: no cover
                     raise Exception(var[2])
                 var = var[1]
                 out.append(
@@ -3885,7 +3973,11 @@ def clamp(self, value=None):
         """
         # Set value
         if value is None:
-            value = self.state_value() if self._is_state else self._rhs.eval()
+            self.model().validate()  # Model must be valid before evaluations
+            if self._is_state:
+                value = self.initial_value(True)
+            else:
+                value = self._rhs.eval()
         else:
             value = float(value)
 
@@ -3942,6 +4034,10 @@ def _clone2(self, v, lhs_map):
         for k in self.variables():
             k._clone2(v[k.name()], lhs_map)
 
+        # Note: initial values are stored inside the model, and the state order
+        # depends on the order in which promoting occurs, so states are only
+        # created by Model.clone()
+
     def _code(self, b, t):
         """
         Create the code for this variable and any child variables.
@@ -4081,7 +4177,12 @@ def convert_unit(self, new_unit, helpers=None):
 
         # For states, update the current/initial value
         if self._is_state:
-            self.set_state_value(self.state_value() * float(fw))
+            # Number? Then just multiply. Else use expression.
+            value = self.initial_value()
+            if isinstance(value, myokit.Number):
+                self.set_initial_value(float(value) * float(fw))
+            else:
+                self.set_initial_value(myokit.Multiply(value, fw))
 
         # Update all references to the variable
         old_ref = myokit.Name(self)
@@ -4159,6 +4260,29 @@ def _delete(self, recursive=False, ignore_siblings=False):
                     ' can not be removed: it is used by ' + ' and '.join(
                         ['<' + v.qname() + '>' for v in refs]) + '.')
 
+        # Third check: Do initial values depend on this variable?
+        # Note that, instead of using a cached set in every variable, this
+        # reference is just checked by scanning all init expressions (which
+        # contain a cached set of references).
+        # Note that we don't optimise by checking if this variable is constant,
+        # as it's possible to create (invalid) models where non-constants are
+        # referenced in initial values (but validate() will pick this up!).
+        refs = set()
+        m = self.model()
+        n = myokit.Name(self)
+        for v, e in zip(m._state_vars, m._state_init):
+            if n in e.references():
+                refs.add(v)
+        if ignore_siblings:
+            # Refs from sibling variables are allowed
+            refs = refs.difference(
+                set([x for x in refs if x.has_ancestor(self._parent)]))
+        if refs:
+            raise myokit.IntegrityError(
+                'Variable <' + self.qname() + '> can not be removed: it is'
+                ' used in the inital value(s) for ' + ' and '.join(
+                    ['<' + v.qname() + '>' for v in refs]) + '.')
+
         # At this point it's OK to delete. Rest of the code makes changes,
         # shouldn't raise errors.
 
@@ -4183,7 +4307,6 @@ def _delete(self, recursive=False, ignore_siblings=False):
             self.set_label(None)
 
             # Remove any aliases
-            m = self.parent(Model)
             for c in m.components():
                 c.remove_aliases_for(self)
 
@@ -4218,10 +4341,10 @@ def demote(self):
         model = self.model()
         try:
             # Remove initial value
-            del model._current_state[self._indice]
+            del model._state_init[self._indice]
 
             # Remove this variable from the state
-            del model._state[self._indice]
+            del model._state_vars[self._indice]
 
             # Set lhs to name expression
             self._lhs = myokit.Name(self)
@@ -4269,6 +4392,25 @@ def indice(self):
             raise Exception('Only state variables have initial values.')
         return self._indice
 
+    def initial_value(self, as_float=False):
+        """
+        Returns a state variable's initial value, or raises an exception when
+        called on a non-state variable.
+
+        By default, a :class:`myokit.Expression` is returned. To evaluate and
+        return a float set ``as_float=True``.
+        """
+        if not self._is_state:
+            raise Exception('Only state variables have initial values.')
+
+        model = self.model()
+        expr = model._state_init[self._indice]
+        if not as_float:
+            return expr
+        if not expr.is_literal():
+            model.validate()
+        return expr.eval()
+
     def is_bound(self):
         """
         Returns ``True`` if a binding label has been added to this variable.
@@ -4341,13 +4483,25 @@ def lhs(self):
         """
         return self._lhs
 
-    def promote(self, state_value=0):
+    def promote(self, initial_value=0, state_value=None):
         """
-        Turns this variable into a state variable with a current state value
-        given by ``state_value``.
+        Turns this variable into a state variable with an initial value given
+        by ``initial_value``.
 
-        This will reset the validation status of the model this variable
-        belongs to.
+        The new ``initial_value`` should be:
+
+        1. A numerical value.
+        2. A :class:`myokit.Expression`.
+        3. A string which can be parsed to a :class:`myokit.Expression`. Any
+           references to variables must be made using their fully qualified
+           names.
+
+        Note that expressions can contain references to non-nested and
+        constant-valued variables (i.e. their right-hand side is either a
+        literal expression or refers only to constants).
+
+        Calling ``promote`` will reset the validation status of the model this
+        variable belongs to.
         """
         if self._indice is not None:
             raise Exception('Variable is already a state variable')
@@ -4357,26 +4511,40 @@ def promote(self, state_value=0):
             raise Exception(
                 'State variables cannot be bound to an external value.')
 
-        # Check state value argument
-        if isinstance(state_value, myokit.Expression):
-            if not state_value.is_literal():
-                raise myokit.NonLiteralValueError(
-                    'Expressions for state values can not contain references'
-                    ' to other variables.')
+        # Deprecated on 2023-02-22
+        if state_value is not None:
+            if initial_value != 0:
+                raise Exception('Deprecated keyword argument `state_value` can'
+                                ' not be used at the same time as its'
+                                ' replacement `initial_value`.')
+            initial_value = state_value
 
+            import warnings
+            warnings.warn('The keyword argument `state_value` is deprecated.'
+                          ' Please use `initial_value` instead.')
+
+        # Handle string and number rhs's
         model = self.model()
+        if not isinstance(initial_value, myokit.Expression):
+            if isinstance(initial_value, basestring):
+                # Expressions are evaluated in model context
+                initial_value = myokit.parse_expression(
+                    initial_value, context=model)
+            elif initial_value is not None:
+                initial_value = myokit.Number(initial_value)
+
         try:
             # Set lhs to derivative expression
             self._lhs = myokit.Derivative(myokit.Name(self))
 
             # Get new indice
-            self._indice = len(model._state)
+            self._indice = len(model._state_vars)
 
             # Add to list of states
-            model._state.append(self)
+            model._state_vars.append(self)
 
-            # Add value to list of current values
-            model._current_state.append(float(state_value))
+            # Add initial_value to list of current values
+            model._state_init.append(initial_value)
 
             # All references to this variable are now considered references to
             # its state value
@@ -4606,6 +4774,42 @@ def set_binding(self, binding):
             # Reset model validation
             model._reset_validation()
 
+    def set_initial_value(self, value):
+        """
+        Sets the initial value of a state variable, or raises an exception if
+        called on a non-state variable.
+
+        The new value can be passed in as an expression, number, or a string
+        (in which case it will be parsed as an expression). Expressions can
+        refer to variables as long as they are not nested and are constant in
+        time. Variable references in strings must be made using fully qualified
+        names (``component.variable``).
+        """
+        if not self._is_state:
+            raise Exception('Only state variables have state values.')
+        self._set_initial_value(value, True)
+
+    def _set_initial_value(self, value, make_the_change):
+        """ Internal version of `set_initial_value`. """
+        # Handle strings and floats
+        model = self.model()
+        if not isinstance(value, myokit.Expression):
+            if isinstance(value, basestring):
+                value = myokit.parse_expression(value, context=model)
+            else:
+                value = myokit.Number(value)
+
+        # Allow internal calls to parse `value` without making a change
+        if not make_the_change:
+            return value
+
+        # Update
+        try:
+            model._state_init[self._indice] = value
+        finally:
+            # Reset model validation, but not the variable cache
+            model._reset_validation()
+
     def set_label(self, label=None):
         """
         Adds a unique ``label`` for this variable, indicated that its value can
@@ -4693,19 +4897,14 @@ def set_rhs(self, rhs):
 
     def set_state_value(self, value):
         """
-        If this variable is a state variable, its current value will be
-        updated. For all other variables this raises an exception.
+        Deprecated method, use :meth:`set_initial_value` instead.
         """
-        if not self._is_state:
-            raise Exception('Only state variables have state values.')
-        model = self.model()
-        if isinstance(value, myokit.Expression):
-            if not value.is_literal():
-                raise myokit.NonLiteralValueError(
-                    'Expressions for state values can not contain references'
-                    ' to other variables.')
-        model._current_state[self._indice] = float(value)
-        # No need to reset validation status or cache here.
+        # Deprecated since 2023-02-22
+        import warnings
+        warnings.warn(
+            'The method `set_state_value` is deprecated. Please use'
+            ' `set_initial_value` instead.')
+        self.set_initial_value(value)
 
     def set_unit(self, unit=None):
         """
@@ -4723,12 +4922,15 @@ def set_unit(self, unit=None):
 
     def state_value(self):
         """
-        For state variables, this will return their current value.
-        For all other variables, this will raise an exception.
+        Deprecated method, use
+        :meth:`initial_value(as_float=True)<initial_value>` instead.
         """
-        if not self._is_state:
-            raise Exception('Only state variables have initial values.')
-        return self.model()._current_state[self._indice]
+        # Deprecated since 2023-02-22
+        import warnings
+        warnings.warn(
+            'The method `state_value` is deprecated. Please use'
+            ' `initial_value(as_float=True)` instead.')
+        return self.initial_value(as_float=True)
 
     def unit(self, mode=myokit.UNIT_TOLERANT):
         """
@@ -4752,46 +4954,89 @@ def validate(self):
         """
         Attempts to check this variable's validity, raises errors if it isn't.
         """
+        #
         # Validate rhs
+        #
         if self._rhs is None:
             raise myokit.MissingRhsError(self)
         self._rhs.validate()
 
-        # Partial derivatives are not allowed in an RHS
+        # RHS: No PartialDerivative objects
         if self._rhs.contains_type(myokit.PartialDerivative):
             raise myokit.IntegrityError(
                 'Partial derivatives may not appear in expressions set as'
-                ' right-hand side of a variable.')
+                ' right-hand side of a variable: <' + self.qname() + '>.')
 
-        # Initial values are not allowed in an RHS
+        # RHS: No InitialValue objects
         if self._rhs.contains_type(myokit.InitialValue):
             raise myokit.IntegrityError(
                 'Initial value operators may not appear in expressions set as'
-                ' right-hand side of a variable.')
+                ' right-hand side of a variable: <' + self.qname() + '>.')
 
-        # Conditions are not allowed as an RHS
+        # RHS: Can't evaluate to True or False
         if isinstance(self._rhs, myokit.Condition):
             raise myokit.IntegrityError(
                 'The right-hand side expression for a variable can not be a'
-                ' condition.')
+                ' condition: <' + self.qname() + '>.')
 
+        #
         # Check state variables
+        #
         is_state = self._indice is not None
         is_deriv = self.lhs().is_derivative()
         if is_state:
+            # Derivative is set
             if not is_deriv:        # pragma: no cover
                 raise myokit.IntegrityError(
                     'Variable <' + self.qname() + '> is listed as a state'
                     ' variable but its lhs is not a derivative.')
+
+            # Not nested
             if self._is_nested:     # pragma: no cover
                 raise myokit.IntegrityError(
                     'State variables should not be nested: <'
                     + str(self.qname()) + '>.')
+
+            # Index matches model
             m = self.model()
-            if not m._state[self._indice] == self:  # pragma: no cover
+            if not m._state_vars[self._indice] == self:  # pragma: no cover
                 raise myokit.IntegrityError(
                     'State variable not listed in model state vector at'
                     ' correct indice: <' + self.qname() + '>.')
+
+            # Initial value is an expression
+            i = m._state_init[self._indice]
+            if not isinstance(i, myokit.Expression):  # pragma: no cover
+                raise myokit.IntegrityError(
+                    'Initial value for <' + self.qname() + '> is not an'
+                    ' expression.')
+
+            # Init: No PartialDerivative or InitialValue operators
+            if i.contains_type(myokit.PartialDerivative):
+                raise myokit.IntegrityError(
+                    'Partial derivatives may not appear in model expressions:'
+                    ' initial value for <' + self.qname() + '>.')
+            if i.contains_type(myokit.InitialValue):
+                raise myokit.IntegrityError(
+                    'Initial values may not appear in model expressions:'
+                    ' initial value for < ' + self.qname() + '>.')
+
+            # Init: Can't evaluate to True or False
+            if isinstance(i, myokit.Condition):
+                raise myokit.IntegrityError(
+                    'The initial value for a variable can not be a'
+                    ' condition: <' + self.qname() + '>.')
+
+            # Init: No nested variables or non-constants
+            for ref in i.references():
+                var = ref.var()
+                if var.is_nested():
+                    raise myokit.IllegalReferenceInInitialValueError(self, ref)
+                if not var.is_constant():
+                    raise myokit.IntegrityError(
+                        'Initial value for variable <' + self.qname() + '> is'
+                        ' not constant: ' + i.code() + '.')
+
         elif is_deriv:  # pragma: no cover
             raise myokit.IntegrityError(
                 'A derivative was set for <' + self.qname() + '> but this is'
diff --git a/myokit/_myokit_version.py b/myokit/_myokit_version.py
index 8cbdaf12e..d95079a7d 100644
--- a/myokit/_myokit_version.py
+++ b/myokit/_myokit_version.py
@@ -15,7 +15,7 @@
 #    incompatibility
 #  - Changes to revision indicate bugfixes, tiny new features
 #  - There is no significance to odd/even numbers
-__version_tuple__ = 1, 33, 10
+__version_tuple__ = 1, 34, 0
 
 # String version of the version number
 __version__ = '.'.join([str(x) for x in __version_tuple__])
diff --git a/myokit/_parsing.py b/myokit/_parsing.py
index 284943f8b..c47ea8246 100644
--- a/myokit/_parsing.py
+++ b/myokit/_parsing.py
@@ -352,6 +352,7 @@ class ParseInfo(object):
     def __init__(self):
         self.model = None
         self.initial_values = OrderedDict()
+        self.initial_value_tokens = {}
         self.alias_map = {}
         self.user_functions = {}
 
@@ -420,10 +421,10 @@ def parse_model_from_stream(stream, syntax_only=False):
                     'Duplicate initial value', t0[2], t0[3],
                     'A value for <' + name + '> was already specified.')
             expect(next(stream), EQUAL)
-            expr = parse_expression_stream(stream)
+            expr = parse_proto_expression(stream)
             expect(next(stream), EOL)
             info.initial_values[name] = expr
-            reg_token(info, t0, expr)
+            info.initial_value_tokens[name] = t0
 
         token = stream.peek()
 
@@ -439,15 +440,6 @@ def parse_model_from_stream(stream, syntax_only=False):
     if syntax_only:
         return True
 
-    # All initial variables must have been used
-    for qname, e in info.initial_values.items():
-        raise ParseError(
-            'Unused initial value', 0, 0,
-            'An unused initial value was found for "' + str(qname) + '".')
-
-    # Re-order the model state
-    model.reorder_state(state_order)
-
     # Order encountered tokens
     m = model._tokens
     model._tokens = {}
@@ -468,6 +460,23 @@ def parse_model_from_stream(stream, syntax_only=False):
             var.set_rhs(convert_proto_expression(var._proto_rhs, var, info))
         del var._proto_rhs
 
+    # Resolve variable references in initial values
+    for i, var in enumerate(model.states()):
+        e = info.initial_values[var.qname()]
+        expr = convert_proto_expression(e, model, info)
+        var.set_initial_value(expr)
+        del info.initial_values[var.qname()]
+
+    # All initial variables must have been used
+    for qname, e in info.initial_values.items():
+        t = info.initial_value_tokens[qname]
+        raise ParseError(
+            'Unused initial value', t[2], t[3],
+            'An unused initial value was found for "' + qname + '".')
+
+    # Re-order the model state
+    model.reorder_state(state_order)
+
     # Check the semantics of the model
     try:
         model.validate()
@@ -675,14 +684,7 @@ def parse_variable(stream, info, parent):
             raise ParseError(
                 'Missing initial value', line, char,
                 'No initial value was found for "' + var.qname() + '"')
-        state_value = info.initial_values[var.qname()]
-        try:
-            var.promote(state_value)
-        except myokit.NonLiteralValueError as e:
-            t = state_value._token
-            raise ParseError(
-                'Illegal state value', t[2], t[3], str(e), cause=e)
-        del info.initial_values[var.qname()]
+        var.promote()
 
     # Parse definition, quick unit, bind, label and description syntax
     # These token must occur in a fixed order!
diff --git a/myokit/_sim/cable.py b/myokit/_sim/cable.py
index 01e7baf8f..4eec6032f 100644
--- a/myokit/_sim/cable.py
+++ b/myokit/_sim/cable.py
@@ -176,7 +176,7 @@ def __init__(self, model, protocol=None, ncells=50, rl=False):
                 ' next')
 
         # Set state and default state
-        self._state = self._model.state() * ncells
+        self._state = self._model.initial_values(True) * ncells
         self._default_state = list(self._state)
 
         # Unique simulation id
diff --git a/myokit/_sim/cmodel.h b/myokit/_sim/cmodel.h
index be2d1e19f..5481f12c7 100644
--- a/myokit/_sim/cmodel.h
+++ b/myokit/_sim/cmodel.h
@@ -690,7 +690,7 @@ if model.binding('pace') is not None:
 }
 
 /*
- * Updates the state variables to the values given in `states`.
+ * Updates the state variables to the float values given in `states`.
  *
  * If any of the values are changed, the model caches are cleared.
  *
@@ -1222,10 +1222,10 @@ for eq in parameters.values():
         return NULL;
     }
 
-    /* States */
+    /* State values */
 <?
 for var in model.states():
-    print(tab + v(var) + ' = ' + myokit.float.str(var.state_value()) + ';')
+    print(tab + v(var) + ' = ' + myokit.float.str(var.initial_value(True)) + ';')
 ?>
     /*
      * Caching.
diff --git a/myokit/_sim/cvodessim.py b/myokit/_sim/cvodessim.py
index a354dbb85..9abbbf2bc 100644
--- a/myokit/_sim/cvodessim.py
+++ b/myokit/_sim/cvodessim.py
@@ -167,6 +167,18 @@ def __init__(self, model, protocol=None, sensitivities=None, path=None):
         cmodel = myokit.CModel(self._model, sensitivities)
         if cmodel.has_sensitivities:
             self._sensitivities = (cmodel.dependents, cmodel.independents)
+
+            # Check for sensitivities w.r.t. variables used in initial state
+            # expressions. This is not implemented yet.
+            inits = self._model.initial_values()
+            for i in self._sensitivities[1]:  # Expressions
+                if isinstance(i, myokit.Name):
+                    for e in inits:
+                        if e.depends_on(i, deep=True):
+                            raise NotImplementedError(
+                                'Sensitivities with respect to parameters used'
+                                ' in initial conditions is not implemented ('
+                                + e.code() + ' depends on ' + i.code() + ').')
         else:
             self._sensitivities = None
 
@@ -186,7 +198,7 @@ def __init__(self, model, protocol=None, sensitivities=None, path=None):
         del cmodel
 
         # Get state and default state from model
-        self._state = self._model.state()
+        self._state = self._model.initial_values(as_floats=True)
         self._default_state = list(self._state)
 
         # Set state and default state for sensitivities
@@ -903,7 +915,7 @@ def set_constant(self, var, value):
 
     def set_default_state(self, state):
         """
-        Allows you to manually set the default state.
+        Change the default state to ``state``.
         """
         self._default_state = self._model.map_to_state(state)
 
diff --git a/myokit/_sim/fiber_tissue.py b/myokit/_sim/fiber_tissue.py
index 1df6070a3..42a1adff7 100644
--- a/myokit/_sim/fiber_tissue.py
+++ b/myokit/_sim/fiber_tissue.py
@@ -303,8 +303,8 @@ def __init__(
                 ' and the tissue model: ' + str(ucf) + ' vs ' + str(uct) + '.')
 
         # Set state and default state
-        self._statef = self._modelf.state() * self._ntotalf
-        self._statet = self._modelt.state() * self._ntotalt
+        self._statef = self._modelf.initial_values(True) * self._ntotalf
+        self._statet = self._modelt.initial_values(True) * self._ntotalt
         self._default_statef = list(self._statef)
         self._default_statet = list(self._statet)
 
diff --git a/myokit/_sim/jacobian.py b/myokit/_sim/jacobian.py
index d61d8289f..6ed80884a 100644
--- a/myokit/_sim/jacobian.py
+++ b/myokit/_sim/jacobian.py
@@ -223,6 +223,10 @@ def __init__(self, model):
         # Clone model
         self._model = model.clone()
 
+        # Store the initial values (won't be able to access once time binding
+        # is removed).
+        self._state = self._model.initial_values(as_floats=True)
+
         # Unbind all inputs
         for label, var in self._model.bindings():
             var.set_binding(None)
@@ -310,7 +314,7 @@ def newton_root(self, x=None, accuracy=0, max_iter=50, damping=1):
 
         # Get initial state
         if x is None:
-            x = self._model.state()
+            x = self._state
         x = np.array(x)
 
         # Calculate derivatives & jacobian
diff --git a/myokit/_sim/opencl.py b/myokit/_sim/opencl.py
index b63307a71..4aa658edb 100644
--- a/myokit/_sim/opencl.py
+++ b/myokit/_sim/opencl.py
@@ -292,7 +292,7 @@ def format(self):
         b = []
         for i, platform in enumerate(self.platforms):
             pname = 'Platform ' + str(i)
-            if platform.name == selected_platform:
+            if platform.name == selected_platform:  # pragma: no cover
                 pname += '  <-- Selected platform in myokit.ini'
             b.append(pname)
             platform._format(b, selected_platform, selected_device, pre=' ')
@@ -363,7 +363,7 @@ def _format(self, b, p_selected=None, d_selected=None, pre='', name=True):
             b.append(pre + 'Devices:')
             for j, device in enumerate(self.devices):
                 dname = pre + ' Device ' + str(j)
-                if self.name == p_selected and device.name == d_selected:
+                if self.name == p_selected and device.name == d_selected:  # pragma: no cover # noqa
                     dname += '  <-- Selected device in myokit.ini'
                 b.append(dname)
                 device._format(b, pre + '  ')
diff --git a/myokit/_sim/openclsim.py b/myokit/_sim/openclsim.py
index 3cece26db..66f609187 100644
--- a/myokit/_sim/openclsim.py
+++ b/myokit/_sim/openclsim.py
@@ -272,7 +272,7 @@ def __init__(
         self._nstate = self._model.count_states()
 
         # Set state and default state
-        self._state = self._model.state() * self._ntotal
+        self._state = self._model.initial_values(True) * self._ntotal
         self._default_state = list(self._state)
 
         # List of globally logged inputs
diff --git a/myokit/formats/ansic/template/cable.c b/myokit/formats/ansic/template/cable.c
index 1c1219aa4..856b1b553 100644
--- a/myokit/formats/ansic/template/cable.c
+++ b/myokit/formats/ansic/template/cable.c
@@ -155,8 +155,8 @@ static void
 Cell_set_initial_state(Cell *cell)
 {
 <?
-for eq in model.inits():
-    print(tab + w.eq(eq) + ';')
+for var, init in zip(model.states(), model.initial_values(True)):
+    print(tab + v(var) + ' = ' + myokit.float.str(init) + ';')
 ?>}
 
 /*
diff --git a/myokit/formats/ansic/template/euler.c b/myokit/formats/ansic/template/euler.c
index 2b915efad..689b8b0ad 100644
--- a/myokit/formats/ansic/template/euler.c
+++ b/myokit/formats/ansic/template/euler.c
@@ -220,7 +220,7 @@ int main()
     Real state[N_STATE];
 <?
 for var in model.states():
-    print(tab + v(var) + ' = ' + myokit.float.str(var.state_value()) + ';')
+    print(tab + v(var) + ' = ' + myokit.float.str(var.initial_value(True)) + ';')
 
 ?>
 
diff --git a/myokit/formats/ansic/template/sim.c b/myokit/formats/ansic/template/sim.c
index 80003c8f0..9b45a8434 100644
--- a/myokit/formats/ansic/template/sim.c
+++ b/myokit/formats/ansic/template/sim.c
@@ -131,8 +131,8 @@ static void
 default_initial_values(N_Vector y)
 {
 <?
-for eq in model.inits():
-    print(tab + w.eq(eq) + ';')
+for var, init in zip(model.states(), model.initial_values(True)):
+    print(tab + v(var) + ' = ' + myokit.float.str(init) + ';')
 ?>
 }
 
diff --git a/myokit/formats/cellml/v1/_api.py b/myokit/formats/cellml/v1/_api.py
index 9f2d3aa15..939b3a61a 100644
--- a/myokit/formats/cellml/v1/_api.py
+++ b/myokit/formats/cellml/v1/_api.py
@@ -898,7 +898,7 @@ def variable_unit(variable, time_unit):
                 # Promote states and set rhs and initial value
                 elif variable.is_state():
                     v.set_is_state(True)
-                    v.set_initial_value(variable.state_value())
+                    v.set_initial_value(variable.initial_value(True))
                     v.set_rhs(rhs)
 
                 # Store literals (single number) in initial value
diff --git a/myokit/formats/cellml/v2/_api.py b/myokit/formats/cellml/v2/_api.py
index 4ec82841e..e6943b355 100644
--- a/myokit/formats/cellml/v2/_api.py
+++ b/myokit/formats/cellml/v2/_api.py
@@ -802,7 +802,7 @@ def variable_unit(variable, time_unit):
 
                 # Promote states and set rhs and initial value
                 elif variable.is_state():
-                    v.set_initial_value(variable.state_value())
+                    v.set_initial_value(variable.initial_value(True))
                     v.set_equation(myokit.Equation(lhs, rhs))
 
                 # Store literals (single number) in initial value
diff --git a/myokit/formats/cuda/template/kernel.cu b/myokit/formats/cuda/template/kernel.cu
index fc004983e..60c73132b 100644
--- a/myokit/formats/cuda/template/kernel.cu
+++ b/myokit/formats/cuda/template/kernel.cu
@@ -66,6 +66,9 @@ model.reserve_unique_names(
     )
 model.create_unique_names()
 
+# Get initial values (before time binding is removed and model becomes invalid)
+initial_values = model.initial_values(True)
+
 # Process bindings, remove unsupported bindings, get map of bound variables to
 # internal names
 bound_variables = model.prepare_bindings({
@@ -295,10 +298,10 @@ int get_default_initial_state(Real *state)
     if (state == 0) return(-1);
 
 <?
-for var in model.states():
+for var, val in zip(model.states(), initial_values):
     if 'desc' in var.meta:
         print(tab + '// ' + var.meta['desc'])
-    print(tab + v(var) + ' = ' + myokit.float.str(var.state_value()) + ';')
+    print(tab + v(var) + ' = ' + myokit.float.str(val) + ';')
 ?>
 
     return(0);
diff --git a/myokit/formats/easyml/_exporter.py b/myokit/formats/easyml/_exporter.py
index 7aa0f29e5..cb9b8f1d0 100644
--- a/myokit/formats/easyml/_exporter.py
+++ b/myokit/formats/easyml/_exporter.py
@@ -372,8 +372,8 @@ def lhs(e):
 
             # Write initial conditions
             for v in model.states():
-                f.write(lhs(v) + '_init = ' + myokit.float.str(v.state_value())
-                        + eos)
+                f.write(lhs(v) + '_init = '
+                        + myokit.float.str(v.initial_value(True)) + eos)
             f.write(eol)
 
             # Write remaining variables
diff --git a/myokit/formats/latex/_exporter.py b/myokit/formats/latex/_exporter.py
index 0fc320309..5d493f009 100644
--- a/myokit/formats/latex/_exporter.py
+++ b/myokit/formats/latex/_exporter.py
@@ -76,7 +76,7 @@ def model(self, path, model, protocol=None):
                 f.write('\\begin{dmath}\n')
                 f.write(e.ex(myokit.Name(v)))
                 f.write(' = ')
-                f.write(e.ex(myokit.Number(v.state_value())))
+                f.write(e.ex(v.initial_value()))
                 f.write('\\end{dmath}\n')
 
             # Write each component
diff --git a/myokit/formats/matlab/template/main.m b/myokit/formats/matlab/template/main.m
index e3be9468a..d1fd47003 100644
--- a/myokit/formats/matlab/template/main.m
+++ b/myokit/formats/matlab/template/main.m
@@ -5,6 +5,7 @@
 # This file is part of Myokit.
 # See http://myokit.org for copyright, sharing, and licensing details.
 #
+import myokit
 ?>%
 % Driver file for <?= model.name() ?>
 %
@@ -15,8 +16,8 @@
 
 % Initial values
 <?
-for eq in model.inits():
-    print(e(eq) + ';')
+for var, value in zip(model.states(), model.initial_values(True)):
+    print(v(var) + ' = ' + myokit.float.str(value) + ';')
 ?>
 
 % Vector for initial values
diff --git a/myokit/formats/opencl/_exporter.py b/myokit/formats/opencl/_exporter.py
index d90116f4a..e5cd10590 100644
--- a/myokit/formats/opencl/_exporter.py
+++ b/myokit/formats/opencl/_exporter.py
@@ -44,7 +44,7 @@ def _dict(self):
             'kernel.cl': 'kernel.cl',
             'plot.py': 'plot.py',
             'minilog.py': 'minilog.py',
-            'test': 'test',
+            'test.sh': 'test.sh',
         }
 
     def _vars(self, model, protocol):
diff --git a/myokit/formats/opencl/template/cable.c b/myokit/formats/opencl/template/cable.c
index d1fea5e10..3f7ffbc40 100644
--- a/myokit/formats/opencl/template/cable.c
+++ b/myokit/formats/opencl/template/cable.c
@@ -192,7 +192,7 @@ void set_initial_values(int n_cells, Real* s)
     for(i=0; i<n_cells; i++) {
 <?
 for v in model.states():
-    print(tab*2 + '*s = ' + myokit.float.str(v.state_value()) + '; s++; // ' + str(v.qname()))
+    print(tab*2 + '*s = ' + myokit.float.str(v.initial_value(True)) + '; s++; // ' + str(v.qname()))
 ?>
     }
 }
@@ -234,14 +234,14 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     cl_device_id device_id = get_device_id();
     char buffer[65536];
     ocl(clGetDeviceInfo(device_id, CL_DEVICE_NAME, sizeof(buffer), buffer, NULL));
-    printf("Using device: %s\n", buffer);
+    // printf("Using device: %s\n", buffer);
 
     /* Create a context and command queue */
     cl_int err;
     context = clCreateContext(NULL, 1, &device_id, NULL, NULL, &err);
     ocl_check("clCreateContext", err);
     command_queue = clCreateCommandQueue(context, device_id, 0, &err);
-    printf("Created context and command queue.\n");
+    // printf("Created context and command queue.\n");
 
     /* Create state vector */
     int dsize_state = n_cells * n_state * sizeof(Real);
@@ -258,12 +258,12 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     ocl_check("clCreateBuffer mbuf_state", err);
     mbuf_idiff = clCreateBuffer(context, CL_MEM_READ_ONLY, dsize_idiff, NULL, &err);
     ocl_check("clCreateBuffer mbuf_idiff", err);
-    printf("Created buffers.\n");
+    // printf("Created buffers.\n");
 
     /* Copy data into buffers */
     ocl(clEnqueueWriteBuffer(command_queue, mbuf_state, CL_TRUE, 0, dsize_state, state, 0, NULL, NULL));
     ocl(clEnqueueWriteBuffer(command_queue, mbuf_idiff, CL_TRUE, 0, dsize_idiff, idiff, 0, NULL, NULL));
-    printf("Set initial state.\n");
+    // printf("Set initial state.\n");
 
     /* Load and compile the kernel program(s) */
     char* source = read_kernel("kernel.cl");
@@ -285,14 +285,14 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
         abort();
     }
     ocl_check("clBuildProgram", err);
-    printf("Program created and built.\n");
+    // printf("Program created and built.\n");
 
     /* Create the kernels */
     kernel_cell = clCreateKernel(program, "cell_step", &err);
     ocl_check("clCreateKernel", err);
     kernel_diff = clCreateKernel(program, "diff_step", &err);
     ocl_check("clCreateKernel", err);
-    printf("Kernels created.\n");
+    // printf("Kernels created.\n");
 
     /* Pass arguments into kernels */
     ocl(clSetKernelArg(kernel_cell, 0, sizeof(n_cells), &n_cells));
@@ -304,7 +304,7 @@ int run(const int n_cells, const Real time_start, const Real time_end, const Rea
     ocl(clSetKernelArg(kernel_diff, 1, sizeof(g), &g));
     ocl(clSetKernelArg(kernel_diff, 2, sizeof(mbuf_state), &mbuf_state));
     ocl(clSetKernelArg(kernel_diff, 3, sizeof(mbuf_idiff), &mbuf_idiff));
-    printf("Arguments passed into kernels.\n");
+    // printf("Arguments passed into kernels.\n");
 
     /* Add log header */
     printf("engine.time");
@@ -342,7 +342,7 @@ for v in model.states():
 
 done:
     /* Tidy up */
-    printf("Tidying up.\n");
+    // printf("Tidying up.\n");
 
     ocl(clFlush(command_queue));
     ocl(clFinish(command_queue));
@@ -356,7 +356,7 @@ for v in model.states():
     free(state);
     free(idiff);
 
-    printf("Done.\n");
+    // printf("Done.\n");
 }
 
 int main()
diff --git a/myokit/formats/opencl/template/minilog.py b/myokit/formats/opencl/template/minilog.py
index 25eadf915..c82ad720c 100644
--- a/myokit/formats/opencl/template/minilog.py
+++ b/myokit/formats/opencl/template/minilog.py
@@ -69,7 +69,7 @@ def e(line, char, msg):
         # Delimiters
         quote = '"'
         delim = ','
-        with open(filename, 'rb') as f:
+        with open(filename, 'r', newline=None) as f:
             # Read header
             keys = []
             try:
diff --git a/myokit/formats/python/template/sim.py b/myokit/formats/python/template/sim.py
index ae0440aa8..4989aed5e 100644
--- a/myokit/formats/python/template/sim.py
+++ b/myokit/formats/python/template/sim.py
@@ -96,9 +96,9 @@ def _constants(self):
         Resets the state variables to their initial values
         """
 <?
-    for eq in model.inits():
-        if eq.lhs.var().parent() == comp:
-            print(tab2 + e(eq))
+    for var, value in zip(model.states(), model.initial_values(True)):
+        if var.parent() == comp:
+            print(tab2 + v(var, comp) + ' = ' + myokit.float.str(value))
     ?>    def update(self):
         """
         Re-calculates all values for the current time and state
diff --git a/myokit/formats/sbml/_api.py b/myokit/formats/sbml/_api.py
index 24c3177aa..d4f865f00 100644
--- a/myokit/formats/sbml/_api.py
+++ b/myokit/formats/sbml/_api.py
@@ -1298,7 +1298,7 @@ def set_rhs_sizes(sbml_model, variable_references, expression_references):
                         'with an initial value. Default is set to 1.')
 
                 # Promote size to state variable
-                var.promote(state_value=state_value)
+                var.promote(state_value)
 
             # Set RHS
             # (assignmentRule overwrites initialAssignment)
@@ -1383,7 +1383,7 @@ def set_rhs_reactants(
                         'value. Default is set to 0.')
 
                 # Promote size to state variable
-                var.promote(state_value=state_value)
+                var.promote(state_value)
                 var.set_rhs(myokit.Number(0))
 
             if (not var.rhs().is_literal()) or var.rhs().eval():
@@ -1473,7 +1473,7 @@ def set_rhs_products(
                         'value. Default is set to 0.')
 
                 # Promote size to state variable
-                var.promote(state_value=state_value)
+                var.promote(state_value)
                 var.set_rhs(myokit.Number(0))
 
             if (not var.rhs().is_literal()) or var.rhs().eval():
@@ -1546,7 +1546,7 @@ def set_rhs_species(
                         'value. Default is set to 0.')
 
                 # Promote size to state variable
-                var.promote(state_value=state_value)
+                var.promote(state_value)
 
             # Set RHS (reactions are dealt with elsewhere)
             expr = species.value()
@@ -1607,7 +1607,7 @@ def set_rhs_parameters(
                         ' initial value. Default is set to 1.')
 
                 # Promote size to state variable
-                var.promote(state_value=state_value)
+                var.promote(state_value)
 
             # Set RHS
             # (assignmentRule overwrites initialAssignment)
@@ -1673,7 +1673,7 @@ def set_rhs_stoichiometries(
                             ' 1.')
 
                     # Promote size to state variable
-                    var.promote(state_value=state_value)
+                    var.promote(state_value)
 
                 # Set RHS
                 # (assignmentRule overwrites initialAssignment)
diff --git a/myokit/lib/hh.py b/myokit/lib/hh.py
index e24dcd649..a620a8420 100644
--- a/myokit/lib/hh.py
+++ b/myokit/lib/hh.py
@@ -29,7 +29,8 @@ class HHModel(object):
     The model variables to treat as parameter are specified by the user when
     the model is created. Any other variables, for example state variables such
     as intercellular calcium or constants such as temperature, are fixed when
-    the current model is created and can no longer be changed.
+    the current model is created and can no longer be changed. Initial values
+    written as expressions are evaluated when the model is made.
 
     The current variable is optional.
 
@@ -83,7 +84,7 @@ class HHModel(object):
         # Load a model from disk
         model = myokit.load_model('some-model.mmt')
 
-        # Extract a markov model
+        # Extract a Hodgkin-Huxley style channel model
         mm = hh.HHModel.from_component(model.get('ina'))
 
     """
@@ -183,15 +184,18 @@ def __init__(self, model, states, parameters=None, current=None, vm=None):
                 'The membrane potential should not be the current variable.')
         del vm
 
+        # Get values of all states
+        # Note: Do this _before_ changing the model!
+        s = self._model.initial_values(True)
+        self._model.set_initial_values(s)  # Remove expressions
+        self._default_state = np.array(
+            [v.initial_value(True) for v in self._states])
+
         #
         # Demote unnecessary states and remove bindings
         #
-        # Get values of all states
-        # Note: Do this _before_ changing the model!
-        self._default_state = np.array([v.state_value() for v in self._states])
 
         # Freeze remaining, non-current-model states
-        s = self._model.state()   # Get state values before changing anything!
         # Note: list() cast is required so that we iterate over a static list,
         # otherwise we can get issues because the iterator depends on the model
         # (which we're changing).
@@ -346,19 +350,19 @@ def current(self):
 
     def default_membrane_potential(self):
         """
-        Returns this markov model's default membrane potential value.
+        Returns this HH model's default membrane potential value.
         """
         return self._default_inputs[0]
 
     def default_parameters(self):
         """
-        Returns this markov model's default parameter values
+        Returns this HH model's default parameter values
         """
         return list(self._default_inputs[1:])
 
     def default_state(self):
         """
-        Returns this markov model's default state values.
+        Returns this HH model's default state values.
         """
         return list(self._default_state)
 
diff --git a/myokit/lib/markov.py b/myokit/lib/markov.py
index 150ba697a..a51379bd0 100644
--- a/myokit/lib/markov.py
+++ b/myokit/lib/markov.py
@@ -37,7 +37,8 @@ class LinearModel(object):
     The model variables to treat as parameter are specified by the user when
     the model is created. Any other variables, for example state variables such
     as intercellular calcium or constants such as temperature, are fixed when
-    the markov model is created and can no longer be changed.
+    the markov model is created and can no longer be changed. Initial values
+    written as expressions are evaluated when the model is made.
 
     To create a :class:`Markov`, pass in a :class:`myokit.Model` and select a
     list of states. All other states will be fixed at their current value and
@@ -205,15 +206,17 @@ def __init__(self, model, states, parameters=None, current=None, vm=None):
                 'The membrane potential should not be the current variable.')
         del vm
 
+        # Get values of all states
+        # Note: Do this _before_ changing the model!
+        s = self._model.initial_values(True)
+        self._default_state = np.array(
+            [v.initial_value(True) for v in self._states])
+
         #
         # Demote unnecessary states, remove bindings and validate model.
         #
-        # Get values of all states
-        # Note: Do this _before_ changing the model!
-        self._default_state = np.array([v.state_value() for v in self._states])
 
         # Freeze remaining, non-markov-model states
-        s = self._model.state()   # Get state values before changing anything!
         for k, state in enumerate(self._model.states()):
             if state not in self._states:
                 state.demote()
diff --git a/myokit/lib/plots.py b/myokit/lib/plots.py
index 79d80299a..c0c01916f 100644
--- a/myokit/lib/plots.py
+++ b/myokit/lib/plots.py
@@ -350,7 +350,10 @@ def cumulative_current(
             colors.extend(colors)
     else:
         # Colormap
-        cmap = matplotlib.cm.get_cmap(name='tab20')
+        try:
+            cmap = matplotlib.colormaps['tab20']
+        except AttributeError:  # pragma: no cover
+            cmap = matplotlib.cm.get_cmap(name='tab20')
         colors = [cmap(i) for i in range(nc)]
 
     # Line drawing keyword arguments


Test patch:
diff --git a/myokit/formats/opencl/template/test b/myokit/formats/opencl/template/test.sh
similarity index 100%
rename from myokit/formats/opencl/template/test
rename to myokit/formats/opencl/template/test.sh
diff --git a/myokit/tests/data/beeler-1977-model-compare-b.mmt b/myokit/tests/data/beeler-1977-model-compare-b.mmt
index be7bf6618..567886408 100644
--- a/myokit/tests/data/beeler-1977-model-compare-b.mmt
+++ b/myokit/tests/data/beeler-1977-model-compare-b.mmt
@@ -2,7 +2,7 @@
 name: beeler-1977-with-differences
 desc: """
     The 1997 Beeler Reuter model of the AP in ventricular myocytes.
-    
+
     This test version has some differences
 
     Reference:
@@ -18,7 +18,7 @@ membrane.V  = -84.622
 calcium.Cai = 2e-7
 ina.m       = 0.01
 ina.h       = 1.0
-ina.j       = 0.98
+ina.j       = 0.97 + 0.01
 isiz.d       = 0.003
 isiz.f       = 1.0
 
diff --git a/myokit/tests/data/clancy-1999-fitting.mmt b/myokit/tests/data/clancy-1999-fitting.mmt
index e8d684077..b67ada3a9 100644
--- a/myokit/tests/data/clancy-1999-fitting.mmt
+++ b/myokit/tests/data/clancy-1999-fitting.mmt
@@ -143,6 +143,7 @@ dot(O)  = -O  * (a2  + b13     ) + C1 * a13 + IF * b2
 dot(IF) = -IF * (a4  + b2  + a3) + IS * b4  + O  * a2  + C1 * b3
 dot(IS) = -IS * (b4            ) + IF * a4
 i = gmax * O * (V - E)
+p = 1e-4
 
 [[protocol]]
 # Level  Start    Length   Period   Multiplier
diff --git a/myokit/tests/test_aux.py b/myokit/tests/test_aux.py
index c7ef5ae8a..c6079a7fc 100755
--- a/myokit/tests/test_aux.py
+++ b/myokit/tests/test_aux.py
@@ -163,6 +163,7 @@ def test_model_comparison(self):
             '[1] Missing User function <g(1)>.',
             '[x] Mismatched Time variable: [1]<engine.time> [2]<engine.toim>',
             '[x] Mismatched Initial value for <ina.h>',
+            '[x] Mismatched Initial value for <ina.j>',
             '[x] Mismatched State at position 5: [1]<isi.d> [2]<isiz.d>',
             '[x] Mismatched State at position 6: [1]<isi.f> [2]<isiz.f>',
             '[2] Missing state at position 7',
@@ -193,9 +194,6 @@ def test_model_comparison(self):
             '  ' + str(len(differences)) + ' differences found',
         ]
 
-        # Show massive diff messages
-        self.maxDiff = None
-
         caught_differences = set(capture.text().splitlines())
         live = set(live)
         self.assertEqual(live, caught_differences)
@@ -276,8 +274,6 @@ def test_run(self):
             'import myokit',
             'm = get_model()',  # Test magic methods
             'p = get_protocol()',
-            's = myokit.Simulation(m, p)',
-            's.run(200)',
         ])
         with myokit.tools.capture():
             myokit.run(m, p, x)
@@ -513,6 +509,9 @@ def almost_equal(x, y):
         # Test positive/negative zero comparison
         m1 = myokit.Model()
         c = m1.add_component('c')
+        t = c.add_variable('t')
+        t.set_binding('time')
+        t.set_rhs(0)
         x = c.add_variable('x')
         x.promote(1)
         x.set_rhs('-0.0')
diff --git a/myokit/tests/test_cellml_v1_api.py b/myokit/tests/test_cellml_v1_api.py
index 048aa7bbe..7cc0621ff 100755
--- a/myokit/tests/test_cellml_v1_api.py
+++ b/myokit/tests/test_cellml_v1_api.py
@@ -730,8 +730,8 @@ def test_evaluating_states(self):
         # Recreate myokit model and test states
         mm = cm.myokit_model()
         mm.validate()
-        state_1 = m.state()
-        state_2 = mm.state()
+        state_1 = m.initial_values(True)
+        state_2 = mm.initial_values(True)
         states_1 = [x.name() for x in m.states()]
         states_2 = [x.name() for x in mm.states()]
         state_2 = [state_2[states_2.index(x)] for x in states_1]
@@ -830,7 +830,7 @@ def test_c2m_basic(self):
 
         # Check state
         self.assertTrue(mx.is_state())
-        self.assertEqual(mx.state_value(), 0.123)
+        self.assertEqual(mx.initial_value(True), 0.123)
 
         # Check binding
         self.assertEqual(mt.binding(), 'time')
diff --git a/myokit/tests/test_cellml_v2_api.py b/myokit/tests/test_cellml_v2_api.py
index c464b2de0..c992196f1 100755
--- a/myokit/tests/test_cellml_v2_api.py
+++ b/myokit/tests/test_cellml_v2_api.py
@@ -909,8 +909,8 @@ def test_evaluating_states(self):
         # Recreate myokit model and test states
         mm = cm.myokit_model()
         mm.validate()
-        state_1 = m.state()
-        state_2 = mm.state()
+        state_1 = m.initial_values(True)
+        state_2 = mm.initial_values(True)
         states_1 = [x.name() for x in m.states()]
         states_2 = [x.name() for x in mm.states()]
         state_2 = [state_2[states_2.index(x)] for x in states_1]
@@ -1013,7 +1013,7 @@ def test_c2m_basic(self):
 
         # Check state
         self.assertTrue(mx.is_state())
-        self.assertEqual(mx.state_value(), 0.123)
+        self.assertEqual(mx.initial_value(True), 0.123)
 
         # Check binding
         self.assertEqual(mt.binding(), 'time')
diff --git a/myokit/tests/test_component.py b/myokit/tests/test_component.py
index aa972e057..cffc40fdc 100755
--- a/myokit/tests/test_component.py
+++ b/myokit/tests/test_component.py
@@ -58,7 +58,7 @@ def test_move_variable(self):
                 myokit.Name(y)
             )
         ))
-        x.set_state_value(10)
+        x.set_initial_value(10)
         y.set_rhs(myokit.Plus(
             myokit.Multiply(
                 myokit.PrefixMinus(myokit.Name(c)), myokit.Name(y)
@@ -67,7 +67,7 @@ def test_move_variable(self):
                 myokit.Multiply(myokit.Name(d), myokit.Name(x)), myokit.Name(y)
             )
         ))
-        y.set_state_value(5)
+        y.set_initial_value(5)
         Z = m.add_component('Z')
         t = Z.add_variable('total')
         t.set_rhs(myokit.Plus(myokit.Name(x), myokit.Name(y)))
@@ -117,92 +117,92 @@ def test_remove_variable(self):
         m = myokit.Model('LotkaVolterra')
 
         # Add a variable 'a'
-        X = m.add_component('X')
+        z = m.add_component('z')
 
         # Simplest case
-        a = X.add_variable('a')
-        self.assertEqual(X.count_variables(), 1)
-        X.remove_variable(a)
-        self.assertEqual(X.count_variables(), 0)
-        self.assertRaises(Exception, X.remove_variable, a)
+        a = z.add_variable('a')
+        self.assertEqual(z.count_variables(), 1)
+        z.remove_variable(a)
+        self.assertEqual(z.count_variables(), 0)
+        self.assertRaises(Exception, z.remove_variable, a)
 
         # Test re-adding
-        a = X.add_variable('a')
+        a = z.add_variable('a')
         a.set_rhs(myokit.Number(5))
-        self.assertEqual(X.count_variables(), 1)
+        self.assertEqual(z.count_variables(), 1)
 
         # Test deleting dependent variables
-        b = X.add_variable('b')
-        self.assertEqual(X.count_variables(), 2)
+        b = z.add_variable('b')
+        self.assertEqual(z.count_variables(), 2)
         b.set_rhs(myokit.Plus(myokit.Number(3), myokit.Name(a)))
 
         # Test blocking of removal
-        self.assertRaises(myokit.IntegrityError, X.remove_variable, a)
-        self.assertEqual(X.count_variables(), 2)
+        self.assertRaises(myokit.IntegrityError, z.remove_variable, a)
+        self.assertEqual(z.count_variables(), 2)
 
         # Test removal in the right order
-        X.remove_variable(b)
-        self.assertEqual(X.count_variables(), 1)
-        X.remove_variable(a)
-        self.assertEqual(X.count_variables(), 0)
+        z.remove_variable(b)
+        self.assertEqual(z.count_variables(), 1)
+        z.remove_variable(a)
+        self.assertEqual(z.count_variables(), 0)
 
         # Test reference to current state variable values
-        a = X.add_variable('a')
+        a = z.add_variable('a')
         a.set_rhs(myokit.Number(5))
         a.promote()
-        b = X.add_variable('b')
+        b = z.add_variable('b')
         b.set_rhs(myokit.Plus(myokit.Number(3), myokit.Name(a)))
-        self.assertRaises(myokit.IntegrityError, X.remove_variable, a)
-        X.remove_variable(b)
-        X.remove_variable(a)
-        self.assertEqual(X.count_variables(), 0)
+        self.assertRaises(myokit.IntegrityError, z.remove_variable, a)
+        z.remove_variable(b)
+        z.remove_variable(a)
+        self.assertEqual(z.count_variables(), 0)
 
         # Test reference to current state variable values with "self"-ref
-        a = X.add_variable('a')
+        a = z.add_variable('a')
         a.promote()
         a.set_rhs(myokit.Name(a))
-        X.remove_variable(a)
+        z.remove_variable(a)
 
         # Test it doesn't interfere with normal workings
-        a = X.add_variable('a')
+        a = z.add_variable('a')
         a.promote()
         a.set_rhs(myokit.Name(a))
-        b = X.add_variable('b')
+        b = z.add_variable('b')
         b.set_rhs(myokit.Name(a))
-        self.assertRaises(myokit.IntegrityError, X.remove_variable, a)
-        X.remove_variable(b)
-        X.remove_variable(a)
+        self.assertRaises(myokit.IntegrityError, z.remove_variable, a)
+        z.remove_variable(b)
+        z.remove_variable(a)
 
         # Test reference to dot
-        a = X.add_variable('a')
+        a = z.add_variable('a')
         a.set_rhs(myokit.Number(5))
         a.promote()
-        b = X.add_variable('b')
+        b = z.add_variable('b')
         b.set_rhs(myokit.Derivative(myokit.Name(a)))
-        self.assertRaises(myokit.IntegrityError, X.remove_variable, a)
-        X.remove_variable(b)
-        X.remove_variable(a)
+        self.assertRaises(myokit.IntegrityError, z.remove_variable, a)
+        z.remove_variable(b)
+        z.remove_variable(a)
 
         # Test if orphaned
         self.assertIsNone(b.parent())
 
         # Test deleting variable with nested variables
-        a = X.add_variable('a')
+        a = z.add_variable('a')
         b = a.add_variable('b')
         b.set_rhs(myokit.Plus(myokit.Number(2), myokit.Number(2)))
         a.set_rhs(myokit.Multiply(myokit.Number(3), myokit.Name(b)))
-        self.assertRaises(myokit.IntegrityError, X.remove_variable, a)
+        self.assertRaises(myokit.IntegrityError, z.remove_variable, a)
         self.assertEqual(a.count_variables(), 1)
-        self.assertEqual(X.count_variables(), 1)
+        self.assertEqual(z.count_variables(), 1)
 
         # Test recursive deleting
-        X.remove_variable(a, recursive=True)
+        z.remove_variable(a, recursive=True)
         self.assertEqual(a.count_variables(), 0)
-        self.assertEqual(X.count_variables(), 0)
+        self.assertEqual(z.count_variables(), 0)
 
         # Test deleting variable with nested variables that depend on each
         # other
-        a = X.add_variable('a')
+        a = z.add_variable('a')
         b = a.add_variable('b')
         c = a.add_variable('c')
         d = a.add_variable('d')
@@ -211,12 +211,12 @@ def test_remove_variable(self):
         b.set_rhs('2 * a - d')
         c.set_rhs('a + b + d')
         d.set_rhs('3 * a')
-        self.assertRaises(myokit.IntegrityError, X.remove_variable, a)
+        self.assertRaises(myokit.IntegrityError, z.remove_variable, a)
         self.assertEqual(a.count_variables(), 3)
-        self.assertEqual(X.count_variables(), 1)
-        X.remove_variable(a, recursive=True)
+        self.assertEqual(z.count_variables(), 1)
+        z.remove_variable(a, recursive=True)
         self.assertEqual(a.count_variables(), 0)
-        self.assertEqual(X.count_variables(), 0)
+        self.assertEqual(z.count_variables(), 0)
 
         # Test if removed from model's label and binding lists
         m = myokit.Model()
@@ -236,6 +236,68 @@ def test_remove_variable(self):
         self.assertIs(m.binding('time'), None)
         self.assertIs(m.label('membrane_potential'), None)
 
+        # Test variable can't be removed if used in initial expression
+        m = myokit.Model()
+        c = m.add_component('c')
+        p = c.add_variable('p')
+        p.set_rhs(2)
+        x = c.add_variable('x')
+        x.set_rhs(1)
+        x.promote('1 + c.p')
+        self.assertRaises(myokit.IntegrityError, c.remove_variable, p)
+        q = c.add_variable('q')
+        q.set_rhs(10)
+        x.set_initial_value('1 / c.q')
+        c.remove_variable(p)
+        self.assertRaises(myokit.IntegrityError, c.remove_variable, q)
+        x.set_initial_value('sqrt(2)')
+        c.remove_variable(q)
+
+        # But it's OK if we're deleting a whole component
+        m = myokit.Model()
+        c = m.add_component('c')
+        p = c.add_variable('p')
+        p.set_rhs(2)
+        x = c.add_variable('x')
+        x.set_rhs(1)
+        x.promote('1 + c.p')
+        m.remove_component('c')
+
+    def test_remove_variable_with_alias(self):
+        # Test cloning of a variable with an alias after an add / remove event.
+
+        m = myokit.Model('AddRemoveClone')
+        c = m.add_component('c')
+        p = c.add_variable('p')
+        p.set_binding('time')
+        p.set_rhs(0)
+        q = c.add_variable('q')
+        q.set_rhs(12)
+        m.validate()    # Raises error if not ok
+        m.clone()       # Raises error if not ok
+        d = m.add_component('d')
+        d.add_alias('bert', p)
+        e = d.add_variable('e')
+        e.set_rhs('10 * bert')
+        m.validate()
+        m.clone()
+        d.add_alias('ernie', q)
+        m.validate()
+        m.clone()
+        c.remove_variable(q)
+        m.validate()
+        m.clone()   # Will raise error if alias isn't deleted
+
+    def test_resolve(self):
+        # Test if an error is raised when a variable can't be resolved.
+
+        m = myokit.Model('Resolve')
+        c = m.add_component('c')
+        p = c.add_variable('p')
+        q = c.add_variable('q')
+        p.set_rhs('10 * q')
+        self.assertRaises(myokit.ParseError, q.set_rhs, '10 * r')
+
     def test_sequence_interface(self):
         # Test the sequence interface implementation
 
diff --git a/myokit/tests/test_config.py b/myokit/tests/test_config.py
index 63e37cd17..892bde666 100755
--- a/myokit/tests/test_config.py
+++ b/myokit/tests/test_config.py
@@ -17,7 +17,7 @@
 from myokit.tests import TemporaryDirectory, WarningCollector
 
 
-debug = True
+debug = False
 
 
 # Simple test
diff --git a/myokit/tests/test_io.py b/myokit/tests/test_io.py
index ede91b7c6..0579c2a94 100755
--- a/myokit/tests/test_io.py
+++ b/myokit/tests/test_io.py
@@ -340,22 +340,22 @@ def test_load_save_state(self):
         with TemporaryDirectory() as d:
             # Test save and load without model
             f = d.path('state.txt')
-            myokit.save_state(f, m.state())
-            self.assertEqual(myokit.load_state(f), m.state())
+            myokit.save_state(f, m.initial_values(True))
+            self.assertEqual(myokit.load_state(f), m.initial_values(True))
 
             # Test save and load with model argument
-            myokit.save_state(f, m.state(), m)
-            self.assertEqual(myokit.load_state(f, m), m.state())
+            myokit.save_state(f, m.initial_values(True), m)
+            self.assertEqual(myokit.load_state(f, m), m.initial_values(True))
 
             # Save without, load with model
-            myokit.save_state(f, m.state())
-            self.assertEqual(myokit.load_state(f, m), m.state())
+            myokit.save_state(f, m.initial_values(True))
+            self.assertEqual(myokit.load_state(f, m), m.initial_values(True))
 
             # Save with model, load without
             # Loaded version is dict!
-            myokit.save_state(f, m.state(), m)
-            dct = dict(zip([v.qname() for v in m.states()], m.state()))
-            self.assertEqual(myokit.load_state(f), dct)
+            myokit.save_state(f, m.initial_values(True), m)
+            self.assertEqual(myokit.load_state(f), dict(zip(
+                [v.qname() for v in m.states()], m.initial_values(True))))
 
     def test_load_save_state_bin(self):
         # Test loading/saving state in binary format.
@@ -364,13 +364,15 @@ def test_load_save_state_bin(self):
 
             # Test save and load with double precision
             f = d.path('state.bin')
-            myokit.save_state_bin(f, m.state())
-            self.assertEqual(myokit.load_state_bin(f), m.state())
+            myokit.save_state_bin(f, m.initial_values(True))
+            self.assertEqual(myokit.load_state_bin(f), m.initial_values(True))
 
             # Test save and load with single precision
             f = d.path('state.bin')
-            myokit.save_state_bin(f, m.state(), myokit.SINGLE_PRECISION)
-            d = np.array(myokit.load_state_bin(f)) - np.array(m.state())
+            myokit.save_state_bin(
+                f, m.initial_values(True), myokit.SINGLE_PRECISION)
+            d = np.array(myokit.load_state_bin(f))
+            d -= np.array(m.initial_values(True))
             self.assertTrue(np.all(np.abs(d) < 1e-5))   # Not very precise!
 
 
diff --git a/myokit/tests/test_jacobian_calculator.py b/myokit/tests/test_jacobian_calculator.py
index 4ebb7c1a1..571d64194 100755
--- a/myokit/tests/test_jacobian_calculator.py
+++ b/myokit/tests/test_jacobian_calculator.py
@@ -37,16 +37,16 @@ def test_simple(self):
 
         # Test if still runs with initial x all zero
         # (But does cause a linear algebra issue in this case)
-        x = np.zeros(len(m.state()))
+        x = np.zeros(m.count_states())
         c.newton_root(x, damping=0.01, max_iter=50)
 
         # Test if still works with a single zero (Enno's bug)
-        x = np.array(m.state())
+        x = np.array(m.initial_values(True))
         x[1] = 0
         x, f, j, e = c.newton_root(x, damping=0.01, max_iter=50)
 
         # Test quick return
-        x = np.array(m.state())
+        x = np.array(m.initial_values(True))
         x[0] = 0
         x2, f, j, e = c.newton_root(damping=0.01, max_iter=1)
         self.assertTrue(np.sum((x2 - x)**2) > 10)
@@ -58,12 +58,12 @@ def test_simple(self):
             ValueError, 'Damping', c.newton_root, damping=1.1)
 
         # Missing a state
-        x = m.state()[:-1]
+        x = m.initial_values(True)[:-1]
         self.assertRaisesRegex(
             ValueError, 'must have length', c.calculate, x)
 
         # Non-numbers in state
-        x = m.state()
+        x = m.initial_values(True)
         x[0] = 'Hello'
         self.assertRaisesRegex(
             ValueError, 'floats', c.calculate, x)
diff --git a/myokit/tests/test_lib_hh.py b/myokit/tests/test_lib_hh.py
index 99abd9e5a..7f8ed25ea 100755
--- a/myokit/tests/test_lib_hh.py
+++ b/myokit/tests/test_lib_hh.py
@@ -365,13 +365,13 @@ def test_rush_larsen_conversion(self):
         a = m2.get('ikr.a')
         rl = hh.get_rl_expression(a, myokit.Name(dt))
         a1 = rl.eval()
-        a2 = a.state_value() + dt.eval() * a.rhs().eval()
+        a2 = a.initial_value(True) + dt.eval() * a.rhs().eval()
         self.assertAlmostEqual(a1, a2)
         # And for r
         r = m2.get('ikr.r')
         rl = hh.get_rl_expression(r, myokit.Name(dt))
         r1 = rl.eval()
-        r2 = r.state_value() + dt.eval() * r.rhs().eval()
+        r2 = r.initial_value(True) + dt.eval() * r.rhs().eval()
         self.assertAlmostEqual(r1, r2)
 
         # Dt must be an expression
@@ -608,7 +608,7 @@ def test_hh_model_steady_state(self):
         # Test if derivatives are zero
         for k, x in enumerate(['ina.m', 'ina.h', 'ina.j']):
             x = model.get(x)
-            x.set_state_value(ss[k])
+            x.set_initial_value(ss[k])
             self.assertAlmostEqual(x.eval(), 0)
 
         # Test arguments
@@ -645,6 +645,20 @@ def test_automatic_creation_with_v_independence(self):
 
         self.assertEqual(len(m.states()), 3)
 
+    def test_initial_value_conversion(self):
+        # Tests that initial value expressions are converted to floats
+
+        model = myokit.parse_model(MODEL)
+        model.get('ikr.a').set_initial_value('1 / sqrt(7)')
+        model.get('ikr.r').set_initial_value('log(binding.koff, 10)')
+        m = hh.HHModel.from_component(model.get('ikr'))
+        x0 = m.default_state()
+        self.assertEqual(len(x0), 2)
+        self.assertIsInstance(x0[0], float)
+        self.assertIsInstance(x0[1], float)
+        self.assertAlmostEqual(x0[0], 0.377964473)
+        self.assertAlmostEqual(x0[1], -5)
+
 
 class AnalyticalSimulationTest(unittest.TestCase):
     """
diff --git a/myokit/tests/test_lib_markov.py b/myokit/tests/test_lib_markov.py
index ac4c3dff3..362268286 100755
--- a/myokit/tests/test_lib_markov.py
+++ b/myokit/tests/test_lib_markov.py
@@ -155,8 +155,8 @@ def test_manual_creation(self):
 
         # States must sum to 1
         m2 = model.clone()
-        m2.get(states[0]).set_state_value(0.6)
-        m2.get(states[1]).set_state_value(0.6)
+        m2.get(states[0]).set_initial_value(0.6)
+        m2.get(states[1]).set_initial_value(0.6)
         self.assertRaisesRegex(
             markov.LinearModelError, 'sum of states',
             markov.LinearModel, m2, states, parameters, current)
@@ -283,7 +283,7 @@ def test_linear_model_steady_state_1(self):
 
         # Check that derivatives with ss are close to zero
         ss = list(ss)
-        model.set_state(ss + ss)    # Model has 2 ina's
+        model.set_initial_values(ss + ss)    # Model has 2 ina's
         derivs = model.evaluate_derivatives()
         for i in range(len(ss)):
             self.assertAlmostEqual(0, derivs[i])
@@ -310,7 +310,7 @@ def test_linear_model_steady_state_2(self):
         self.assertTrue(np.all(ss <= 1))
 
         # Check that derivatives with ss are close to zero
-        model.set_state(ss)
+        model.set_initial_values(ss)
         derivs = model.evaluate_derivatives()
         for i in range(len(ss)):
             self.assertAlmostEqual(0, derivs[i])
@@ -331,6 +331,22 @@ def test_rates(self):
             ValueError, 'Illegal parameter vector size',
             m.rates, parameters=[0.01] * 22)
 
+    def test_initial_value_conversion(self):
+        # Tests that initial value expressions are converted to floats
+
+        fname = os.path.join(DIR_DATA, 'clancy-1999-fitting.mmt')
+        model = myokit.load_model(fname)
+        model.get('ina.C3').set_initial_value('1 / sqrt(7)')
+        model.get('ina.C2').set_initial_value('-1 / log(ina_ref.p)')
+        model.get('ina.C1').set_initial_value(0.5134619065149598)
+        m = markov.LinearModel.from_component(model.get('ina'))
+        x0 = m.default_state()
+        self.assertEqual(len(x0), 6)
+        self.assertIsInstance(x0[0], float)
+        self.assertIsInstance(x0[1], float)
+        self.assertAlmostEqual(x0[0], 0.3779644730092272)
+        self.assertAlmostEqual(x0[1], 0.10857362047581297)
+
 
 class AnalyticalSimulationTest(unittest.TestCase):
     """
diff --git a/myokit/tests/test_model.py b/myokit/tests/test_model.py
index 091c18bec..1fad446dd 100755
--- a/myokit/tests/test_model.py
+++ b/myokit/tests/test_model.py
@@ -301,6 +301,26 @@ def test_clone(self):
         self.assertTrue(m1.has_parse_info())
         self.assertFalse(m2.has_parse_info())
 
+        # Test initial value expressions are maintained
+        m1 = myokit.Model('moo')
+        component = m1.add_component('c')
+        p = component.add_variable('p')
+        q = component.add_variable('q')
+        r = component.add_variable('r')
+        s = component.add_variable('s')
+        z = component.add_variable('z')
+        z.set_rhs(123)
+        p.promote(1)
+        q.promote('1 + exp(3)')
+        r.promote('c.z')
+        s.promote('1 [g] + c.z / 2 [1/g]')
+        m2 = m1.clone()
+        self.assertEqual(m2.get('c.p').initial_value(), myokit.Number(1))
+        self.assertEqual(m2.get('c.q').initial_value().code(), '1 + exp(3)')
+        self.assertEqual(m2.get('c.r').initial_value(), m2.get('c.z').lhs())
+        self.assertEqual(m2.get('c.s').initial_value().code(),
+                         '1 [g] + c.z / 2 [1/g]')
+
     def test_code(self):
         # Test :meth:`Model.code()`.
 
@@ -355,44 +375,41 @@ def test_code(self):
             '13 d = comp1.a\n'
         )
 
-    def test_is_similar(self):
-        # Check that equality takes both code() and unames into account
-
-        # Test without custom reserved names
-        m1 = myokit.load_model('example')
-        m2 = m1.clone()
-        self.assertIsInstance(m2, myokit.Model)
-        self.assertFalse(m1 is m2)
-        self.assertNotEqual(m1, m2)
-        self.assertNotEqual(m2, m1)
-        self.assertTrue(m1.is_similar(m2, False))
-        self.assertTrue(m1.is_similar(m2, True))
-        self.assertTrue(m2.is_similar(m1, False))
-        self.assertTrue(m2.is_similar(m1, True))
-        self.assertTrue(m1.is_similar(m1, True))
-        self.assertTrue(m2.is_similar(m2, False))
-
-        # Test with none-model
-        self.assertFalse(m1.is_similar(None))
-        self.assertFalse(m1.is_similar(m1.code()))
-
-        # Add reserved names
-        m1.reserve_unique_names('bertie')
-        self.assertFalse(m1.is_similar(m2))
-        m1.reserve_unique_names('clair')
-        self.assertFalse(m1.is_similar(m2))
-        m2.reserve_unique_names('clair', 'bertie')
-        self.assertTrue(m1.is_similar(m2))
-
-        # Add reserved name prefixes
-        m1.reserve_unique_name_prefix('aa', 'bb')
-        m1.reserve_unique_name_prefix('cc', 'dd')
-        self.assertFalse(m1.is_similar(m2))
-        m2.reserve_unique_name_prefix('aa', 'bb')
-        m2.reserve_unique_name_prefix('cc', 'ee')
-        self.assertFalse(m1.is_similar(m2))
-        m2.reserve_unique_name_prefix('cc', 'dd')
-        self.assertTrue(m1.is_similar(m2))
+        # Test initial value expressions are maintained
+        model = myokit.Model('moo')
+        component = model.add_component('c')
+        p = component.add_variable('p')
+        q = component.add_variable('q')
+        r = component.add_variable('r')
+        s = component.add_variable('s')
+        z = component.add_variable('z')
+        z.set_rhs(123)
+        p.promote(1)
+        q.promote('1 + exp(3)')
+        r.promote('c.z')
+        s.promote('1 [g] + c.z / 2 [1/g]')
+
+        p.set_rhs(1)
+        q.set_rhs(2)
+        r.set_rhs(3)
+        s.set_rhs('4 [g/s]')
+
+        self.assertEqual(model.code(line_numbers=False), '\n'.join([
+            '[[model]]',
+            'name: moo',
+            '# Initial values',
+            'c.p = 1',
+            'c.q = 1 + exp(3)',
+            'c.r = c.z',
+            'c.s = 1 [g] + c.z / 2 [1/g]',
+            '',
+            '[c]',
+            'dot(p) = 1',
+            'dot(q) = 2',
+            'dot(r) = 3',
+            'dot(s) = 4 [g/s]',
+            'z = 123\n\n',
+        ]))
 
     def test_evaluate_derivatives(self):
         # Test Model.evaluate_derivatives().
@@ -494,7 +511,7 @@ def test_format_state(self):
 
         # Test with invalid state argument
         self.assertRaisesRegex(
-            ValueError, r'list of \(8\)', m.format_state, [1, 2, 3])
+            ValueError, r'sequence of \(8\)', m.format_state, [1, 2, 3])
 
         # Test with precision argument
         state1 = [1, 2, 3, 4, 5, 6, 7, 8]
@@ -526,9 +543,26 @@ def test_format_state(self):
 
         # Test with invalid second state argument
         self.assertRaisesRegex(
-            ValueError, r'list of \(8\)', m.format_state,
+            ValueError, r'sequence of \(8\)', m.format_state,
             [1, 2, 3, 4, 5, 6, 7, 8], [1, 2, 3])
 
+        # Test with model that has expressions in its initial state.
+        # These will be evaluated (the method handles states, not initial state
+        # expressions).
+        m.get('ina.m').set_initial_value('sqrt(4) / 4')
+        m.get('ina.h').set_initial_value('ina.gNa / 20')
+        self.assertEqual(
+            m.format_state(),
+            'membrane.V = -84.5286\n'
+            'ina.m      = 0.5\n'
+            'ina.h      = 0.8\n'
+            'ina.j      = 0.995484\n'
+            'ica.d      = 3e-06\n'
+            'ica.f      = 1.0\n'
+            'ik.x       = 0.0057\n'
+            'ica.Ca_i   = 0.0002'
+        )
+
     def test_format_state_derivatives(self):
         # Test Model.format_state_derivatives().
 
@@ -585,7 +619,7 @@ def almost_equal(x, y):
 
         # Test with invalid state argument
         self.assertRaisesRegex(
-            ValueError, r'list of \(8\)',
+            ValueError, r'sequence of \(8\)',
             m.format_state_derivatives, [1, 2, 3])
 
         # Test with state and precision argument
@@ -636,9 +670,28 @@ def almost_equal(x, y):
 
         # Test with invalid derivs argument
         self.assertRaisesRegex(
-            ValueError, r'list of \(8\)', m.format_state_derivatives,
+            ValueError, r'sequence of \(8\)', m.format_state_derivatives,
             [1, 2, 3, 4, 5, 6, 7, 8], [1, 2, 3])
 
+        # Test without expressions in the initial state (these will be
+        # evaluated)
+        m.get('ina.m').set_initial_value('sqrt(0.25)')
+        m.get('ina.h').set_initial_value('ina.gNa / 20')
+        x = m.format_state_derivatives().splitlines()
+        y = [
+'membrane.V = -84.5286                   dot =  2.46843981754470434e+02', # noqa
+'ina.m      = 0.5                        dot = -8.68225658924664856e+01', # noqa
+'ina.h      = 0.8                        dot =  4.89677127030544446e-02', # noqa
+'ina.j      = 0.995484                   dot = -3.70409866928434243e-04', # noqa
+'ica.d      = 3e-06                      dot =  3.68067721821794798e-04', # noqa
+'ica.f      = 1.0                        dot = -3.55010150519739432e-07', # noqa
+'ik.x       = 0.0057                     dot = -2.04613933160084307e-07', # noqa
+'ica.Ca_i   = 0.0002                     dot = -6.99430692442154227e-06'  # noqa
+        ]
+        for a, b in zip(x, y):
+            self.assertTrue(almost_equal(a, b))
+        self.assertEqual(len(x), len(y))
+
     def test_get(self):
         # Test Model.get().
 
@@ -725,10 +778,12 @@ def test_import_component(self):
         # Source model, to import stuff from
         ms = myokit.parse_model('''
             [[model]]
+            n.c = 1 + n.p
+            n.b = sqrt(0.25)
             p.b = 0.2
             q.e = 0.2
-            x.a = 0.2
             y.e = 0.2
+            x.a = 0.2
 
             [e]
             t = 0 [s] bind time
@@ -778,14 +833,20 @@ def test_import_component(self):
                 sub_e = 2 * e
                     in [m]
 
-            # another group component but this isn't independant
+            # Another group component but this isn't independent
             [z]
             use x.d
             use y.e
             use e.h
-
             f = d * e * h
                 in [m/s]
+
+            # A component with expressions in initial states
+            [n]
+            a = 1
+            dot(b) = 1 [1/s]
+            dot(c) = 2 [1/s]
+            p = 0.1
         ''')
         ms.validate()
         ms.check_units(myokit.UNIT_STRICT)
@@ -916,6 +977,19 @@ def test_import_component(self):
         self.assertEqual(m1['x'].code(), ms['x'].code())
         self.assertEqual(m1['y'].code(), ms['y'].code())
         self.assertTrue(ms.is_similar(ms_unaltered, True))
+        # Check that state order is preserved
+        self.assertLess(m1.get('y.e').indice(), m1.get('x.a').indice())
+
+        # Import multiple components and rename
+        component_list = [ms['x'], ms['y']]
+        m1.import_component(component_list, new_name=['xx', 'yy'])
+        self.assertTrue(m1.has_component('xx'))
+        self.assertTrue(m1.has_component('yy'))
+        self.assertFalse(m1['xx'] is ms['x'])
+        self.assertFalse(m1['yy'] is ms['y'])
+        self.assertTrue(ms.is_similar(ms_unaltered, True))
+        # Check that state order is preserved
+        self.assertLess(m1.get('yy.e').indice(), m1.get('xx.a').indice())
 
         # Import 1 component in list
         m1.import_component([ms['p']], new_name='p3')
@@ -927,6 +1001,26 @@ def test_import_component(self):
         self.assertEqual(cs, c1)
         self.assertTrue(ms.is_similar(ms_unaltered, True))
 
+        # Import state variables with expressions in their initial values
+        m1_unaltered = m1.clone()
+        m1.import_component(ms['n'])
+        self.assertEqual(m1.get('n.b').initial_value().code(), 'sqrt(0.25)')
+        self.assertEqual(
+            m1.get('n.c').initial_value(),
+            myokit.Plus(myokit.Number(1), myokit.Name(m1.get('n.p'))))
+        # Check that state order is preserved
+        self.assertLess(m1.get('n.c').indice(), m1.get('n.b').indice())
+        m1 = m1_unaltered
+
+        # ...including references to parameters in other components
+        ms2 = ms.clone()
+        ms2.get('n.c').set_initial_value('1 + e.h')
+        ms2_unaltered = ms2.clone()
+        self.assertRaises(
+            myokit.VariableMappingError, m1.import_component, ms2['n'])
+        self.assertTrue(m1.is_similar(m1_unaltered, True))
+        self.assertTrue(ms2.is_similar(ms2_unaltered, True))
+
         # Try and fail to import r without a mapping
         m1_unaltered = m1.clone()
         self.assertRaises(
@@ -1154,7 +1248,7 @@ def test_import_component_units(self):
         self.assertTrue(m1.get('r.f').is_state())
         self.assertEqual(m1.get('r.f').unit(), ms.get('r.f').unit())
         self.assertEqual(
-            m1.get('r.f').state_value(), ms.get('r.f').state_value())
+            m1.get('r.f').initial_value(), ms.get('r.f').initial_value())
         self.assertEqual(
             m1.get('r.f').rhs().code(),
             myokit.Multiply(ms.get('r.f').rhs(), s2ms).code())
@@ -1167,12 +1261,12 @@ def test_import_component_units(self):
         self.assertFalse(m1['y'] is ms['y'])
         self.assertEqual(m1.get('x.a').unit(), ms.get('x.a').unit())
         self.assertEqual(
-            m1.get('x.a').state_value(), ms.get('x.a').state_value())
+            m1.get('x.a').initial_value(), ms.get('x.a').initial_value())
         self.assertEqual(
             m1.get('x.a').rhs().code(),
             myokit.Multiply(ms.get('x.a').rhs(), s2ms).code())
         self.assertEqual(
-            m1.get('y.e').state_value(), ms.get('y.e').state_value())
+            m1.get('y.e').initial_value(), ms.get('y.e').initial_value())
         self.assertEqual(
             m1.get('y.e').rhs().code(),
             myokit.Multiply(ms.get('y.e').rhs(), s2ms).code())
@@ -1316,6 +1410,113 @@ def test_import_component_units(self):
             m1.import_component, ms['q'], var_map=vm, convert_units=True)
         self.assertTrue(ms.is_similar(ms_unaltered, True))
 
+    def test_initial_values(self):
+        # Tests :meth:`Model.initial_values`.
+
+        m = myokit.parse_model('''
+            [[model]]
+            c.x = 2
+            c.y = 1 + 2
+            c.z = 2 * c.p
+
+            [c]
+            t = 0 bind time
+            dot(x) = 0
+            dot(y) = 0
+            dot(z) = 0
+            p = 3
+            q = 4
+        ''')
+
+        # Test expression version
+        p = m.get('c.p')
+        x0 = m.initial_values()
+        self.assertEqual(x0[0], myokit.Number(2))
+        self.assertEqual(
+            x0[1], myokit.Plus(myokit.Number(1), myokit.Number(2)))
+        self.assertEqual(x0[2], myokit.Multiply(myokit.Number(2), p.lhs()))
+
+        # Test float version
+        self.assertEqual(m.initial_values(True), [2, 3, 6])
+
+        # Test deprecated alias
+        with WarningCollector() as w:
+            self.assertEqual(m.state(), m.initial_values(True))
+        self.assertIn('deprecated', w.text())
+
+        # Test cycles are detected before evaluation
+        p.set_rhs('1 / q')
+        m.get('c.q').set_rhs('1 + p')
+        m.initial_values()  # No evaluation, so no error
+        self.assertRaises(
+            myokit.CyclicalDependencyError, m.initial_values, as_floats=True)
+
+    def test_inits(self):
+        # Tests :meth:`Model.inits`
+
+        m = myokit.parse_model('''
+            [[model]]
+            c.x = 2
+            c.y = 1 + 2
+            c.z = 2 * c.p
+
+            [c]
+            t = 0 bind time
+            dot(x) = 0
+            dot(y) = 0
+            dot(z) = 0
+            p = 3
+        ''')
+        with WarningCollector() as w:
+            eqs = m.inits()
+        self.assertIn('deprecated', w.text())
+
+        E, N, M = myokit.Equation, myokit.Number, myokit.Multiply
+        self.assertEqual(list(eqs), [
+            E(myokit.Name(m.get('c.x')), N(2)),
+            E(myokit.Name(m.get('c.y')), myokit.Plus(N(1), N(2))),
+            E(myokit.Name(m.get('c.z')), M(N(2), m.get('c.p').lhs()))
+        ])
+
+    def test_is_similar(self):
+        # Check that equality takes both code() and unames into account
+
+        # Test without custom reserved names
+        m1 = myokit.load_model('example')
+        m2 = m1.clone()
+        self.assertIsInstance(m2, myokit.Model)
+        self.assertFalse(m1 is m2)
+        self.assertNotEqual(m1, m2)
+        self.assertNotEqual(m2, m1)
+        self.assertTrue(m1.is_similar(m2, False))
+        self.assertTrue(m1.is_similar(m2, True))
+        self.assertTrue(m2.is_similar(m1, False))
+        self.assertTrue(m2.is_similar(m1, True))
+        self.assertTrue(m1.is_similar(m1, True))
+        self.assertTrue(m2.is_similar(m2, False))
+
+        # Test with none-model
+        self.assertFalse(m1.is_similar(None))
+        self.assertFalse(m1.is_similar(m1.code()))
+
+        # Add reserved names
+        m1.reserve_unique_names('bertie')
+        self.assertFalse(m1.is_similar(m2))
+        m1.reserve_unique_names('clair')
+        self.assertFalse(m1.is_similar(m2))
+        m2.reserve_unique_names('clair', 'bertie')
+        self.assertTrue(m1.is_similar(m2))
+
+        # Add reserved name prefixes
+        m1.reserve_unique_name_prefix('aa', 'bb')
+        m1.reserve_unique_name_prefix('cc', 'dd')
+        self.assertFalse(m1.is_similar(m2))
+        m2.reserve_unique_name_prefix('aa', 'bb')
+        m2.reserve_unique_name_prefix('cc', 'ee')
+        self.assertFalse(m1.is_similar(m2))
+        m2.reserve_unique_name_prefix('cc', 'dd')
+        self.assertTrue(m1.is_similar(m2))
+
     def test_item_at_text_position(self):
         # Test :meth:`Model.item_at_text_position()`.
 
@@ -1454,17 +1655,21 @@ def test_load_save_state(self):
         # Test :meth:`Model.save_state()` and :meth:`Model.load_state()`.
 
         m = myokit.load_model('example')
-        s1 = m.state()
+        s1 = m.initial_values(True)
         with TemporaryDirectory() as d:
             path = d.path('state.csv')
-            m.save_state(path)
-            self.assertEqual(m.state(), s1)
+            with WarningCollector() as w:
+                m.save_state(path)
+            self.assertIn('deprecated', w.text())
+            self.assertEqual(m.initial_values(True), s1)
             sx = list(s1)
             sx[0] = 10
-            m.set_state(sx)
-            self.assertNotEqual(m.state(), s1)
-            m.load_state(path)
-            self.assertEqual(m.state(), s1)
+            m.set_initial_values(sx)
+            self.assertNotEqual(m.initial_values(True), s1)
+            with WarningCollector() as w:
+                m.load_state(path)
+            self.assertIn('deprecated', w.text())
+            self.assertEqual(m.initial_values(True), s1)
 
     def test_map_to_state(self):
         # Test :meth:`Model.map_to_state()`.
@@ -1624,7 +1829,7 @@ def test_remove_component(self):
     def test_remove_derivative_references(self):
         # Test the remove_derivative_references() method.
 
-        m0 = myokit.parse_model("""
+        m0 = myokit.parse_model('''
             [[model]]
             c.x = 0
             c.y = 1
@@ -1643,8 +1848,8 @@ def test_remove_derivative_references(self):
 
             [d]
             z = 3 * dot(c.x) / dot(c.y)
-            """)
-        m2 = myokit.parse_model("""
+            ''')
+        m2 = myokit.parse_model('''
             [[model]]
             c.x = 0
             c.y = 1
@@ -1665,7 +1870,7 @@ def test_remove_derivative_references(self):
 
             [d]
             z = 3 * c.dot_x_1 / c.dot_y
-            """)
+            ''')
 
         # Remove derivatives from m1
         m1 = m0.clone()
@@ -1700,31 +1905,6 @@ def test_remove_derivative_references(self):
         m1.remove_derivative_references()
         self.assertEqual(m1.get('c.dot_y').unit(), myokit.parse_unit('mV/ms'))
 
-    def test_remove_variable_with_alias(self):
-        # Test cloning of a variable with an alias after an add / remove event.
-
-        m = myokit.Model('AddRemoveClone')
-        c = m.add_component('c')
-        p = c.add_variable('p')
-        p.set_binding('time')
-        p.set_rhs(0)
-        q = c.add_variable('q')
-        q.set_rhs(12)
-        m.validate()    # Raises error if not ok
-        m.clone()       # Raises error if not ok
-        d = m.add_component('d')
-        d.add_alias('bert', p)
-        e = d.add_variable('e')
-        e.set_rhs('10 * bert')
-        m.validate()
-        m.clone()
-        d.add_alias('ernie', q)
-        m.validate()
-        m.clone()
-        c.remove_variable(q)
-        m.validate()
-        m.clone()   # Will raise error if alias isn't deleted
-
     def test_reorder_state(self):
         # Test :meth:`Model.reorder_state()`.
 
@@ -1830,6 +2010,70 @@ def test_sequence_interface(self):
         c = model['membrane']
         self.assertEqual(c.name(), 'membrane')
 
+    def test_set_initial_values(self):
+        # Tests :meth:`Model.set_initial_values()`.
+        m = myokit.load_model('example')
+        states = [v.qname() for v in m.states()]
+        values = m.initial_values(as_floats=True)
+
+        # Set with floats
+        values[0] += 2
+        m.set_initial_values(values)
+        for val, ref in zip(m.initial_values(True), values):
+            self.assertEqual(val, ref)
+
+        # Set with expressions
+        values[1] -= 3
+        m.set_initial_values([myokit.Number(v) for v in values])
+        for val, ref in zip(m.initial_values(True), values):
+            self.assertEqual(val, ref)
+
+        # Set with mix, including strings and expressions
+        values[2] += 10
+        mixed = list(values)
+        mixed[2] = myokit.Number(mixed[2])
+        mixed[3] = str(values[3])
+        m.set_initial_values(mixed)
+        for val, ref in zip(m.initial_values(True), values):
+            self.assertEqual(val, ref)
+
+        # Set with dict of floats
+        values = m.initial_values(as_floats=True)
+        s = dict(zip(states, values))
+        m.set_initial_values(s)
+        for val, ref in zip(m.initial_values(True), values):
+            self.assertEqual(val, ref)
+
+        # Set with one big string
+        values[0] += 5
+        s = dict(zip(states, values))
+        s = '\n'.join([str(a) + '=' + str(b) for a, b in s.items()])
+        m.set_initial_values(s)
+        for val, ref in zip(m.initial_values(True), values):
+            self.assertEqual(val, ref)
+
+        # Test that errors are detected in validation
+        values[0] = myokit.Name(m.get('membrane.V'))
+        m.set_initial_values(values)
+        self.assertRaisesRegex(myokit.IntegrityError, 'not const', m.validate)
+        values[0] = 'sqrt(3)'
+        m.set_initial_values(values)
+        m.validate()
+
+        # Deprecated alias
+        values = m.initial_values(as_floats=True)
+        values[0] += 10
+        with WarningCollector() as w:
+            m.set_state(values)
+        self.assertIn('deprecated', w.text())
+        for val, ref in zip(m.initial_values(True), values):
+            self.assertEqual(val, ref)
+
+        # Wrong number of elements
+        self.assertRaisesRegex(
+            ValueError, 'Wrong number of initial values',
+            m.set_initial_values, [1, 2, 3])
+
     def test_show_evaluation_of(self):
         # Test :meth:`Model.show_evaluation_of(variable)`.
         # Depends mostly on `references()`, and `code()` methods.
@@ -1867,6 +2111,16 @@ def test_show_evaluation_of(self):
         self.assertIn('Literal constant', e)
         self.assertEqual(len(e.splitlines()), 7)
 
+        # Test with expressions in initial values
+        m.get('ina.m').set_initial_value('1 / sqrt(4)')
+        e = m.show_evaluation_of('ina.m')
+        self.assertIn('Initial value = 1 / sqrt(4)', e)
+        self.assertEqual(len(e.splitlines()), 16)
+        m.get('ina.m').set_initial_value('1 / sqrt(ina.gNa)')
+        e = m.show_evaluation_of('ina.m')
+        self.assertIn('Initial value = 1 / sqrt(ina.gNa)', e)
+        self.assertIn('              = 2.5000000', e)
+
         # Test with nothing similar
         m = myokit.Model()
         self.assertRaises(Exception, m.show_evaluation_of, 'Hello')
@@ -1886,7 +2140,7 @@ def test_show_line_of(self):
         m = myokit.load_model('example')
         v = m.get('ina.INa')
         e = m.show_line_of(v)
-        self.assertIn('Defined on line 91', e)
+        self.assertIn('Defined on line 90', e)
         self.assertIn('Intermediary variable', e)
         self.assertEqual(len(e.splitlines()), 4)
 
@@ -1899,7 +2153,7 @@ def test_show_line_of(self):
         self.assertEqual(len(e.splitlines()), 3)
 
         # 'raw' version
-        self.assertEqual(m.show_line_of(v, raw=True), 91)
+        self.assertEqual(m.show_line_of(v, raw=True), 90)
         self.assertIsNone(m2.show_line_of(v2, raw=True))
 
     def test_str(self):
@@ -2051,6 +2305,56 @@ def test_unique_names_2(self):
             ValueError, 'prepend cannot start with prefix',
             m.reserve_unique_name_prefix, 'x', 'x')
 
+    def test_unused_and_cycles(self):
+        # Test unused variable and cycle detection.
+
+        m = myokit.Model('LotkaVolterra')
+        c0 = m.add_component('c0')
+        t = c0.add_variable('time')
+        t.set_rhs(myokit.Number(0))
+        t.set_binding('time')
+        c1 = m.add_component('c1')
+        m.add_component('c2')
+        c1_a = c1.add_variable('a')
+        c1_b = c1.add_variable('b')
+        c1_a.promote(1.0)
+        c1_a.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Number(0.5)))
+        c1_b.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Number(1.0)))
+        # b is unused, test if found
+        m.validate()
+        w = m.warnings()
+        self.assertEqual(len(w), 1)
+        self.assertEqual(type(w[0]), myokit.UnusedVariableError)
+        # b is used by c, c is unused, test if found
+        c1_c = c1.add_variable('c')
+        c1_c.set_rhs(myokit.Name(c1_b))
+        m.validate()
+        w = m.warnings()
+        self.assertEqual(len(w), 2)
+        self.assertEqual(type(w[0]), myokit.UnusedVariableError)
+        self.assertEqual(type(w[1]), myokit.UnusedVariableError)
+        # Test 1:1 cycle
+        c1_b.set_rhs(myokit.Name(c1_b))
+        self.assertRaises(myokit.CyclicalDependencyError, m.validate)
+        # Test longer cycles
+        c1_b.set_rhs(myokit.Multiply(myokit.Number(10), myokit.Name(c1_c)))
+        self.assertRaises(myokit.CyclicalDependencyError, m.validate)
+        # Reset
+        c1_b.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Number(1.0)))
+        m.validate()
+        # Test cycle involving state variable
+        c1_a.set_rhs(myokit.Name(c1_b))
+        m.validate()
+        c1_b.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Name(c1_b)))
+        self.assertRaises(myokit.CyclicalDependencyError, m.validate)
+        c1_b.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Name(c1_c)))
+        c1_c.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Number(3)))
+        m.validate()
+        w = m.warnings()
+        self.assertEqual(len(w), 0)
+        c1_c.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Name(c1_b)))
+        self.assertRaises(myokit.CyclicalDependencyError, m.validate)
+
     def test_validate_and_remove_unused_variables(self):
         # Test :class:`Model.validate` with ``remove_unused_variables=True``.
 
diff --git a/myokit/tests/test_model_building.py b/myokit/tests/test_model_building.py
index a2f85deec..7d21de6d1 100755
--- a/myokit/tests/test_model_building.py
+++ b/myokit/tests/test_model_building.py
@@ -2,6 +2,9 @@
 #
 # Tests the model building API.
 #
+# NOTE: THESE TESTS SHOULD SLOWLY BE MERGED INTO test_model, test_variable, etc
+#
+#
 # This file is part of Myokit.
 # See http://myokit.org for copyright, sharing, and licensing details.
 #
@@ -33,13 +36,23 @@ def test_model_creation(self):
         # Create a model
         m = myokit.Model('LotkaVolterra')
 
-        # Add the first component
+        # Add engine component
+        E = m.add_component('engine')
+        self.assertEqual(len(m), 1)
+        time = E.add_variable('time')
+        time.set_rhs(0)
+        self.assertIsNone(time.binding())
+        time.set_binding('time')
+        self.assertIsNotNone(time.binding())
+
+        # Add the first test component
         X = m.add_component('X')
+        self.assertNotEqual(E, X)
         self.assertEqual(X.qname(), 'X')
         self.assertEqual(X.parent(), m)
         self.assertIsInstance(X, myokit.Component)
         self.assertIn(X.qname(), m)
-        self.assertEqual(len(m), 1)
+        self.assertEqual(len(m), 2)
 
         # Add variable a
         self.assertFalse(X.has_variable('a'))
@@ -133,10 +146,32 @@ def test_model_creation(self):
         self.assertFalse(x.is_constant())
         self.assertEqual(x.lhs(), myokit.Derivative(myokit.Name(x)))
 
+        # Set number initial value
+        x.demote()
+        x.promote(1)
+        self.assertEqual(x.initial_value(as_float=True), 1)
+
+        # Set constant expression initial value
+        x.set_initial_value(myokit.Name(b))
+        self.assertEqual(x.initial_value(as_float=True), 2)
+        self.assertEqual(x.initial_value(), myokit.Name(b))
+
+        # Non-constant expression initial value can be set, but isn't valid
+        m.validate()
+        b.promote()
+        self.assertRaisesRegex(
+            myokit.IntegrityError, 'not constant', m.validate)
+        b.demote()
+
+        # Set literal valued expression initial value
+        x.demote()
+        x.promote('1 + 2')
+        self.assertEqual(x.initial_value(as_float=True), 3)
+
         # Add second component, variables
         Y = m.add_component('Y')
         self.assertNotEqual(X, Y)
-        self.assertEqual(len(m), 2)
+        self.assertEqual(len(m), 3)
         c = Y.add_variable('c')
         c.set_rhs(myokit.Minus(myokit.Name(a), myokit.Number(1)))
         d = Y.add_variable('d')
@@ -152,7 +187,7 @@ def test_model_creation(self):
                 myokit.Name(y)
             )
         ))
-        x.set_state_value(10)
+        x.set_initial_value(10)
         self.assertEqual(x.rhs().code(), 'X.a * X.x - X.b * X.x * Y.y')
         y.set_rhs(myokit.Plus(
             myokit.Multiply(
@@ -163,14 +198,14 @@ def test_model_creation(self):
                 myokit.Name(y)
             )
         ))
-        y.set_state_value(5)
+        y.set_initial_value(5)
         self.assertEqual(y.rhs().code(), '-Y.c * Y.y + Y.d * X.x * Y.y')
 
         # Add ano component, variables
         Z = m.add_component('Z')
         self.assertNotEqual(X, Z)
         self.assertNotEqual(Y, Z)
-        self.assertEqual(len(m), 3)
+        self.assertEqual(len(m), 4)
         t = Z.add_variable('total')
         self.assertEqual(t.name(), 'total')
         self.assertEqual(t.qname(), 'Z.total')
@@ -185,18 +220,6 @@ def test_model_creation(self):
         self.assertEqual(t.rhs().code(Y), 'X.x + y')
         self.assertEqual(t.rhs().code(Z), 'X.x + Y.y')
 
-        # Add engine component
-        E = m.add_component('engine')
-        self.assertNotEqual(X, E)
-        self.assertNotEqual(Y, E)
-        self.assertNotEqual(Z, E)
-        self.assertEqual(len(m), 4)
-        time = E.add_variable('time')
-        time.set_rhs(0)
-        self.assertIsNone(time.binding())
-        time.set_binding('time')
-        self.assertIsNotNone(time.binding())
-
         # Check state
         state = [i for i in m.states()]
         self.assertEqual(len(state), 2)
@@ -355,16 +378,6 @@ def has(var, *dps):
         code2 = m.clone().code()
         self.assertEqual(code1, code2)
 
-    def test_resolve(self):
-        # Test if an error is raised when a variable can't be resolved.
-
-        m = myokit.Model('Resolve')
-        c = m.add_component('c')
-        p = c.add_variable('p')
-        q = c.add_variable('q')
-        p.set_rhs('10 * q')
-        self.assertRaises(myokit.ParseError, q.set_rhs, '10 * r')
-
     def test_scope(self):
         # Test if illegal references are detected.
 
@@ -436,56 +449,6 @@ def test_invalid_names(self):
         self.assertRaises(myokit.InvalidNameError, v1.add_variable, 'not')
         self.assertRaises(myokit.InvalidNameError, v1.add_variable, 'in')
 
-    def test_unused_and_cycles(self):
-        # Test unused variable and cycle detection.
-
-        m = myokit.Model('LotkaVolterra')
-        c0 = m.add_component('c0')
-        t = c0.add_variable('time')
-        t.set_rhs(myokit.Number(0))
-        t.set_binding('time')
-        c1 = m.add_component('c1')
-        m.add_component('c2')
-        c1_a = c1.add_variable('a')
-        c1_b = c1.add_variable('b')
-        c1_a.promote(1.0)
-        c1_a.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Number(0.5)))
-        c1_b.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Number(1.0)))
-        # b is unused, test if found
-        m.validate()
-        w = m.warnings()
-        self.assertEqual(len(w), 1)
-        self.assertEqual(type(w[0]), myokit.UnusedVariableError)
-        # b is used by c, c is unused, test if found
-        c1_c = c1.add_variable('c')
-        c1_c.set_rhs(myokit.Name(c1_b))
-        m.validate()
-        w = m.warnings()
-        self.assertEqual(len(w), 2)
-        self.assertEqual(type(w[0]), myokit.UnusedVariableError)
-        self.assertEqual(type(w[1]), myokit.UnusedVariableError)
-        # Test 1:1 cycle
-        c1_b.set_rhs(myokit.Name(c1_b))
-        self.assertRaises(myokit.CyclicalDependencyError, m.validate)
-        # Test longer cycles
-        c1_b.set_rhs(myokit.Multiply(myokit.Number(10), myokit.Name(c1_c)))
-        self.assertRaises(myokit.CyclicalDependencyError, m.validate)
-        # Reset
-        c1_b.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Number(1.0)))
-        m.validate()
-        # Test cycle involving state variable
-        c1_a.set_rhs(myokit.Name(c1_b))
-        m.validate()
-        c1_b.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Name(c1_b)))
-        self.assertRaises(myokit.CyclicalDependencyError, m.validate)
-        c1_b.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Name(c1_c)))
-        c1_c.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Number(3)))
-        m.validate()
-        w = m.warnings()
-        self.assertEqual(len(w), 0)
-        c1_c.set_rhs(myokit.Multiply(myokit.Name(c1_a), myokit.Name(c1_b)))
-        self.assertRaises(myokit.CyclicalDependencyError, m.validate)
-
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/myokit/tests/test_parsing.py b/myokit/tests/test_parsing.py
index 4bc081058..1d91bf261 100755
--- a/myokit/tests/test_parsing.py
+++ b/myokit/tests/test_parsing.py
@@ -223,6 +223,7 @@ class PhasedParseTest(unittest.TestCase):
     """
     Tests several phases of parsing.
     """
+
     def test_segment_parsing(self):
         # Test :meth:`parse_model()`.
         from myokit._parsing import parse
@@ -1229,7 +1230,12 @@ def test_format_parse_error(self):
 
 
 class ModelParseTest(unittest.TestCase):
+    """
+    Test parsing models.
+    """
+
     def test_model_creation(self):
+        # Basic model creation via the parser
         m = myokit.load_model(os.path.join(DIR_DATA, 'lr-1991.mmt'))
 
         # Test components
@@ -1263,36 +1269,14 @@ def test_model_creation(self):
             1,
             0.0057,
             0.0002]
+        self.assertEqual(len(values), m.count_states())
+
+        # Test initial value parsing
         out = ', '.join([str(x) for x in m.states()])
         ref = ', '.join(states)
         self.assertEqual(ref, out)
-        for k, eq in enumerate(m.inits()):
-            self.assertEqual(eq.rhs.eval(), values[k])
-
-        # Test state parsing / setting
-        m.set_state(values)
-        for k, eq in enumerate(m.inits()):
-            self.assertEqual(eq.rhs.eval(), values[k])
-        s = dict(zip(states, values))
-        m.set_state(s)
-        for k, eq in enumerate(m.inits()):
-            self.assertEqual(eq.rhs.eval(), values[k])
-        s = '\n'.join([str(a) + '=' + str(b) for a, b in s.items()])
-        m.set_state(s)
-        for k, eq in enumerate(m.inits()):
-            self.assertEqual(eq.rhs.eval(), values[k])
-
-        # Test cloning
-        try:
-            m2 = m.clone()
-        except Exception as e:
-            s = m.code(line_numbers=True)
-            print('\n')
-            print(s)
-            print('-' * 80)
-            print(myokit.format_parse_error(e, s.splitlines()))
-            raise e
-        self.assertEqual(m.code(), m2.code())
+        for val, ref in zip(m.initial_values(as_floats=True), values):
+            self.assertEqual(val, ref)
 
     def test_unresolved_reference_error(self):
         # Test unresolved reference errors.
@@ -1354,19 +1338,21 @@ def test_piecewise(self):
         # Test evaluation
         x = m.get('test.x')
         y = m.get('test.y')
-        s = m.state()
+        s = m.initial_values(as_floats=True)
         i = m.get('membrane.V').indice()
         # Test x, xo, y, yo
+
+        v = myokit.Name(m.get('membrane.V'))
         s[i] = -80
-        m.set_state(s)
+        m.set_initial_values(s)
         self.assertEqual(x.rhs().eval(), 3)
         self.assertEqual(y.rhs().eval(), 2)
         s[i] = -10
-        m.set_state(s)
+        m.set_initial_values(s)
         self.assertEqual(x.rhs().eval(), 2)
         self.assertEqual(y.rhs().eval(), 2)
         s[i] = 30
-        m.set_state(s)
+        m.set_initial_values(s)
         self.assertEqual(x.rhs().eval(), 1)
         self.assertEqual(y.rhs().eval(), 1)
         # Test code() output by cloning
@@ -1387,13 +1373,50 @@ def test_initial_values(self):
             dot(p) = 1
             dot(q) = 2
             """
-        myokit.parse(code)
+        m = myokit.parse_model(code)
+        self.assertEqual(m.get('c.p').initial_value(), myokit.Number(1))
+        self.assertEqual(m.get('c.q').initial_value().code(), '10 * 2')
+
+        # initial expression
+        code = """
+            [[model]]
+            c.q = 10 * 2 + c.p
+            c.r = c.p + c.a
+
+            [engine]
+            time = 0 bind time
+
+            [c]
+            p = 0.5
+            dot(q) = 2
+            dot(r) = 1
+            a = p + 1
+            """
+        m = myokit.parse_model(code)
+        self.assertEqual(
+            m.get('c.q').initial_value().code(), '10 * 2 + c.p')
+
+        # initial expression with just literals
+        code = """
+            [[model]]
+            c.q = 10 * exp(2)
+
+            [engine]
+            time = 0 bind time
 
-        # Non-literal value
+            [c]
+            p = 0.5
+            dot(q) = 2
+            """
+        m = myokit.parse_model(code)
+        self.assertEqual(
+            m.get('c.q').initial_value().code(), '10 * exp(2)')
+
+        # non-existent variable in initial value
         code = """
             [[model]]
             c.p = 1.0
-            c.q = 10 * 2 + b
+            c.q = 10 * 2 + c.b
 
             [engine]
             time = 0 bind time
@@ -1402,7 +1425,82 @@ def test_initial_values(self):
             dot(p) = 1
             dot(q) = 2
             """
-        self.assertRaises(myokit.ParseError, myokit.parse, code)
+        self.assertRaisesRegex(
+            myokit.ParseError, 'Unresolved reference', myokit.parse, code)
+
+        # initial expression using child variables should fail
+        code = """
+            [[model]]
+            c.q = a
+
+            [engine]
+            time = 0 bind time
+
+            [c]
+            p = 0.5
+            dot(q) = 2
+                a = 10 * 2 + p
+            """
+        self.assertRaisesRegex(
+            myokit.ParseError, 'Unresolved reference', myokit.parse, code)
+
+        # non-constant initial value
+        code = """
+            [[model]]
+            c.p = engine.time
+
+            [engine]
+            time = 0 bind time
+
+            [c]
+            dot(p) = 1
+            """
+        self.assertRaisesRegex(
+            myokit.ParseError, 'is not constant', myokit.parse, code)
+
+        # indirect non-constant initial value
+        code = """
+            [[model]]
+            c.p = c.r
+            c.q = 1
+
+            [engine]
+            time = 0 bind time
+
+            [c]
+            dot(p) = 1
+            dot(q) = 2
+            r = 2 * q
+            """
+        self.assertRaisesRegex(
+            myokit.ParseError, 'is not constant', myokit.parse, code)
+
+        # looks like it shouldn't be constant, but it is!
+        code = """
+            [[model]]
+            c.p = dot(c.p)
+
+            [engine]
+            time = 0 bind time
+
+            [c]
+            dot(p) = 0
+            """
+        # Treat all state variables as non-constant
+        self.assertRaisesRegex(
+            myokit.ParseError, 'is not constant', myokit.parse, code)
+
+        # Initial value cycles are caught in validation
+        code = """
+            [[model]]
+            x.x = x.x
+
+            [x]
+            dot(x) = x
+            t = 0 bind time
+            """
+        self.assertRaisesRegex(
+            myokit.ParseError, 'is not constant', myokit.parse, code)
 
     def test_aliases(self):
         code = """
diff --git a/myokit/tests/test_sbml_api.py b/myokit/tests/test_sbml_api.py
index 579d1af90..2a3b0e232 100755
--- a/myokit/tests/test_sbml_api.py
+++ b/myokit/tests/test_sbml_api.py
@@ -1364,7 +1364,7 @@ def test_parameter_values_rate(self):
         pp = mm.get('myokit.param')
         self.assertTrue(pp.is_state())
         self.assertEqual(pp.rhs(), myokit.Number(3.2))
-        self.assertEqual(pp.state_value(), 0)
+        self.assertEqual(pp.initial_value(as_float=True), 0)
 
         # With initial value
         p.set_initial_value(myokit.Number(1))
@@ -1372,7 +1372,7 @@ def test_parameter_values_rate(self):
         pp = mm.get('myokit.param')
         self.assertTrue(pp.is_state())
         self.assertEqual(pp.rhs(), myokit.Number(3.2))
-        self.assertEqual(pp.state_value(), 1)
+        self.assertEqual(pp.initial_value(as_float=True), 1)
 
         # Test units
         self.assertIsNone(pp.unit())
@@ -1542,13 +1542,13 @@ def test_species_value_rate(self):
         ms = mm.get('comp.spec_1_amount')
         self.assertTrue(ms.is_state())
         self.assertEqual(ms.rhs(), myokit.Number(3))
-        self.assertEqual(ms.state_value(), 0)
+        self.assertEqual(ms.initial_value(as_float=True), 0)
         s1.set_initial_value(myokit.Number(7))
         mm = m.myokit_model()
         ms = mm.get('comp.spec_1_amount')
         self.assertTrue(ms.is_state())
         self.assertEqual(ms.rhs(), myokit.Number(3))
-        self.assertEqual(ms.state_value(), 7)
+        self.assertEqual(ms.initial_value(as_float=True), 7)
 
         # Species in concentration
         s2 = m.add_species(c, 'spec_2', is_amount=False)
@@ -1562,7 +1562,7 @@ def test_species_value_rate(self):
         sa = mm.get('comp.spec_2_amount')
         self.assertTrue(sa.is_state())
         self.assertEqual(sa.rhs().code(), '4 * comp.size')
-        self.assertEqual(sa.state_value(), 0)
+        self.assertEqual(sa.initial_value(as_float=True), 0)
 
         # I: Set compartment size
         m = sbml.Model()
@@ -1577,13 +1577,13 @@ def test_species_value_rate(self):
         ms = mm.get('comp.spec_1_amount')
         self.assertTrue(ms.is_state())
         self.assertEqual(ms.rhs(), myokit.Number(3))
-        self.assertEqual(ms.state_value(), 0)
+        self.assertEqual(ms.initial_value(as_float=True), 0)
         s1.set_initial_value(myokit.Number(7))
         mm = m.myokit_model()
         ms = mm.get('comp.spec_1_amount')
         self.assertTrue(ms.is_state())
         self.assertEqual(ms.rhs(), myokit.Number(3))
-        self.assertEqual(ms.state_value(), 7)
+        self.assertEqual(ms.initial_value(as_float=True), 7)
 
         # Species in concentration
         s2 = m.add_species(c, 'spec_2', is_amount=False)
@@ -1599,7 +1599,7 @@ def test_species_value_rate(self):
         sa = mm.get('comp.spec_2_amount')
         self.assertTrue(sa.is_state())
         self.assertEqual(sa.rhs().code(), '4 * comp.size')
-        self.assertEqual(sa.state_value(), 6 * 2)
+        self.assertEqual(sa.initial_value(as_float=True), 6 * 2)
 
     def test_species_units(self):
         # Tests whether species units are set properly.
@@ -1968,7 +1968,7 @@ def test_reaction_stoichiometry_parameter(self):
         self.assertEqual(var.eval(), 5)
 
         var = mm.get('c.sr2')
-        self.assertEqual(var.state_value(), 0)
+        self.assertEqual(var.initial_value(as_float=True), 0)
         self.assertEqual(var.eval(), 3.82)
 
     def test_reaction_stoichiometries_exist(self):
@@ -2061,7 +2061,7 @@ def test_reaction_stoichiometry_values(self):
         self.assertEqual(var.eval(), 15.23)
 
         var = mm.get('c.sr2')
-        self.assertEqual(var.state_value(), 3.5)
+        self.assertEqual(var.initial_value(as_float=True), 3.5)
         self.assertEqual(var.eval(), 9.23)
 
         # Bad value
diff --git a/myokit/tests/test_simulation_1d.py b/myokit/tests/test_simulation_1d.py
index e64d756bd..50978451e 100755
--- a/myokit/tests/test_simulation_1d.py
+++ b/myokit/tests/test_simulation_1d.py
@@ -31,6 +31,7 @@ class Simulation1dTest(unittest.TestCase):
     """
     Test the non-parallel 1d simulation.
     """
+
     def test_basic(self):
         # Test basic usage.
 
@@ -42,23 +43,24 @@ def test_basic(self):
         s = myokit.Simulation1d(m, p, ncells=ncells)
         s.set_step_size(0.05)
 
+        x0 = m.initial_values(True)
         self.assertEqual(s.time(), 0)
-        self.assertEqual(s.state(0), m.state())
-        self.assertEqual(s.default_state(0), m.state())
+        self.assertEqual(s.state(0), x0)
+        self.assertEqual(s.default_state(0), x0)
         d = s.run(5, log_interval=1)
         self.assertEqual(s.time(), 5)
-        self.assertNotEqual(s.state(0), m.state())
-        self.assertEqual(s.default_state(0), m.state())
+        self.assertNotEqual(s.state(0), x0)
+        self.assertEqual(s.default_state(0), x0)
 
         # Test full state getting and reset
-        self.assertEqual(s.default_state(), m.state() * ncells)
-        self.assertNotEqual(s.state(), m.state() * ncells)
+        self.assertEqual(s.default_state(), x0 * ncells)
+        self.assertNotEqual(s.state(), x0 * ncells)
         s.reset()
         self.assertEqual(s.state(), s.default_state())
 
         # Test pre updates the default state.
         s.pre(1)
-        self.assertNotEqual(s.default_state(0), m.state())
+        self.assertNotEqual(s.default_state(0), x0)
 
         # Test running without a protocol
         s.set_protocol(None)
@@ -103,6 +105,31 @@ def test_basic(self):
         s.set_time(100)
         self.assertEqual(s.time(), 100)
 
+    def test_initial_value_expressions(self):
+        # Test if initial value expressions are converted to floats
+        m = myokit.parse_model('''
+            [[model]]
+            c.x = 1 + sqrt(3)
+            c.y = 1 / c.p
+            c.z = 3
+
+            [c]
+            t = 0 bind time
+            dot(x) = 1 label membrane_potential
+            dot(y) = 2
+            dot(z) = 3
+            p = log(3)
+            q = 0 bind diffusion_current
+        ''')
+        s = myokit.Simulation1d(m, ncells=2)
+        x = s.state()
+        self.assertIsInstance(x[0], float)
+        self.assertIsInstance(x[1], float)
+        self.assertIsInstance(x[2], float)
+        self.assertEqual(x[:3], x[3:])
+        self.assertEqual(x, m.initial_values(True) * 2)
+        self.assertEqual(x, s.default_state())
+
     def test_with_progress_reporter(self):
         # Test running with a progress reporter.
         m, p, _ = myokit.load(os.path.join(DIR_DATA, 'lr-1991.mmt'))
@@ -130,7 +157,8 @@ def test_set_state(self):
         n = 4
 
         s = myokit.Simulation1d(m, p, n)
-        self.assertEqual(s.state(), m.state() * n)
+        x0 = m.initial_values(True)
+        self.assertEqual(s.state(), x0 * n)
 
         # Test setting a full state
         sx = [0] * 8 * n
@@ -142,8 +170,8 @@ def test_set_state(self):
         sx = [0] * 8
         s.set_state(sx)
         self.assertEqual(s.state(), sx * n)
-        s.set_state(m.state())
-        self.assertEqual(s.state(), m.state() * n)
+        s.set_state(x0)
+        self.assertEqual(s.state(), x0 * n)
 
         # Test setting a single state
         j = 1
@@ -152,7 +180,7 @@ def test_set_state(self):
             if i == j:
                 self.assertEqual(s.state(i), sx)
             else:
-                self.assertEqual(s.state(i), m.state())
+                self.assertEqual(s.state(i), x0)
 
         # Invalid cell index
         s.set_state(sx, 0)
@@ -167,7 +195,7 @@ def test_set_default_state(self):
         n = 4
 
         s = myokit.Simulation1d(m, p, n)
-        self.assertEqual(s.state(), m.state() * n)
+        self.assertEqual(s.state(), m.initial_values(True) * n)
 
         # Test setting a full state
         sx = [0] * 8 * n
@@ -179,8 +207,9 @@ def test_set_default_state(self):
         sx = [0] * 8
         s.set_default_state(sx)
         self.assertEqual(s.default_state(), sx * n)
-        s.set_default_state(m.state())
-        self.assertEqual(s.default_state(), m.state() * n)
+        sy = m.initial_values(True)
+        s.set_default_state(sy)
+        self.assertEqual(s.default_state(), sy * n)
 
         # Test setting a single state
         j = 1
@@ -189,7 +218,7 @@ def test_set_default_state(self):
             if i == j:
                 self.assertEqual(s.default_state(i), sx)
             else:
-                self.assertEqual(s.default_state(i), m.state())
+                self.assertEqual(s.default_state(i), sy)
 
         # Invalid cell index
         s.set_default_state(sx, 0)
diff --git a/myokit/tests/test_simulation_cvodes.py b/myokit/tests/test_simulation_cvodes.py
index 097252e47..eed42eb86 100755
--- a/myokit/tests/test_simulation_cvodes.py
+++ b/myokit/tests/test_simulation_cvodes.py
@@ -154,6 +154,30 @@ def test_in_parts(self):
         self.assertNotEqual(e['engine.time'][n - 1], e['engine.time'][n])
         self.assertGreater(e['engine.time'][n], e['engine.time'][n - 1])
 
+    def test_initial_value_expressions(self):
+        # Test if initial value expressions are converted to floats
+
+        m = myokit.parse_model('''
+            [[model]]
+            c.x = 1 + sqrt(3)
+            c.y = 1 / c.p
+            c.z = 3
+
+            [c]
+            t = 0 bind time
+            dot(x) = 1
+            dot(y) = 2
+            dot(z) = 3
+            p = log(3)
+        ''')
+        s = myokit.Simulation(m)
+        x = s.state()
+        self.assertIsInstance(x[0], float)
+        self.assertIsInstance(x[1], float)
+        self.assertIsInstance(x[2], float)
+        self.assertEqual(x, m.initial_values(True))
+        self.assertEqual(x, s.default_state())
+
     def test_pacing_values_at_event_transitions(self):
         # Tests the value of the pacing signal at event transitions
 
@@ -335,12 +359,12 @@ def test_eval_derivatives(self):
         self.sim.set_state(s1)
         self.assertEqual(d1, self.sim.eval_derivatives())
 
-    def test_sensitivites_initial(self):
+    def test_sensitivities_initial(self):
         # Test setting initial sensitivity values.
 
         m = myokit.parse_model('''
             [[model]]
-            e.y = 2.3
+            e.y = 1 + 1.3
 
             [e]
             t = 0 bind time
@@ -352,6 +376,26 @@ def test_sensitivites_initial(self):
         #TODO: Test results
         s = myokit.Simulation(m, sensitivities=(['e.y'], ['e.p', 'init(e.y)']))
 
+        # Warn if a parameter sensitivity won't be picked up.
+        m = myokit.parse_model('''
+            [[model]]
+            c.x = 1 / c.p
+
+            [c]
+            t = 0 bind time
+            p = 1 / q
+            q = 5
+            r = 3
+            dot(x) = 2 + r
+            ''')
+        m.validate()
+        s = (['c.x'], ['c.q'])
+        self.assertRaisesRegex(
+            NotImplementedError, 'respect to parameters used in initial',
+            myokit.Simulation, m, sensitivities=s)
+        s = (['c.x'], ['c.r'])
+        s = myokit.Simulation(m, sensitivities=s)
+
         # Test bad initial matrix
         self.assertRaisesRegex(
             ValueError, 'None or a list',
@@ -539,8 +583,8 @@ def test_sensitivities_with_derivatives(self):
             ''')
         m.validate()
 
-        x0 = m.get('e.x').state_value()
-        y0 = m.get('e.y').state_value()
+        x0 = m.get('e.x').initial_value(True)
+        y0 = m.get('e.y').initial_value(True)
         p = m.get('e.p').eval()
         q = m.get('e.q').eval()
         r = m.get('e.r').eval()
diff --git a/myokit/tests/test_simulation_fiber_tissue.py b/myokit/tests/test_simulation_fiber_tissue.py
index 0cec629d5..6a453e718 100755
--- a/myokit/tests/test_simulation_fiber_tissue.py
+++ b/myokit/tests/test_simulation_fiber_tissue.py
@@ -535,7 +535,7 @@ def test_fiber_state(self):
         m = self.mf.count_states()
         nx, ny = self.nfx, self.nfy
 
-        sm = self.mf.state()
+        sm = self.mf.initial_values(True)
         for i in range(nx):
             for j in range(ny):
                 self.assertEqual(sm, self.s1.fiber_state(i, j))
@@ -622,6 +622,42 @@ def test_fiber_state(self):
         finally:
             self.s1.reset()
 
+    def test_initial_value_expressions(self):
+        # Test if initial value expressions are converted to floats
+        m = myokit.parse_model('''
+            [[model]]
+            c.x = 1 + sqrt(3)
+            c.y = 1 / c.p
+            c.z = 3
+
+            [c]
+            t = 0 bind time
+            dot(x) = 1 label membrane_potential
+                in [mV]
+            dot(y) = 2
+            dot(z) = 3
+            p = log(3)
+            q = 0 bind diffusion_current
+                in [A/F]
+        ''')
+        s = myokit.FiberTissueSimulation(
+            m, m, ncells_fiber=(1, 1), ncells_tissue=(3, 1))
+        x = s.fiber_state()
+        self.assertIsInstance(x[0], float)
+        self.assertIsInstance(x[1], float)
+        self.assertIsInstance(x[2], float)
+        self.assertEqual(x, m.initial_values(True))
+        self.assertEqual(x, s.default_fiber_state())
+
+        x = s.tissue_state()
+        self.assertIsInstance(x[0], float)
+        self.assertIsInstance(x[1], float)
+        self.assertIsInstance(x[2], float)
+        self.assertEqual(x[:3], x[3:6])
+        self.assertEqual(x[:3], x[6:])
+        self.assertEqual(x, m.initial_values(True) * 3)
+        self.assertEqual(x, s.default_tissue_state())
+
     def test_tissue_state(self):
         # Test the set_tissue_state and set_default_tissue_state methods
 
@@ -629,7 +665,7 @@ def test_tissue_state(self):
         m = self.mt.count_states()
         nx, ny = self.ntx, self.nty
 
-        sm = self.mt.state()
+        sm = self.mt.initial_values(True)
         for i in range(nx):
             for j in range(ny):
                 self.assertEqual(sm, self.s1.tissue_state(i, j))
diff --git a/myokit/tests/test_simulation_opencl.py b/myokit/tests/test_simulation_opencl.py
index 60c908330..75708c849 100755
--- a/myokit/tests/test_simulation_opencl.py
+++ b/myokit/tests/test_simulation_opencl.py
@@ -666,6 +666,31 @@ def test_diffusion_free(self):
         self.assertTrue(np.all(d['membrane.V', 0] == d['membrane.V', 1]))
         self.assertTrue(np.all(d['membrane.V', 0] == d['membrane.V', 2]))
 
+    def test_initial_value_expressions(self):
+        # Test if initial value expressions are converted to floats
+        m = myokit.parse_model('''
+            [[model]]
+            c.x = 1 + sqrt(3)
+            c.y = 1 / c.p
+            c.z = 3
+
+            [c]
+            t = 0 bind time
+            dot(x) = 1 label membrane_potential
+            dot(y) = 2
+            dot(z) = 3
+            p = log(3)
+            q = 0 bind diffusion_current
+        ''')
+        s = myokit.SimulationOpenCL(m, ncells=2)
+        x = s.state()
+        self.assertIsInstance(x[0], float)
+        self.assertIsInstance(x[1], float)
+        self.assertIsInstance(x[2], float)
+        self.assertEqual(x[:3], x[3:])
+        self.assertEqual(x, m.initial_values(True) * 2)
+        self.assertEqual(x, s.default_state())
+
     def test_neighbours_0d(self):
         # Test listing neighbours in a 0d simulation
 
@@ -1180,7 +1205,7 @@ def test_set_state_1d(self):
         m = 8
         n = 10
         self.s1.reset()
-        sm = self.m.state()
+        sm = self.m.initial_values(True)
         ss = [self.s1.state(x) for x in range(n)]
         for si in ss:
             self.assertEqual(sm, si)
@@ -1277,7 +1302,7 @@ def test_set_state_2d(self):
         m = 8
         nx, ny = 4, 3
         self.s2.reset()
-        sm = self.m.state()
+        sm = self.m.initial_values(True)
         for i in range(nx):
             for j in range(ny):
                 self.assertEqual(sm, self.s2.state(i, j))
diff --git a/myokit/tests/test_variable.py b/myokit/tests/test_variable.py
index 3510aef82..60e973cd6 100755
--- a/myokit/tests/test_variable.py
+++ b/myokit/tests/test_variable.py
@@ -8,10 +8,14 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 
-import myokit
-import numpy as np
 import unittest
 
+import numpy as np
+
+import myokit
+
+from myokit.tests import WarningCollector
+
 
 # Unit testing in Python 2 and 3
 try:
@@ -24,6 +28,7 @@ class VariableTest(unittest.TestCase):
     """
     Tests parts of :class:`myokit.Variable`.
     """
+
     def test_clamp(self):
         # Tests clamping a variable to a fixed value
         m = myokit.parse_model("""
@@ -92,7 +97,7 @@ def test_clamp(self):
         # Clamp simple state
         g = m.get('z.g')
         self.assertTrue(g.is_state())
-        self.assertTrue(g.state_value() != 10)
+        self.assertNotEqual(g.initial_value(True), 10)
         self.assertEqual(g.unit(), myokit.units.m)
         self.assertEqual(g.rhs().unit(), myokit.units.m / myokit.units.s)
         self.assertEqual(m.count_states(), 3)
@@ -158,6 +163,7 @@ def test_convert_unit(self):
         m = myokit.parse_model("""
             [[model]]
             membrane.V = -83
+            membrane.w = 10
 
             [env]
             t = 0 [ms] bind time
@@ -168,6 +174,8 @@ def test_convert_unit(self):
                 in [mV]
             dotv = 5 [mV/ms] + dot(V)
                 in [mV/ms]
+            dot(w) = 0.1 [mg/ms]
+                in [mg]
 
             [cell]
             Cm = 0.123 [uF]
@@ -194,7 +202,7 @@ def test_convert_unit(self):
         v.convert_unit('mV')
         self.assertEqual(code, m.code())
 
-        # Convert state
+        # Convert state (note: initial value is treated as PrefixMinus expr!)
         vdot = v.rhs().eval()
         self.assertNotEqual(v.unit(), myokit.units.V)
         v.convert_unit('V')
@@ -205,6 +213,28 @@ def test_convert_unit(self):
         self.assertEqual(vdot, v.rhs().eval())
         m.check_units(myokit.UNIT_STRICT)
 
+        # Check units of the initial value have been converted
+        self.assertEqual(v.initial_value(True), -0.083)
+        # Not sure if this one is desirable, but -83 is parsed as the
+        # expression -(83)
+        self.assertEqual(
+            v.initial_value(),
+            myokit.Multiply(
+                myokit.PrefixMinus(myokit.Number(83)),
+                myokit.Number(0.001, myokit.units.V / myokit.units.mV)))
+
+        # Convert state with number initial value
+        w = m.get('membrane.w')
+        wdot = w.rhs().eval()
+        w.convert_unit('ug')
+        self.assertEqual(w.unit(), myokit.units.g / 1e6)
+        self.assertEqual(wdot * 1000, w.rhs().eval())
+        m.check_units(myokit.UNIT_STRICT)
+        self.assertEqual(w.initial_value(True), 10000)
+        self.assertEqual(w.initial_value(), myokit.Number(10000))
+        #self.assertEqual(
+        #    w.initial_value(), myokit.Number(10000, myokit.units.g / 1e6))
+
         # Convert non-state
         i = m.get('ina.i')
         self.assertNotEqual(i.unit(), myokit.parse_unit('uA/uF'))
@@ -226,6 +256,49 @@ def test_convert_unit(self):
         self.assertEqual(vdot, v.rhs().eval())
         m.check_units(myokit.UNIT_STRICT)
 
+    def test_initial_value(self):
+        # Tests :meth:`Variable.initial_value`.
+
+        m = myokit.parse_model('''
+            [[model]]
+            c.x = 2
+
+            [c]
+            t = 0 bind time
+            dot(x) = 0
+            p = 3
+            q = 4
+        ''')
+
+        # Test expression version
+        x, p = m.get('c.x'), m.get('c.p')
+        p = m.get('c.p')
+        self.assertEqual(x.initial_value(), myokit.Number(2))
+        x.set_initial_value(3)
+        self.assertEqual(x.initial_value(), myokit.Number(3))
+        x.set_initial_value('3 [kg]')
+        self.assertEqual(x.initial_value(), myokit.Number(3, myokit.units.kg))
+        x.set_initial_value('sqrt(c.p) / log(c.p)')
+        self.assertEqual(x.initial_value(), myokit.Divide(
+            myokit.Sqrt(p.lhs()), myokit.Log(p.lhs())))
+
+        # Test float version
+        p.set_rhs(9)
+        x.set_initial_value('1 + sqrt(c.p)')
+        self.assertEqual(x.initial_value(True), 4)
+
+        # Test deprecated alias
+        with WarningCollector() as w:
+            self.assertEqual(x.state_value(), x.initial_value(True))
+        self.assertIn('deprecated', w.text())
+
+        # Test cycles are detected
+        p.set_rhs('1 - 2 * q')
+        m.get('c.q').set_rhs('7 * sqrt(p)')
+        x.initial_value()  # No evaluation, so no error
+        self.assertRaises(
+            myokit.CyclicalDependencyError, x.initial_value, as_float=True)
+
     def test_is_referenced(self):
         # Test :meth:`Variable.is_referenced().
 
@@ -276,6 +349,8 @@ def test_promote_demote(self):
 
         m = myokit.Model()
         c = m.add_component('c')
+        p = c.add_variable('p')
+        p.set_rhs(1)
         v = c.add_variable('v')
         v.set_rhs(3)
 
@@ -286,7 +361,7 @@ def test_promote_demote(self):
         self.assertEqual(v.lhs(), myokit.Name(v))
         self.assertRaises(Exception, v.demote)
         self.assertRaises(Exception, v.indice)
-        self.assertRaises(Exception, v.state_value)
+        self.assertRaises(Exception, v.initial_value)
 
         v.promote(3)
         self.assertFalse(v.is_literal())
@@ -295,7 +370,7 @@ def test_promote_demote(self):
         self.assertTrue(v.is_state())
         self.assertEqual(v.lhs(), myokit.Derivative(myokit.Name(v)))
         self.assertEqual(v.indice(), 0)
-        self.assertEqual(v.state_value(), 3)
+        self.assertEqual(v.initial_value(), myokit.Number(3))
 
         v.demote()
         self.assertTrue(v.is_literal())
@@ -305,7 +380,7 @@ def test_promote_demote(self):
         self.assertEqual(v.lhs(), myokit.Name(v))
         self.assertRaises(Exception, v.demote)
         self.assertRaises(Exception, v.indice)
-        self.assertRaises(Exception, v.state_value)
+        self.assertRaises(Exception, v.initial_value)
 
         # Test errors
         v.promote(3)
@@ -317,6 +392,19 @@ def test_promote_demote(self):
         self.assertRaisesRegex(
             Exception, 'only be added to Components', w.promote, 4)
 
+        # Test that promote accepts expressions and strings
+        v.set_binding(None)
+        v.remove_variable(w)
+        v.promote(myokit.Number(3, myokit.units.kg))
+        self.assertEqual(v.initial_value(), myokit.Number(3, myokit.units.kg))
+        v.demote()
+        v.promote('1 + c.p')
+        self.assertEqual(v.initial_value().code(), '1 + c.p')
+        v.demote()
+
+        # Test that promote only works with global syntax references
+        self.assertRaisesRegex(myokit.ParseError, 'Unresolved', v.promote, 'p')
+
         # Test we can't demote a variable with references to its derivative
         m = myokit.Model()
         c = m.add_component('c')
@@ -330,6 +418,15 @@ def test_promote_demote(self):
         y.set_rhs('1 + x')
         x.demote()
 
+        # Test deprecated keyword argument
+        with WarningCollector() as w:
+            x.promote(state_value=3)
+        self.assertIn('deprecated', w.text())
+        self.assertEqual(x.initial_value(), myokit.Number(3))
+        x.demote()
+        self.assertRaisesRegex(
+            Exception, 'at the same time as', x.promote, 3, state_value=3)
+
     def test_pyfunc(self):
         # Test :meth:`Variable.pyfunc().
 
@@ -638,31 +735,72 @@ def test_rename(self):
         self.assertEqual(v.name(), 'w')
         self.assertEqual(v.qname(), 'c.w')
 
-    def test_set_state_value(self):
-        # Test :meth:`Variable.set_state_value()`.
+    def test_set_initial_value(self):
+        # Test :meth:`Variable.set_initial_value()`.
 
         m = myokit.Model()
         c = m.add_component('c')
+        t = c.add_variable('t')
+        t.set_rhs(0)
+        t.set_binding('time')
         v = c.add_variable('v')
+        v.set_rhs(3)
         w = c.add_variable('w')
+        w.set_rhs('1 + 2')
+        m.validate()
 
         # Test basic functionality
         v.promote(10)
-        self.assertEqual(v.state_value(), 10)
-        v.set_state_value(12)
-        self.assertEqual(v.state_value(), 12)
+        self.assertEqual(v.initial_value(), myokit.Number(10))
+        v.set_initial_value(myokit.Number(12, myokit.units.g))
+        self.assertEqual(v.initial_value(), myokit.Number(12, myokit.units.g))
+
+        # Test setting expressions
+        v.set_initial_value('1 + 11')
+        #self.assertEqual(v.initial_value(),
+        #                 myokit.Plus(myokit.Number(1), myokit.Number(11)))
+        v.set_initial_value('1 + c.w')
+        self.assertEqual(v.initial_value(),
+                         myokit.Plus(myokit.Number(1), myokit.Name(w)))
+
+        # Strings must use global context
+        self.assertRaisesRegex(myokit.ParseError, 'No component specified for',
+                               v.set_initial_value, '1 / w')
 
         # Only states have this option
         v.demote()
         self.assertRaisesRegex(
-            Exception, 'Only state variables', v.set_state_value, 3)
+            Exception, 'Only state variables', v.set_initial_value, 3)
         self.assertRaisesRegex(
-            Exception, 'Only state variables', w.set_state_value, 3)
-
-        # State values must be literals
+            Exception, 'Only state variables', w.set_initial_value, 3)
         v.promote(3)
+        m.validate()
+
+        # Initial values must be constant, but this is only picked up in
+        # validation
+        w.promote(1)
+        v.set_initial_value(myokit.Name(w))
+        self.assertRaises(myokit.IntegrityError, m.validate)
+        w.demote()
+        m.validate()
+
+        # Nested variables may not be used, but this is only picked up when
+        # validating
+        x = v.add_variable('x')
+        x.set_rhs(1)
+        v.set_initial_value(x.lhs())
         self.assertRaises(
-            myokit.NonLiteralValueError, v.set_state_value, w.lhs())
+            myokit.IllegalReferenceInInitialValueError, m.validate)
+        v.set_initial_value(x)
+
+        # InitialValue and PartialDerivative can not be used
+        v.set_initial_value(myokit.PartialDerivative(v.lhs(), x.lhs()))
+
+        # Deprecated alias
+        with WarningCollector() as w:
+            v.set_state_value(1)
+        self.assertIn('deprecated', w.text())
+        self.assertEqual(v.initial_value(), myokit.Number(1))
 
     def test_set_unit(self):
         # Test :meth:`Variable.set_unit()`.
@@ -745,6 +883,59 @@ def test_validate(self):
         p.set_rhs('q == 3')
         self.assertRaisesRegex(
             myokit.IntegrityError, 'can not be a condition', p.validate)
+        p.set_rhs(3)
+        p.validate()
+
+        # Initial value can not reference nested variables
+        r = c.add_variable('r')
+        r1 = r.add_variable('r1')
+        r1.set_rhs(1)
+        r.set_rhs('1 + r1')
+        p.validate()
+        p.set_initial_value(myokit.Plus(myokit.Number(1), myokit.Name(r1)))
+        self.assertRaises(
+            myokit.IllegalReferenceInInitialValueError, p.validate)
+
+        # Initial value must be constant
+        p.set_initial_value('1 / c.q')
+        q.promote(3)
+        self.assertRaisesRegex(myokit.IntegrityError, 'not const', p.validate)
+        q.demote()
+        p.validate()
+        q.set_rhs('1 + c.p')
+        self.assertRaisesRegex(myokit.IntegrityError, 'not const', p.validate)
+        q.set_rhs(10)
+        p.validate()
+        q.set_binding('hello')
+        self.assertRaisesRegex(myokit.IntegrityError, 'not const', p.validate)
+        # Don't test everything, is_const() should work!
+
+        # Initial value cannot be a partial, init, or condition
+        p.set_initial_value(myokit.PartialDerivative(p.lhs(), q.lhs()))
+        self.assertRaisesRegex(
+            myokit.IntegrityError, 'Partial derivatives', p.validate)
+        p.set_initial_value(myokit.InitialValue(myokit.Name(p)))
+        self.assertRaisesRegex(
+            myokit.IntegrityError, 'Initial value', p.validate)
+        p.set_initial_value('c.q == 3')
+        self.assertRaisesRegex(
+            myokit.IntegrityError, 'can not be a condition', p.validate)
+
+        # Initial value cycles are caught in validation, don't cause loops
+        m = myokit.parse_model('''
+            [[model]]
+            x.x = 0
+            [x]
+            dot(x) = x
+            t = 0 bind time
+            ''')
+        x = m.get('x.x')
+        x.set_initial_value('x.x')
+        self.assertRaisesRegex(myokit.IntegrityError, 'not const', m.validate)
+        x.set_initial_value(1)
+        m.validate()
+        x.set_initial_value('1 + 1 / x.x')
+        self.assertRaisesRegex(myokit.IntegrityError, 'not const', m.validate)
 
     def test_value(self):
         # Test :meth:`Variable.value()`.