diff --git a/bluepyopt/ephys/create_hoc.py b/bluepyopt/ephys/create_hoc.py
index 1f0ae4d..f85df55 100644
--- a/bluepyopt/ephys/create_hoc.py
+++ b/bluepyopt/ephys/create_hoc.py
@@ -76,7 +76,9 @@ def _generate_parameters(parameters):
             continue
         for param in param_locations[loc]:
             if isinstance(param, NrnRangeParameter):
-                if isinstance(param.value_scaler, NrnSegmentSomaDistanceScaler):
+                if isinstance(
+                        param.value_scaler,
+                        NrnSegmentSomaDistanceScaler):
                     value = param.value_scaler.inst_distribution
                     value = re.sub(r'math\.', '', value)
                     value = re.sub('{distance}', FLOAT_FORMAT, value)
@@ -97,27 +99,41 @@ def _generate_parameters(parameters):
     return global_params, ordered_section_params, range_params
 
 
-def create_hoc(mechs, parameters, morphology=None, ignored_globals=(),
-               replace_axon=None, template_name='CCell',
-               template='cell_template.jinja2', disable_banner=None):
+def create_hoc(
+        mechs,
+        parameters,
+        morphology=None,
+        ignored_globals=(),
+        replace_axon=None,
+        template_name='CCell',
+        template_filename='cell_template.jinja2',
+        disable_banner=None,
+        template_dir=None):
     '''return a string containing the hoc template
 
     Args:
-        mechs(): All the mechs for the hoc template
-        parameters(): All the parameters in the hoc template
-        morpholgy(str): Name of morphology
-        ignored_globals(iterable str): HOC coded is added for each
+        mechs (): All the mechs for the hoc template
+        parameters (): All the parameters in the hoc template
+        morpholgy (str): Name of morphology
+        ignored_globals (iterable str): HOC coded is added for each
         NrnGlobalParameter
         that exists, to test that it matches the values set in the parameters.
         This iterable contains parameter names that aren't checked
-        replace_axon(str): String replacement for the 'replace_axon' command.
+        replace_axon (str): String replacement for the 'replace_axon' command.
         Must include 'proc replace_axon(){ ... }
-        template(str): name of the template to use 'cell_template.jinja2',
+        template (str): file name of the jinja2 template
+        template_dir (str): dir name of the jinja2 template
     '''
-    templates_basepath = os.path.abspath(os.path.dirname(__file__))
-    template = os.path.join(templates_basepath, 'templates', template)
-    with open(template) as fd:
-        template = fd.read()
+
+    if template_dir is None:
+        template_dir = os.path.abspath(
+            os.path.join(
+                os.path.dirname(__file__),
+                'templates'))
+
+    template_path = os.path.join(template_dir, template_filename)
+    with open(template_path) as template_file:
+        template = template_file.read()
         template = jinja2.Template(template)
 
     channels = _generate_channels_by_location(mechs)
diff --git a/bluepyopt/ephys/models.py b/bluepyopt/ephys/models.py
index f31069f..e7c4389 100644
--- a/bluepyopt/ephys/models.py
+++ b/bluepyopt/ephys/models.py
@@ -260,7 +260,8 @@ def check_nonfrozen_params(self, param_names):  # pylint: disable=W0613
 
     def create_hoc(self, param_values,
                    ignored_globals=(), template='cell_template.jinja2',
-                   disable_banner=False):
+                   disable_banner=False,
+                   template_dir=None):
         """Create hoc code for this model"""
 
         to_unfreeze = []
@@ -282,7 +283,8 @@ def create_hoc(self, param_values,
                                     ignored_globals=ignored_globals,
                                     replace_axon=replace_axon,
                                     template_name=template_name,
-                                    template=template,
+                                    template_filename=template,
+                                    template_dir=template_dir,
                                     disable_banner=disable_banner)
 
         self.unfreeze(to_unfreeze)
diff --git a/bluepyopt/tests/test_ephys/test_create_hoc.py b/bluepyopt/tests/test_ephys/test_create_hoc.py
index b4e83f2..53cd63f 100644
--- a/bluepyopt/tests/test_ephys/test_create_hoc.py
+++ b/bluepyopt/tests/test_ephys/test_create_hoc.py
@@ -2,6 +2,8 @@
 
 # pylint: disable=W0212
 
+import os
+
 import utils
 from bluepyopt.ephys import create_hoc
 
@@ -49,3 +51,23 @@ def test_create_hoc():
     nt.ok_('CCell' in hoc)
     nt.ok_('begintemplate' in hoc)
     nt.ok_('endtemplate' in hoc)
+
+
+@attr('unit')
+def test_create_hoc_filename():
+    """ephys.create_hoc: Test create_hoc template_filename"""
+    mech = utils.make_mech()
+    parameters = utils.make_parameters()
+
+    hoc = create_hoc.create_hoc([mech, ],
+                                parameters, template_name='CCell',
+                                template_filename='test.jinja2',
+                                template_dir=os.path.join(
+                                    os.path.dirname(__file__),
+                                    'testdata'))
+    nt.ok_('load_file' in hoc)
+    nt.ok_('CCell' in hoc)
+    nt.ok_('begintemplate' in hoc)
+    nt.ok_('endtemplate' in hoc)
+    nt.ok_('Test template' in hoc)
+
diff --git a/bluepyopt/tests/test_ephys/testdata/test.jinja2 b/bluepyopt/tests/test_ephys/testdata/test.jinja2
new file mode 100644
index 0000000..f26cd58
--- /dev/null
+++ b/bluepyopt/tests/test_ephys/testdata/test.jinja2
@@ -0,0 +1,206 @@
+/* 
+Test template
+
+{%- if banner %}
+{{banner}} 
+{%- endif %}
+*/
+{load_file("stdrun.hoc")}
+{load_file("import3d.hoc")}
+
+{%- if global_params %}
+/*
+ * Check that global parameters are the same as with the optimization
+ */
+proc check_parameter(/* name, expected_value, value */){
+  strdef error
+  if($2 != $3){
+    sprint(error, "Parameter %s has different value %f != %f", $s1, $2, $3)
+    execerror(error)
+  }
+}
+proc check_simulator() {
+  {%- for param, value in global_params.items() %}
+  check_parameter("{{param}}", {{value}}, {{param}})
+  {%- endfor %}
+}
+{%- endif %}
+{%- if ignored_global_params %}
+/* The following global parameters were set in BluePyOpt
+{%- for param, value in ignored_global_params.items() %}
+ * {{param}} = {{value}}
+{%- endfor %}
+ */
+{%- endif %}
+
+begintemplate {{template_name}}
+  public init, morphology, geom_nseg_fixed, geom_nsec
+  public soma, dend, apic, axon, myelin
+  create soma[1], dend[1], apic[1], axon[1], myelin[1]
+
+  objref this, CellRef, segCounts
+
+  public all, somatic, apical, axonal, basal, myelinated, APC
+  objref all, somatic, apical, axonal, basal, myelinated, APC
+
+proc init(/* args: morphology_dir, morphology_name */) {
+  all = new SectionList()
+  apical = new SectionList()
+  axonal = new SectionList()
+  basal = new SectionList()
+  somatic = new SectionList()
+  myelinated = new SectionList()
+
+  //For compatibility with BBP CCells
+  CellRef = this
+
+  forall delete_section()
+
+  if(numarg() >= 2) {
+    load_morphology($s1, $s2)
+  } else {
+  {%- if morphology %}
+    load_morphology($s1, "{{morphology}}")
+  {%- else %}
+    execerror("Template {{template_name}} requires morphology name to instantiate")
+  {%- endif %}
+  }
+
+  geom_nseg()
+  {%- if replace_axon %}
+    replace_axon()
+  {%- endif %}
+  insertChannel()
+  biophys()
+  re_init_rng()
+}
+
+proc load_morphology(/* morphology_dir, morphology_name */) {localobj morph, import, sf, extension
+  strdef morph_path
+  sprint(morph_path, "%s/%s", $s1, $s2)
+
+  sf = new StringFunctions()
+  extension = new String()
+
+  sscanf(morph_path, "%s", extension.s)
+  sf.right(extension.s, sf.len(extension.s)-4)
+
+  if( strcmp(extension.s, ".asc") == 0 ) {
+    morph = new Import3d_Neurolucida3()
+  } else if( strcmp(extension.s, ".swc" ) == 0) {
+    morph = new Import3d_SWC_read()
+  } else {
+    printf("Unsupported file format: Morphology file has to end with .asc or .swc" )
+    quit()
+  }
+
+  morph.quiet = 1
+  morph.input(morph_path)
+
+  import = new Import3d_GUI(morph, 0)
+  import.instantiate(this)
+}
+
+/*
+ * Assignment of mechanism values based on distance from the soma
+ * Matches the BluePyOpt method
+ */
+proc distribute_distance(){local x localobj sl
+  strdef stmp, distfunc, mech
+
+  sl = $o1
+  mech = $s2
+  distfunc = $s3
+  this.soma[0] distance(0, 0.5)
+  sprint(distfunc, "%%s %s(%%f) = %s", mech, distfunc)
+  forsec sl for(x, 0) {
+    sprint(stmp, distfunc, secname(), x, distance(x))
+    execute(stmp)
+  }
+}
+
+proc geom_nseg() {
+  this.geom_nsec() //To count all sections
+  //TODO: geom_nseg_fixed depends on segCounts which is calculated by
+  //  geom_nsec.  Can this be collapsed?
+  this.geom_nseg_fixed(40)
+  this.geom_nsec() //To count all sections
+}
+
+proc insertChannel() {
+  {%- for location, names in channels.items() %}
+  forsec this.{{location}} {
+  {%- for channel in names %}
+    insert {{channel}}
+  {%- endfor %}
+  }
+  {%- endfor %}
+}
+
+proc biophys() {
+  {% for loc, parameters in section_params %}
+  forsec CellRef.{{ loc }} {
+  {%- for param in parameters %}
+    {{ param.name }} = {{ param.value }}
+  {%- endfor %}
+  }
+  {% endfor %}
+  {%- for location, param_name, value in range_params %}
+  distribute_distance(CellRef.{{location}}, "{{param_name}}", "{{value}}")
+  {%- endfor %}
+}
+
+func sec_count(/* SectionList */) { local nSec
+  nSec = 0
+  forsec $o1 {
+      nSec += 1
+  }
+  return nSec
+}
+
+/*
+ * Iterate over the section and compute how many segments should be allocate to
+ * each.
+ */
+proc geom_nseg_fixed(/* chunkSize */) { local secIndex, chunkSize
+  chunkSize = $1
+  soma area(.5) // make sure diam reflects 3d points
+  secIndex = 0
+  forsec all {
+    nseg = 1 + 2*int(L/chunkSize)
+    segCounts.x[secIndex] = nseg
+    secIndex += 1
+  }
+}
+
+/*
+ * Count up the number of sections
+ */
+proc geom_nsec() { local nSec
+  nSecAll = sec_count(all)
+  nSecSoma = sec_count(somatic)
+  nSecApical = sec_count(apical)
+  nSecBasal = sec_count(basal)
+  nSecMyelinated = sec_count(myelinated)
+  nSecAxonalOrig = nSecAxonal = sec_count(axonal)
+
+  segCounts = new Vector()
+  segCounts.resize(nSecAll)
+  nSec = 0
+  forsec all {
+    segCounts.x[nSec] = nseg
+    nSec += 1
+  }
+}
+
+/*
+ * Replace the axon built from the original morphology file with a stub axon
+ */
+{%- if replace_axon %}
+    {{replace_axon}}
+{%- endif %}
+
+
+{{re_init_rng}}
+
+endtemplate {{template_name}}