diff --git a/brian2/equations/codestrings.py b/brian2/equations/codestrings.py
index 39fc07a03..7e6f1644e 100644
--- a/brian2/equations/codestrings.py
+++ b/brian2/equations/codestrings.py
@@ -3,13 +3,18 @@
 information about its namespace. Only serves as a parent class, its subclasses
 `Expression` and `Statements` are the ones that are actually used.
 '''
-
+import re
+import string
 from collections.abc import Hashable
+from typing import Sequence
+import numbers
 
 import sympy
+import numpy as np
 
 from brian2.utils.logger import get_logger
 from brian2.utils.stringtools import get_identifiers
+from brian2.utils.topsort import topsort
 from brian2.parsing.sympytools import str_to_sympy, sympy_to_str
 
 __all__ = ['Expression', 'Statements']
@@ -35,6 +40,7 @@ def __init__(self, code):
 
         # : Set of identifiers in the code string
         self.identifiers = get_identifiers(code)
+        self.template_identifiers = get_identifiers(code, template=True)
 
     code = property(lambda self: self._code,
                     doc='The code string')
@@ -77,6 +83,11 @@ class Statements(CodeString):
     pass
 
 
+class Default(dict):
+    def __missing__(self, key):
+        return f'{{{key}}}'
+
+
 class Expression(CodeString):
     '''
     Class for representing an expression.
@@ -100,7 +111,7 @@ def __init__(self, code=None, sympy_expression=None):
 
         if code is None:
             code = sympy_to_str(sympy_expression)
-        else:
+        elif '{' not in code:
             # Just try to convert it to a sympy expression to get syntax errors
             # for incorrect expressions
             str_to_sympy(code)
@@ -193,6 +204,95 @@ def __ne__(self, other):
     def __hash__(self):
         return hash(self.code)
 
+    def _do_substitution(self, to_replace, replacement):
+        # Replacements can be lists, deal with single replacements
+        # as single-element lists
+        replaced_name = False
+        replaced_placeholder = False
+        if not isinstance(replacement, Sequence) or isinstance(replacement, str):
+            replacement = [replacement]
+        replacement_strs = []
+        for one_replacement in replacement:
+            if isinstance(one_replacement, str):
+                if any(c not in string.ascii_letters + '_{}'
+                       for c in one_replacement):
+                    # Check whether the replacement can be interpreted as an expression
+                    try:
+                        expr = Expression(one_replacement)
+                        replacement_strs.append(expr.code)
+                    except SyntaxError:
+                        raise SyntaxError(f'Replacement \'{one_replacement}\' for'
+                                          f'\'{to_replace}\' is neither a name nor a '
+                                          f'valid expression.')
+                else:
+                    replacement_strs.append(one_replacement)
+            elif isinstance(one_replacement, (numbers.Number, np.ndarray)):
+                if not getattr(one_replacement, 'shape', ()) == ():
+                    raise TypeError(f'Cannot replace variable \'{to_replace}\' with an '
+                                    f'array of values.')
+                replacement_strs.append(repr(one_replacement))
+            elif isinstance(one_replacement, Expression):
+                replacement_strs.append(one_replacement.code)
+            else:
+                raise TypeError(f'Cannot replace \'{to_replace}\' with an object of type '
+                                f'\'{type(one_replacement)}\'.')
+
+        if len(replacement_strs) == 1:
+            replacement_str = replacement_strs[0]
+            # Be careful if the string is more than just a name/number
+            if any(c not in string.ascii_letters + string.digits + '_.{}'
+                   for c in replacement_str):
+                replacement_str = '(' + replacement_str + ')'
+        else:
+            replacement_str = '(' + (' + '.join(replacement_strs)) + ')'
+
+        new_expr = self
+        if to_replace in new_expr.identifiers:
+            code = new_expr.code
+            new_expr = Expression(re.sub(r'(?<!\w|{)' + to_replace + r'(?!\w|})',
+                                         replacement_str, code))
+            replaced_name = True
+        if to_replace in new_expr.template_identifiers:
+            code = new_expr.code
+            new_expr = Expression(code.replace('{' + to_replace + '}',
+                                               replacement_str))
+            replaced_placeholder = True
+        if not (replaced_name or replaced_placeholder):
+            raise KeyError(f'Replacement argument \'{to_replace}\' does not correspond '
+                           f'to any name or placeholder in the equations.')
+        if replaced_name and replaced_placeholder:
+            logger.warn(f'Replacement argument \'{to_replace}\' replaced both a name '
+                        f'and a placeholder \'{{{to_replace}}}\'.',
+                        name_suffix='ambiguous_replacement')
+        return new_expr
+
+    def __call__(self, **replacements):
+        if len(replacements) == 0:
+            return self
+
+        # Figure out in which order elements should be substituted
+        dependencies = {}
+        for to_replace, replacement in replacements.items():
+            if not isinstance(replacement, Sequence) or isinstance(replacement, str):
+                replacement = [replacement]
+            for one_replacement in replacement:
+                dependencies[to_replace] = set()
+                if not isinstance(one_replacement, (numbers.Number, np.ndarray, str, Expression)):
+                    raise TypeError(f'Cannot use an object of type \'{type(one_replacement)}\''
+                                    f'to replace \'{to_replace}\' in an expression.')
+                if isinstance(one_replacement, Expression):
+                    dependencies[to_replace] |= one_replacement.identifiers | one_replacement.template_identifiers
+        # We only care about dependencies to values that are replaced at the same time
+        for dep_key, deps in dependencies.items():
+            dependencies[dep_key] = {d for d in deps if d in dependencies}
+
+        replacements_in_order = topsort(dependencies)[::-1]
+        expr = self
+        for to_replace in replacements_in_order:
+            replacement = replacements[to_replace]
+            expr = expr._do_substitution(to_replace, replacement)
+        return expr
+
 
 def is_constant_over_dt(expression, variables, dt_value):
     '''
diff --git a/brian2/equations/equations.py b/brian2/equations/equations.py
index aae2fe041..5a0c80c37 100644
--- a/brian2/equations/equations.py
+++ b/brian2/equations/equations.py
@@ -2,10 +2,12 @@
 '''
 Differential equations for Brian models.
 '''
+import numbers
 from collections.abc import Mapping, Hashable
 import keyword
 import re
 import string
+from typing import Sequence
 
 from pyparsing import (Group, ZeroOrMore, OneOrMore, Optional, Word, CharsNotIn,
                        Combine, Suppress, restOfLine, LineEnd, ParseException)
@@ -64,8 +66,8 @@
 # combination of letters, numbers and underscores
 # Note that the check_identifiers function later performs more checks, e.g.
 # names starting with underscore should only be used internally
-IDENTIFIER = Word(string.ascii_letters + '_',
-                  string.ascii_letters + string.digits + '_').setResultsName('identifier')
+IDENTIFIER = Word(string.ascii_letters + '_' + '{',
+                  string.ascii_letters + string.digits + '_' + '{' + '}').setResultsName('identifier')
 
 # very broad definition here, expression will be analysed by sympy anyway
 # allows for multi-line expressions, where each line can have comments
@@ -324,7 +326,7 @@ def dimensions_and_type_from_string(unit_string):
 
 
 @cached
-def parse_string_equations(eqns):
+def parse_string_equations(eqns, template=False):
     """
     parse_string_equations(eqns)
 
@@ -335,7 +337,9 @@ def parse_string_equations(eqns):
     eqns : str
         The (possibly multi-line) string defining the equations. See the
         documentation of the `Equations` class for details.
-    
+    template : bool
+        Whether the equations should be parsed as a template
+
     Returns
     -------
     equations : dict
@@ -363,7 +367,7 @@ def parse_string_equations(eqns):
                                 'variable "%s": %s' % (identifier, ex))
 
         expression = eq_content.get('expression', None)
-        if not expression is None:
+        if expression is not None:
             # Replace multiple whitespaces (arising from joining multiline
             # strings) with single space
             p = re.compile(r'\s{2,}')
@@ -441,10 +445,21 @@ def __init__(self, type, varname, dimensions, var_type=FLOAT, expr=None,
                            if not self.expr is None else set([]),
                            doc='All identifiers in the RHS of this equation.')
 
+    template_identifiers = property(lambda self: (self.expr.template_identifiers
+                                                  if not self.expr is None else set([])) |
+                                                 get_identifiers(self.varname, template=True),
+                                    doc='All template identifiers used in this equation (including as part of the variable '
+                                        'name.')
+
+    template = property(lambda self: len(self.template_identifiers) > 0)
+
     stochastic_variables = property(lambda self: {variable for variable in self.identifiers
                                                   if variable =='xi' or variable.startswith('xi_')},
                                     doc='Stochastic variables in the RHS of this equation')
 
+    template = property(lambda self: self.varname.startswith('{') or
+                                     (self.expr is not None and '{' in self.expr.code))
+
     def __eq__(self, other):
         if not isinstance(other, SingleEquation):
             return NotImplemented
@@ -555,7 +570,7 @@ class Equations(Hashable, Mapping):
 
     def __init__(self, eqns, **kwds):
         if isinstance(eqns, str):
-            self._equations = parse_string_equations(eqns)
+            self._equations = parse_string_equations(eqns, template=True)
             # Do a basic check for the identifiers
             self.check_identifiers()
         else:
@@ -569,6 +584,7 @@ def __init__(self, eqns, **kwds):
                                         eq.varname)
                 self._equations[eq.varname] = eq
 
+        self._orig_equations = self._equations
         self._equations = self._substitute(kwds)
 
         # Check for special symbol xi (stochastic term)
@@ -593,62 +609,149 @@ def __init__(self, eqns, **kwds):
         #: Cache for equations with the subexpressions substituted
         self._substituted_expressions = None
 
-    def _substitute(self, replacements):
-        if len(replacements) == 0:
-            return self._equations
-
+    def _do_substitution(self, equations, to_replace, replacement):
+        # Replacements can be lists, deal with single replacements
+        # as single-element lists
         new_equations = {}
-        for eq in self.values():
-            # Replace the name of a model variable (works only for strings)
-            if eq.varname in replacements:
-                new_varname = replacements[eq.varname]
-                if not isinstance(new_varname, str):
-                    raise ValueError(('Cannot replace model variable "%s" '
-                                      'with a value') % eq.varname)
-                if new_varname in self or new_varname in new_equations:
-                    raise EquationError(
-                        ('Cannot replace model variable "%s" '
-                         'with "%s", duplicate definition '
-                         'of "%s".' % (eq.varname, new_varname,
-                                       new_varname)))
-                # make sure that the replacement is a valid identifier
-                Equations.check_identifier(new_varname)
+        additional_equations = {}
+        replaced_name = False
+        replaced_placeholder = False
+        if not isinstance(replacement, Sequence) or isinstance(replacement, str):
+            replacement = [replacement]
+        replacement_strs = []
+        for one_replacement in replacement:
+            if isinstance(one_replacement, str):
+                if any(c not in string.ascii_letters + '_{}'
+                       for c in one_replacement):
+                    # Check whether the replacement can be interpreted as an expression
+                    try:
+                        expr = Expression(one_replacement)
+                        replacement_strs.append(expr.code)
+                    except SyntaxError:
+                        raise SyntaxError(f'Replacement \'{one_replacement}\' for'
+                                          f'\'{to_replace}\' is neither a name nor a '
+                                          f'valid expression.')
+                else:
+                    replacement_strs.append(one_replacement)
+            elif isinstance(one_replacement, (numbers.Number, Quantity)):
+                if not getattr(one_replacement, 'shape', ()) == ():
+                    raise TypeError(f'Cannot replace variable \'{to_replace}\' with an '
+                                    f'array of values.')
+                replacement_strs.append(repr(one_replacement))
+            elif isinstance(one_replacement, Expression):
+                replacement_strs.append(one_replacement.code)
+            elif isinstance(one_replacement, Equations):
+                replacement_strs.append(list(one_replacement.keys())[0])  # name of first equation
+                for additional_eq in one_replacement:
+                    if additional_eq in equations or additional_eq in additional_equations:
+                        raise SyntaxError(f'Adding equations to replace \'{to_replace}\' leads '
+                                          f'to duplicated definition of variable \'{additional_eq}\'')
+                    additional_equations[additional_eq] = one_replacement[additional_eq]
             else:
-                new_varname = eq.varname
-
-            if eq.type in [SUBEXPRESSION, DIFFERENTIAL_EQUATION]:
-                # Replace values in the RHS of the equation
-                new_code = eq.expr.code
-                for to_replace, replacement in replacements.items():
-                    if to_replace in eq.identifiers:
-                        if isinstance(replacement, str):
-                            # replace the name with another name
-                            new_code = re.sub('\\b' + to_replace + '\\b',
-                                              replacement, new_code)
-                        else:
-                            # replace the name with a value
-                            new_code = re.sub('\\b' + to_replace + '\\b',
-                                              '(' + repr(replacement) + ')',
-                                              new_code)
-                        try:
-                            Expression(new_code)
-                        except ValueError as ex:
-                            raise ValueError(
-                                ('Replacing "%s" with "%r" failed: %s') %
-                                (to_replace, replacement, ex))
-                new_equations[new_varname] = SingleEquation(eq.type, new_varname,
-                                                            dimensions=eq.dim,
-                                                            var_type=eq.var_type,
-                                                            expr=Expression(new_code),
-                                                            flags=eq.flags)
+                raise TypeError(f'Cannot replace \'{to_replace}\' with an object of type '
+                                f'\'{type(one_replacement)}\'.')
+        if len(replacement_strs) == 1:
+            replacement_str = replacement_strs[0]
+            # Be careful if the string is more than just a name/number
+            if any(c not in string.ascii_letters + string.digits + '_.{}'
+                   for c in replacement_str):
+                replacement_str = '(' + replacement_str + ')'
+        else:
+            replacement_str = '(' + (' + '.join(replacement_strs)) + ')'
+
+        for eq in equations.values():
+            # Check whether the variable name itself (or part of it) will be replaced
+            if eq.varname == to_replace:
+                if not len(replacement) == 1:
+                    raise TypeError(f'Cannot replace variable name \'{to_replace}\' with '
+                                    f'a list of values.')
+                if not isinstance(replacement[0], str):
+                    raise TypeError(f'Cannot replace variable name \'{to_replace}\' with '
+                                    f'an object of type \'{type(replacement[0])}\'.')
+                new_varname = replacement[0]
+                replaced_name = True
+            elif '{' + to_replace + '}' in eq.varname:
+                if not len(replacement) == 1:
+                    raise TypeError(f'Cannot replace \'{{{to_replace}}}\' as a part of a variable'
+                                    f'name with a list of values.')
+                if not isinstance(replacement[0], str):
+                    raise TypeError(f'Cannot replace \'{{{to_replace}}}\' as a part of a variable'
+                                    f'name with an object of type \'{type(replacement[0])}\'.')
+                new_varname = eq.varname.replace('{' + to_replace + '}', replacement[0])
+                replaced_placeholder = True
             else:
-                new_equations[new_varname] = SingleEquation(eq.type, new_varname,
-                                                            dimensions=eq.dim,
-                                                            var_type=eq.var_type,
-                                                            flags=eq.flags)
+                new_varname = eq.varname
 
+            # Check whether the new variable name is still valid
+            if '{' not in new_varname:
+                Equations.check_identifier(new_varname)
+                if new_varname != eq.varname:
+                    if new_varname in equations:
+                        raise EquationError(f'Cannot replace \'{eq.varname}\' by \'{new_varname}\', ' \
+                                            'this name is already used for another variable.')
+            # Replace occurrences in the RHS of equations
+            new_expr = eq.expr
+            if to_replace in eq.identifiers and eq.expr is not None:
+                code = eq.expr.code
+                new_expr = Expression(re.sub(r'(?<!\w|{)' + to_replace + r'(?!\w|})',
+                                             replacement_str, code))
+                replaced_name = True
+            if to_replace in eq.template_identifiers and new_expr is not None:
+                code = new_expr.code
+                new_expr = Expression(code.replace('{'+to_replace+'}',
+                                                   replacement_str))
+                replaced_placeholder = True
+
+            new_equations[new_varname] = SingleEquation(eq.type,
+                                                        new_varname,
+                                                        dimensions=eq.dim,
+                                                        var_type=eq.var_type,
+                                                        expr=new_expr,
+                                                        flags=eq.flags)
+        if not (replaced_name or replaced_placeholder):
+            raise KeyError(f'Replacement argument \'{to_replace}\' does not correspond '
+                           f'to any name or placeholder in the equations.')
+        if replaced_name and replaced_placeholder:
+            logger.warn(f'Replacement argument \'{to_replace}\' replaced both a name '
+                        f'and a placeholder \'{{{to_replace}}}\'.',
+                        name_suffix='ambiguous_replacement')
+        new_equations.update(additional_equations)
         return new_equations
 
+    def _substitute(self, replacements):
+        if len(replacements) == 0:
+            return self._equations
+
+        # Figure out in which order elements should be substituted
+        dependencies = {}
+        for to_replace, replacement in replacements.items():
+            if not isinstance(replacement, Sequence) or isinstance(replacement, str):
+                replacement = [replacement]
+            for one_replacement in replacement:
+                dependencies[to_replace] = set()
+                if not isinstance(one_replacement, (numbers.Number, Quantity, str, Expression, Equations)):
+                    raise TypeError(f'Cannot use an object of type \'{type(one_replacement)}\''
+                                    f'to replace \'{to_replace}\'')
+                if isinstance(one_replacement, (Expression, Equations)):
+                    dependencies[to_replace] |= one_replacement.identifiers | one_replacement.template_identifiers
+        # We only care about dependencies to values that are replaced at the same time
+        for dep_key, deps in dependencies.items():
+            dependencies[dep_key] = {d for d in deps if d in dependencies}
+
+        replacements_in_order = topsort(dependencies)[::-1]
+
+        new_equations = dict(self)
+        # The replacement code below generates a full set of equations for each replacement
+        # which is unnecessary most of the time, but simplify things a lot
+        for to_replace in replacements_in_order:
+            replacement = replacements[to_replace]
+            new_equations = self._do_substitution(new_equations, to_replace, replacement)
+
+        return Equations([eq for eq in new_equations.values()])
+
+    def __call__(self, **kwds):
+        return Equations(list(self._substitute(kwds).values()))
+
     def substitute(self, **kwds):
         return Equations(list(self._substitute(kwds).values()))
 
@@ -879,6 +982,11 @@ def _get_stochastic_type(self):
                            doc=('Set of all identifiers used in the equations, '
                                 'excluding the variables defined in the equations'))
 
+    template_identifiers = property(lambda self: set().union(*[eq.template_identifiers for
+                                                               eq in self._equations.values()]),
+                                    doc=('Set of all template identifiers (placeholders) used in the'
+                                         'equations, including as part of the variable names.'))
+
     stochastic_variables = property(lambda self: {variable for variable in self.identifiers
                                                   if variable =='xi' or variable.startswith('xi_')})
 
diff --git a/brian2/groups/neurongroup.py b/brian2/groups/neurongroup.py
index a48e4bb3a..3a3b692f8 100644
--- a/brian2/groups/neurongroup.py
+++ b/brian2/groups/neurongroup.py
@@ -503,6 +503,13 @@ def __init__(self, N, model,
         if not isinstance(model, Equations):
             raise TypeError(('model has to be a string or an Equations '
                              'object, is "%s" instead.') % type(model))
+        if len(model.template_identifiers) > 0:
+            identifiers = [f'\'{identifier}\''
+                           for identifier in model.template_identifiers]
+            identifier_list = ', '.join(sorted(identifiers))
+            raise TypeError('The model equations contain placeholders, substitute '
+                            'names/values for the ' f'following before passing '
+                            f'them to {self.__class__.__name__}: {identifier_list}.')
 
         # Check flags
         model.check_flags({DIFFERENTIAL_EQUATION: ('unless refractory',),
diff --git a/brian2/only.py b/brian2/only.py
index f73167cc3..cf374d76a 100644
--- a/brian2/only.py
+++ b/brian2/only.py
@@ -81,7 +81,8 @@ def restore_initial_state():
     'DEFAULT_FUNCTIONS', 'Function', 'implementation', 'declare_types',
     'PreferenceError', 'BrianPreference', 'prefs', 'brian_prefs',
     'Clock', 'defaultclock',
-    'Equations', 'Expression', 'Statements',
+    'Equations', 'Expression',
+    'Statements',
     'BrianObject',
     'BrianObjectException',
     'Network', 'profiling_summary', 'scheduling_summary',
diff --git a/brian2/synapses/synapses.py b/brian2/synapses/synapses.py
index 17f507b32..b57e8a020 100644
--- a/brian2/synapses/synapses.py
+++ b/brian2/synapses/synapses.py
@@ -740,6 +740,14 @@ def __init__(self, source, target=None, model=None, on_pre=None,
             raise TypeError(('model has to be a string or an Equations '
                              'object, is "%s" instead.') % type(model))
 
+        if len(model.template_identifiers) > 0:
+            identifiers = [f'\'{identifier}\''
+                           for identifier in model.template_identifiers]
+            identifier_list = ', '.join(sorted(identifiers))
+            raise TypeError('The model equations contain placeholders, substitute '
+                            'names/values for the ' f'following before passing '
+                            f'them to {self.__class__.__name__}: {identifier_list}.')
+
         # Check flags
         model.check_flags({DIFFERENTIAL_EQUATION: ['event-driven', 'clock-driven'],
                            SUBEXPRESSION: ['summed', 'shared',
diff --git a/brian2/tests/test_equations.py b/brian2/tests/test_equations.py
index 0c614750c..814a806b4 100644
--- a/brian2/tests/test_equations.py
+++ b/brian2/tests/test_equations.py
@@ -221,13 +221,17 @@ def test_wrong_replacements():
                   ''', v='x')
 
     # Replacing a model variable name with a value
-    with pytest.raises(ValueError):
+    with pytest.raises(TypeError):
         Equations('dv/dt = -v / tau : 1', v=3 * mV)
 
     # Replacing with an illegal value
-    with pytest.raises(SyntaxError):
+    with pytest.raises(TypeError):
         Equations('dv/dt = -v/tau : 1', tau=np.arange(5))
 
+    # Replacing something that does not exist
+    with pytest.raises(KeyError):
+        Equations('dv/dt = -v/tau : 1', x='y')
+
 
 @pytest.mark.codegen_independent
 def test_substitute():
diff --git a/brian2/utils/stringtools.py b/brian2/utils/stringtools.py
index c37239280..ccf17922f 100644
--- a/brian2/utils/stringtools.py
+++ b/brian2/utils/stringtools.py
@@ -152,11 +152,12 @@ def replace(s, substitutions):
 KEYWORDS = {'and', 'or', 'not', 'True', 'False'}
 
 
-def get_identifiers(expr, include_numbers=False):
+def get_identifiers(expr, include_numbers=False, template=False):
     '''
     Return all the identifiers in a given string ``expr``, that is everything
-    that matches a programming language variable like expression, which is
-    here implemented as the regexp ``\\b[A-Za-z_][A-Za-z0-9_]*\\b``.
+    that matches a programming language variable like expression. Placeholder
+    arguments of the form ``{name}`` are handled separately and only returned
+    if the ``template`` argument is ``True``.
 
     Parameters
     ----------
@@ -164,7 +165,9 @@ def get_identifiers(expr, include_numbers=False):
         The string to analyze
     include_numbers : bool, optional
         Whether to include number literals in the output. Defaults to ``False``.
-
+    template : bool, optional
+        Whether to only return template identifiers, i.e. identifiers enclosed
+        by curly braces. Defaults to ``False``
     Returns
     -------
     identifiers : set
@@ -174,20 +177,30 @@ def get_identifiers(expr, include_numbers=False):
     --------
     >>> expr = '3-a*_b+c5+8+f(A - .3e-10, tau_2)*17'
     >>> ids = get_identifiers(expr)
-    >>> print(sorted(list(ids)))
+    >>> print(sorted(ids))
     ['A', '_b', 'a', 'c5', 'f', 'tau_2']
     >>> ids = get_identifiers(expr, include_numbers=True)
-    >>> print(sorted(list(ids)))
+    >>> print(sorted(ids))
     ['.3e-10', '17', '3', '8', 'A', '_b', 'a', 'c5', 'f', 'tau_2']
+    >>> template_expr = '{name}_{suffix} = a*{name}_{suffix} + b'
+    >>> template_ids = get_identifiers(template_expr, template=True)
+    >>> print(sorted(template_ids))
+    ['name', 'suffix']
     '''
-    identifiers = set(re.findall(r'\b[A-Za-z_][A-Za-z0-9_]*\b', expr))
-    if include_numbers:
-        # only the number, not a + or -
-        numbers = set(re.findall(r'(?<=[^A-Za-z_])[0-9]*\.?[0-9]+(?:[eE][-+]?[0-9]+)?|^[0-9]*\.?[0-9]+(?:[eE][-+]?[0-9]+)?',
-                                 expr))
+    if include_numbers and template:
+        raise ValueError('Cannot combine the \'template\' and \'include_numbers\' arguments.')
+
+    if template:
+        return set(re.findall(r'(?:[{])([A-Za-z_][A-Za-z0-9_]*)(?:[}])', expr))
     else:
-        numbers = set()
-    return (identifiers - KEYWORDS) | numbers
+        if include_numbers:
+            # only the number, not a + or -
+            number_regexp = r'(?<=[^A-Za-z_])[0-9]*\.?[0-9]+(?:[eE][-+]?[0-9]+)?|^[0-9]*\.?[0-9]+(?:[eE][-+]?[0-9]+)?'
+            numbers = set(re.findall(number_regexp, expr))
+        else:
+            numbers = set()
+        identifiers = set(re.findall(r'(?<!\w|{)[A-Za-z_][A-Za-z0-9_]*(?!\w|})', expr))
+        return (identifiers - KEYWORDS) | numbers
 
 
 def strip_empty_lines(s):
diff --git a/examples/IF_curve_LIF.py b/examples/IF_curve_LIF.py
index 292579eb4..43a355999 100644
--- a/examples/IF_curve_LIF.py
+++ b/examples/IF_curve_LIF.py
@@ -7,7 +7,8 @@
 The input is set differently for each neuron.
 '''
 from brian2 import *
-
+import brian2tools
+set_device('exporter', debug=True)
 n = 1000
 duration = 1*second
 tau = 10*ms
@@ -23,7 +24,5 @@
 monitor = SpikeMonitor(group)
 
 run(duration)
-plot(group.v0/mV, monitor.count / duration)
-xlabel('v0 (mV)')
-ylabel('Firing rate (sp/s)')
-show()
+
+
diff --git a/examples/frompapers/Rothman_Manis_2003_with_templates.py b/examples/frompapers/Rothman_Manis_2003_with_templates.py
new file mode 100755
index 000000000..6b2f8ec64
--- /dev/null
+++ b/examples/frompapers/Rothman_Manis_2003_with_templates.py
@@ -0,0 +1,157 @@
+#!/usr/bin/env python
+"""
+Cochlear neuron model of Rothman & Manis
+----------------------------------------
+Rothman JS, Manis PB (2003) The roles potassium currents play in
+regulating the electrical activity of ventral cochlear nucleus neurons.
+J Neurophysiol 89:3097-113.
+
+All model types differ only by the maximal conductances.
+
+Adapted from their Neuron implementation by Romain Brette
+"""
+from brian2 import *
+
+#defaultclock.dt=0.025*ms # for better precision
+
+'''
+Simulation parameters: choose current amplitude and neuron type
+(from type1c, type1t, type12, type 21, type2, type2o)
+'''
+neuron_type = 'type1c'
+Ipulse = 250*pA
+
+C = 12*pF
+Eh = -43*mV
+EK = -70*mV  # -77*mV in mod file
+El = -65*mV
+ENa = 50*mV
+nf = 0.85  # proportion of n vs p kinetics
+zss = 0.5  # steady state inactivation of glt
+temp = 22.  # temperature in degree celcius
+q10 = 3. ** ((temp - 22) / 10.)
+# hcno current (octopus cell)
+frac = 0.0
+qt = 4.5 ** ((temp - 33.) / 10.)
+
+# Maximal conductances of different cell types in nS
+maximal_conductances = dict(
+type1c=(1000, 150, 0, 0, 0.5, 0, 2),
+type1t=(1000, 80, 0, 65, 0.5, 0, 2),
+type12=(1000, 150, 20, 0, 2, 0, 2),
+type21=(1000, 150, 35, 0, 3.5, 0, 2),
+type2=(1000, 150, 200, 0, 20, 0, 2),
+type2o=(1000, 150, 600, 0, 0, 40, 2) # octopus cell
+)
+gnabar, gkhtbar, gkltbar, gkabar, ghbar, gbarno, gl = [x * nS for x in maximal_conductances[neuron_type]]
+
+gating_var = Equations('''d{name}/dt = q10*({name}__inf - {name})/tau_{name} : 1
+                          {name}__inf = {rate_expression}                    : 1
+                          tau_{name} = {tau_scale}/({forward_rate} + 
+                                                    {reverse_rate})
+                                       + {tau_base}                          : second''')
+
+pos_sigmoid = Expression('1./(1+exp(-({voltage} - {midpoint}) / {scale}))')
+power_sigmoid = Expression('1./(1+exp(-({voltage} - {midpoint}) / {scale}))**{power}')
+neg_sigmoid = Expression('1./(1+exp(({voltage} - {midpoint}) / {scale}))')
+neg_power_sigmoid = Expression('1./(1+exp(({voltage} - {midpoint}) / {scale}))**{power}')
+shifted_neg_sigmoid = Expression('{shift} + (1.0-{shift})/(1+exp(({voltage} - {midpoint}) / {scale}))')
+exp_voltage_dep = Expression('{magnitude}*exp(({voltage}-{midpoint})/{scale})')
+neg_exp_voltage_dep = Expression('{magnitude}*exp(-({voltage}-{midpoint})/{scale})')
+
+# Classical Na channel
+ina = Equations('ina = gnabar*{m}**3*{h}*(ENa-v) : amp',
+                m=gating_var(name='m',
+                             rate_expression=pos_sigmoid(midpoint=-38., scale=7.),
+                             forward_rate=exp_voltage_dep(magnitude=5., midpoint=-60, scale=18.),
+                             reverse_rate=neg_exp_voltage_dep(magnitude=36., midpoint=-60, scale=25.),
+                             tau_base=0.04*ms, tau_scale=10*ms),
+                h=gating_var(name='h',
+                             rate_expression=neg_sigmoid(midpoint=-65., scale=6.),
+                             forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60., scale=11.),
+                             reverse_rate=neg_exp_voltage_dep(magnitude=10., midpoint=-60., scale=25.),
+                             tau_base=0.6*ms, tau_scale=100*ms))
+
+# KHT channel (delayed-rectifier K+)
+ikht = Equations('ikht = gkhtbar*(nf*{n}**2 + (1-nf)*{p})*(EK-v) : amp',
+                 n=gating_var(name='n',
+                              rate_expression=power_sigmoid(midpoint=-15, scale=5., power=0.5),
+                              forward_rate=exp_voltage_dep(magnitude=11., midpoint=-60, scale=24.),
+                              reverse_rate=neg_exp_voltage_dep(magnitude=21., midpoint=-60, scale=23.),
+                              tau_base=0.7*ms, tau_scale=100*ms),
+                 p=gating_var(name='p',
+                              rate_expression=pos_sigmoid(midpoint=-23., scale=6.),
+                              forward_rate=exp_voltage_dep(magnitude=4., midpoint=-60., scale=32.),
+                              reverse_rate=neg_exp_voltage_dep(magnitude=5., midpoint=-60., scale=22.),
+                              tau_base=5*ms, tau_scale=100*ms))
+
+# Ih channel (subthreshold adaptive, non-inactivating)
+ih = Equations('ih = ghbar*{r}*(Eh - v) : amp',
+               r=gating_var(name='r',
+                            rate_expression=neg_sigmoid(midpoint=-76., scale=7.),
+                            forward_rate=exp_voltage_dep(magnitude=237., midpoint=-60., scale=12.),
+                            reverse_rate=neg_exp_voltage_dep(magnitude=17, midpoint=-60, scale=14.),
+                            tau_scale=100*second, tau_base=25*ms))
+
+# KLT channel (low threshold K+)
+iklt = Equations('iklt = gkltbar*{w}**4*{z}*(EK-v) : amp',
+                 w=gating_var(name='w',
+                              rate_expression=power_sigmoid(midpoint=-48., scale=6., power=0.25),
+                              forward_rate=exp_voltage_dep(magnitude=6., midpoint=-60., scale=6.),
+                              reverse_rate=neg_exp_voltage_dep(magnitude=16, midpoint=-60, scale=45.),
+                              tau_scale=100*ms, tau_base=1.5*ms),
+                 z=gating_var(name='z',
+                              rate_expression=shifted_neg_sigmoid(shift=zss, midpoint=-71., scale=10.),
+                              forward_rate=exp_voltage_dep(magnitude=1., midpoint=-60., scale=20.),
+                              reverse_rate=neg_exp_voltage_dep(magnitude=1., midpoint=-60., scale=8.),
+                              tau_scale=1*second, tau_base=50*ms))
+# Ka channel (transient K+)
+ika = Equations('ika = gkabar*{a}**4*{b}*{c}*(EK-v): amp',
+                a=gating_var(name='a',
+                             rate_expression=power_sigmoid(midpoint=-31., scale=6., power=0.25),
+                             forward_rate=exp_voltage_dep(magnitude=7., midpoint=-60, scale=14.),
+                             reverse_rate=neg_exp_voltage_dep(magnitude=29., midpoint=-60, scale=24.),
+                             tau_scale=100*ms, tau_base=0.1*ms),
+                b=gating_var(name='b',
+                             rate_expression=neg_power_sigmoid(midpoint=-66., scale=7., power=0.5),
+                             forward_rate=exp_voltage_dep(magnitude=14., midpoint=-60., scale=27.),
+                             reverse_rate=neg_exp_voltage_dep(magnitude=29., midpoint=-60., scale=24.),
+                             tau_scale=1*second, tau_base=1*ms),
+                c=gating_var(name='c',
+                             rate_expression=neg_power_sigmoid(midpoint=-66., scale=7., power=0.5),
+                             forward_rate='1',
+                             reverse_rate=neg_exp_voltage_dep(magnitude=1., midpoint=-66., scale=17.),
+                             tau_scale=90*ms, tau_base=10*ms))
+# Leak
+ileak = Equations("ileak = gl*(El-v) : amp")
+
+# h current for octopus cells (gating variables use different equations
+gating_var_octopus = Equations('''d{name}/dt = (hinfno - {name})/tau_{name} : 1
+                                  tau_{name} = beta_{name}/(qt*{tau_scale}*(1+alpha_{name}))*ms : second
+                                  alpha_{name} = exp(1e-3*3*({voltage}-{midpoint})*9.648e4/(8.315*(273.16+temp))) : 1
+                                  beta_{name} = exp(1e-3*3*0.3*({voltage}-{midpoint})*9.648e4/(8.315*(273.16+temp))) : 1 ''')
+
+ihcno = Equations('''ihcno = gbarno*({h1}*frac + {h2}*(1-frac))*(Eh-v) : amp
+                     hinfno = {inf_rate} : 1''',
+                  inf_rate=neg_sigmoid(midpoint=-66., scale=7.),
+                  h1=gating_var_octopus(name='h1', tau_scale=0.008, midpoint=-50.),
+                  h2=gating_var_octopus(name='h2', tau_scale=0.0029, midpoint=-84.))
+
+eqs =Equations("""dv/dt = ({currents} + I)/C : volt
+                  I : amp""",
+               currents=[ileak, ina, ikht, ih, iklt, ika, ihcno],
+               voltage='v/mV')
+neuron = NeuronGroup(1, eqs, method='exponential_euler')
+neuron.v = El
+
+run(50*ms, report='text')  # Go to rest
+
+M = StateMonitor(neuron, 'v', record=0)
+neuron.I = Ipulse
+
+run(100*ms, report='text')
+
+plot(M.t / ms, M[0].v / mV)
+xlabel('t (ms)')
+ylabel('v (mV)')
+show()