diff --git a/lems/parser/LEMS.py b/lems/parser/LEMS.py
index 2cdb725..614b560 100644
--- a/lems/parser/LEMS.py
+++ b/lems/parser/LEMS.py
@@ -9,14 +9,12 @@
 
 from lems.base.base import LEMSBase
 from lems.base.errors import ParseError
-
-from lems.model.fundamental import *
+from lems.base.util import make_id
 from lems.model.component import *
 from lems.model.dynamics import *
-from lems.model.structure import *
+from lems.model.fundamental import *
 from lems.model.simulation import *
-
-from lems.base.util import make_id
+from lems.model.structure import *
 
 
 def get_nons_tag_from_node(node):
@@ -209,9 +207,9 @@ def init_parser(self):
         self.tag_parse_table["eventwriter"] = self.parse_event_writer
         self.tag_parse_table["derivedparameter"] = self.parse_derived_parameter
         self.tag_parse_table["derivedvariable"] = self.parse_derived_variable
-        self.tag_parse_table[
-            "conditionalderivedvariable"
-        ] = self.parse_conditional_derived_variable
+        self.tag_parse_table["conditionalderivedvariable"] = (
+            self.parse_conditional_derived_variable
+        )
         self.tag_parse_table["case"] = self.parse_case
         self.tag_parse_table["dimension"] = self.parse_dimension
         self.tag_parse_table["dynamics"] = self.parse_dynamics
@@ -1072,19 +1070,25 @@ def parse_include(self, node):
 
         :raises ParseError: Raised when the file to be included is not specified.
         """
+        filename = None
+
+        # TODO: remove this hard coding for reading NeuroML includes...
+        if "file" not in node.lattrib:
+            if "href" in node.lattrib:
+                filename = node.lattrib["href"]
+        else:
+            filename = node.lattrib["file"]
+
+        if filename is None:
+            self.raise_error("<Include> must specify the file to be included.")
+
         if not self.include_includes:
             if self.model.debug:
-                print("Ignoring included LEMS file: %s" % node.lattrib["file"])
+                self.model.add_include(Include(filename))
+                print("Not including contents of %s" % filename)
         else:
-            # TODO: remove this hard coding for reading NeuroML includes...
-            if "file" not in node.lattrib:
-                if "href" in node.lattrib:
-                    self.model.include_file(node.lattrib["href"], self.include_dirs)
-                    return
-                else:
-                    self.raise_error("<Include> must specify the file to be included.")
-
-            self.model.include_file(node.lattrib["file"], self.include_dirs)
+            self.model.include_file(filename, self.include_dirs)
+            print("Included contents of %s" % filename)
 
     def parse_kinetic_scheme(self, node):
         """
diff --git a/lems/test/LEMS_NML2_Ex2_Izh.xml b/lems/test/LEMS_NML2_Ex2_Izh.xml
new file mode 100644
index 0000000..b1cefaf
--- /dev/null
+++ b/lems/test/LEMS_NML2_Ex2_Izh.xml
@@ -0,0 +1,78 @@
+<Lems>
+
+        <!-- Example with Bursting Izhikevich cell -->
+
+        <!-- This is a file which can be read and executed by the LEMS Interpreter.
+             It imports the LEMS definitions of the core NeuroML 2 Components, creates
+             a model in "pure" NeuroML 2 and contains some LEMS elements for running
+             a simulation -->
+
+
+<!-- Specify which component to run -->
+<Target component="sim1"/>
+
+<!-- Include core NeuroML2 ComponentType definitions -->
+<Include file="Cells.xml"/>
+<Include file="Networks.xml"/>
+<Include file="Inputs.xml"/>
+<Include file="Simulation.xml"/>
+
+
+        <!-- Main NeuroML2 content. Based on http://www.izhikevich.org/publications/figure1.m -->
+
+<izhikevichCell id="izBurst" v0 = "-70mV" thresh = "30mV" a ="0.02" b = "0.2" c = "-50" d = "2"/>
+
+
+
+<pulseGeneratorDL id="i0" delay="22ms" duration="2000ms" amplitude="15" />
+
+<izhikevichCell id="izTonic" v0 = "-70mV" thresh = "30mV" a ="0.02" b = "0.2" c = "-65" d = "6"/>
+<pulseGeneratorDL id="i1" delay="20ms" duration="2000ms" amplitude="14" />
+
+<izhikevichCell id="izMixed" v0 = "-70mV" thresh = "30mV" a ="0.02" b = "0.2" c = "-55" d = "4"/>
+<pulseGeneratorDL id="i2" delay="20ms" duration="2000ms" amplitude="10" />
+
+<izhikevichCell id="izClass1" v0 = "-60mV" thresh = "30mV" a ="0.02" b = "-0.1" c = "-55" d = "6"/>
+<rampGeneratorDL id="rg0" delay="30ms" duration="170ms" startAmplitude="-32" finishAmplitude="50" baselineAmplitude="-32"/>
+
+<network id="net1">
+    <population id="izpopBurst" component="izBurst" size="1"/>
+    <population id="izpopTonic" component="izTonic" size="1"/>
+    <population id="izpopMixed" component="izMixed" size="1"/>
+    <population id="izpopClass1" component="izClass1" size="1"/>
+
+    <explicitInput target="izpopBurst[0]" input="i0" destination="synapses"/>
+    <explicitInput target="izpopTonic[0]" input="i1" destination="synapses"/>
+    <explicitInput target="izpopMixed[0]" input="i2" destination="synapses"/>
+    <explicitInput target="izpopClass1[0]" input="rg0" destination="synapses"/>
+</network>
+
+        <!-- End of NeuroML2 content -->
+
+
+<Simulation id="sim1" length="200ms" step="0.005ms" target="net1">
+
+    <Display id="d1" title="Ex2: Bursting Izhikevich cell in LEMS" timeScale="1ms" xmin="-20" xmax="220" ymin="-80" ymax="40">
+        <Line id ="v" quantity="izpopBurst[0]/v" scale="1mV"  color="#ee40FF" timeScale="1ms"/>
+        <Line id ="U" quantity="izpopBurst[0]/U" scale="1"  color="#BBA0AA" timeScale="1ms"/>
+        <Line id ="i" quantity="izpopBurst[0]/i0/I" scale="1"  color="#222222" timeScale="1ms"/>
+    </Display>
+
+    <Display id="d2" title="Ex2: Tonic spiking Izhikevich cell in LEMS" timeScale="1ms" xmin="-20" xmax="220" ymin="-80" ymax="40">
+        <Line id ="v" quantity="izpopTonic[0]/v" scale="1mV"  color="#ee40FF" timeScale="1ms"/>
+        <Line id ="U" quantity="izpopTonic[0]/U" scale="1V"  color="#BBA0AA" timeScale="1ms"/>
+        <Line id ="i" quantity="izpopTonic[0]/i1/I" scale="1"  color="#222222" timeScale="1ms"/>
+    </Display>
+
+    <Display id="d3" title="Ex2: Mixed mode Izhikevich cell in LEMS" timeScale="1ms" xmin="-20" xmax="220" ymin="-80" ymax="40">
+        <Line id ="v" quantity="izpopMixed[0]/v" scale="1mV"  color="#ee40FF" timeScale="1ms"/>
+        <Line id ="U" quantity="izpopMixed[0]/U" scale="1V"  color="#BBA0AA" timeScale="1ms"/>
+        <Line id ="i" quantity="izpopMixed[0]/i2/I" scale="1"  color="#222222" timeScale="1ms"/>
+    </Display>
+
+    <!-- See other examples for saving of data traces -->
+
+</Simulation>
+
+
+</Lems>
diff --git a/lems/test/test_load_write.py b/lems/test/test_load_write.py
index bdea25e..3952ac6 100644
--- a/lems/test/test_load_write.py
+++ b/lems/test/test_load_write.py
@@ -3,7 +3,9 @@
 
 """
 
+import os
 import unittest
+
 from lems.model.model import Model
 
 
@@ -30,3 +32,47 @@ def test_load_get_dom(self):
         file_name = "lems/test/hhcell_resaved2.xml"
         model.import_from_file(file_name)
         dom0 = model.export_to_dom()
+
+    def test_include_includes_is_true(self):
+        """Test that include_includes works as expected"""
+        cwd = os.getcwd()
+        os.chdir("lems/test/")
+        model = Model(
+            include_includes=True,
+            fail_on_missing_includes=True,
+        )
+        model.debug = True
+        model.add_include_directory("NeuroML2CoreTypes/")
+
+        model.import_from_file("LEMS_NML2_Ex2_Izh.xml")
+
+        model_string = model.export_to_dom().toprettyxml("    ", "\n")
+        print(model_string)
+        self.assertNotIn("<Include ", model_string)
+
+        self.assertEqual(0, len(model.includes))
+
+        os.chdir(cwd)
+
+    def test_include_includes_is_false(self):
+        """Test that include_includes works as expected"""
+        cwd = os.getcwd()
+        os.chdir("lems/test/")
+        model = Model(
+            include_includes=False,
+            fail_on_missing_includes=True,
+        )
+        model.debug = True
+        model.add_include_directory("NeuroML2CoreTypes/")
+
+        model.import_from_file("LEMS_NML2_Ex2_Izh.xml")
+
+        model_string = model.export_to_dom().toprettyxml("    ", "\n")
+        print(model_string)
+        self.assertIn("<Include ", model_string)
+
+        for include in model.includes:
+            inc_string = include.toxml()
+            self.assertIn("<Include file=", inc_string)
+
+        os.chdir(cwd)