diff --git a/doc/source/api_reference.rst b/doc/source/api_reference.rst
index 3fe6a38bf..17c1f1765 100644
--- a/doc/source/api_reference.rst
+++ b/doc/source/api_reference.rst
@@ -28,4 +28,4 @@ Relationships:
 .. testsetup:: *
 
     from neo import SpikeTrain
-    import quantities as pq
+    import neo.units
diff --git a/doc/source/read_and_analyze.rst b/doc/source/read_and_analyze.rst
index 30ca47842..57f29b5c6 100644
--- a/doc/source/read_and_analyze.rst
+++ b/doc/source/read_and_analyze.rst
@@ -223,7 +223,7 @@ objects directly, e.g.:
 
     In [20]: import numpy as np
 
-    In [21]: from quantities import mV, kHz
+    In [21]: from neo.units import mV, kHz
 
     In [22]: from neo import AnalogSignal
 
@@ -245,7 +245,7 @@ specify the array of times at which the signal was sampled.
 
 .. ipython::
 
-    In [25]: from quantities import ms, nA
+    In [25]: from neo.units import ms, nA
 
     In [26]: from neo import IrregularlySampledSignal
 
@@ -336,7 +336,7 @@ It behaves similarly to :class:`AnalogSignal`, but in 3D rather than 2D.
 
 .. ipython::
 
-    In [38]: from quantities import Hz, micrometer
+    In [38]: from neo.units import Hz, micrometer
 
     In [39]: from neo import ImageSequence
 
@@ -367,7 +367,7 @@ For any metadata not covered by the required or recommended fields, additional a
 
 .. ipython::
 
-    In [43]: from quantities import um as µm
+    In [43]: from neo.units import um as µm
 
     In [44]: signal.annotate(pipette_tip_diameter=1.5 * µm)
 
diff --git a/examples/plot_imageseq.py b/examples/plot_imageseq.py
index a482784c1..bd8bc4296 100644
--- a/examples/plot_imageseq.py
+++ b/examples/plot_imageseq.py
@@ -10,7 +10,7 @@
 from neo.core import ImageSequence
 from neo.core import RectangularRegionOfInterest, CircularRegionOfInterest, PolygonRegionOfInterest
 import matplotlib.pyplot as plt
-import quantities as pq
+import neo.units
 
 import random
 
@@ -33,7 +33,7 @@
 # we then make our image sequence and pull out our results from the
 # image_seq
 
-image_seq = ImageSequence(l, sampling_rate=500 * pq.Hz, spatial_scale="m", units="V")
+image_seq = ImageSequence(l, sampling_rate=500 * neo.units.Hz, spatial_scale="m", units="V")
 
 result = image_seq.signal_from_region(
     CircularRegionOfInterest(image_seq, 50, 50, 25),
diff --git a/examples/plot_multi_tetrode_example.py b/examples/plot_multi_tetrode_example.py
index 2480b88fe..353d1989a 100644
--- a/examples/plot_multi_tetrode_example.py
+++ b/examples/plot_multi_tetrode_example.py
@@ -10,7 +10,7 @@
 
 from itertools import cycle
 import numpy as np
-from quantities import ms, mV, kHz
+from neo.units import ms, mV, kHz
 import matplotlib.pyplot as plt
 from neo import Block, Segment, ChannelView, Group, SpikeTrain, AnalogSignal
 
diff --git a/examples/plot_read_proxy_with_lazy_load.py b/examples/plot_read_proxy_with_lazy_load.py
index 66b2777ff..3e466d338 100644
--- a/examples/plot_read_proxy_with_lazy_load.py
+++ b/examples/plot_read_proxy_with_lazy_load.py
@@ -11,7 +11,7 @@
 
 import urllib
 import neo
-import quantities as pq
+import neo.units
 import numpy as np
 
 
@@ -68,8 +68,8 @@ def apply_my_fancy_average(sig_list):
 # Let's set our limits for both cases. We will
 # use quantities to include time dimensions.
 
-lim_start = -20 * pq.ms  # 20 milliseconds before
-lim_end = +20 * pq.ms  # 20 milliseconds after
+lim_start = -20 * neo.units.ms  # 20 milliseconds before
+lim_end = +20 * neo.units.ms  # 20 milliseconds after
 
 ##################################################
 # We start with eager (where `lazy=False`.) Everything
diff --git a/examples/plot_roi_demo.py b/examples/plot_roi_demo.py
index 653947074..a69fbbbb6 100644
--- a/examples/plot_roi_demo.py
+++ b/examples/plot_roi_demo.py
@@ -15,7 +15,7 @@
 from neo.core import CircularRegionOfInterest, RectangularRegionOfInterest, PolygonRegionOfInterest, ImageSequence
 from numpy.random import rand
 import random
-import quantities as pq
+import neo.units
 
 ##################################################################
 # First we create our image_sequence. Let's generate some data
@@ -32,7 +32,7 @@
             l[frame][y].append(random.randint(0, 50))
 
 # make an ImageSquence in Neo
-image_seq = ImageSequence(l, sampling_rate=500 * pq.Hz, spatial_scale="m", units="V")
+image_seq = ImageSequence(l, sampling_rate=500 * neo.units.Hz, spatial_scale="m", units="V")
 
 
 #################################################################
diff --git a/examples/plot_with_matplotlib.py b/examples/plot_with_matplotlib.py
index 71c30bdd9..d46eae8c0 100644
--- a/examples/plot_with_matplotlib.py
+++ b/examples/plot_with_matplotlib.py
@@ -28,8 +28,8 @@
     ax2 = fig.add_subplot(2, 1, 2)
     ax1.set_title(seg.file_origin)
     ax1.set_ylabel("arbitrary units")
-    mint = 0 * pq.s
-    maxt = np.inf * pq.s
+    mint = 0 * neo.units.s
+    maxt = np.inf * neo.units.s
     for i, asig in enumerate(seg.analogsignals):
         times = asig.times.rescale("s").magnitude
         asig = asig.magnitude
diff --git a/neo/__init__.py b/neo/__init__.py
index 44c463041..581ce4939 100644
--- a/neo/__init__.py
+++ b/neo/__init__.py
@@ -14,3 +14,4 @@
 
 from neo.core import *
 from neo.io import *
+import neo.units
diff --git a/neo/io/asciiimageio.py b/neo/io/asciiimageio.py
index 69a5c0531..2f35d5e4d 100644
--- a/neo/io/asciiimageio.py
+++ b/neo/io/asciiimageio.py
@@ -9,9 +9,9 @@ class AsciiImageIO(BaseIO):
 
     *Usage*:
         >>> from neo import io
-        >>> import quantities as pq
+        >>> import neo.units
         >>> r = io.AsciiImageIO(file_name='File_asciiimage_1.txt',nb_frame=511, nb_row=100,
-        ...                     nb_column=100,units='mm', sampling_rate=1.0*pq.Hz,
+        ...                     nb_column=100,units='mm', sampling_rate=1.0*neo.units.Hz,
         ...                     spatial_scale=1.0*pq.mm)
         >>> block = r.read_block()
         read block
diff --git a/neo/io/elphyio.py b/neo/io/elphyio.py
index 76477281e..998f53185 100644
--- a/neo/io/elphyio.py
+++ b/neo/io/elphyio.py
@@ -52,8 +52,8 @@
 filled with other NEO structures:
 
 >>> import numpy as np
->>> import quantities as pq
->>> a = AnalogSignal( signal=np.random.rand(300), t_start=42*pq.ms)
+>>> import neo.units
+>>> a = AnalogSignal( signal=np.random.rand(300), t_start=42*neo.units.ms)
 >>> s.analogsignals.append( a )
 
 and added to a newly created NEO Block:
diff --git a/neo/io/neomatlabio.py b/neo/io/neomatlabio.py
index 0e575fdd4..75b764cb5 100644
--- a/neo/io/neomatlabio.py
+++ b/neo/io/neomatlabio.py
@@ -156,7 +156,7 @@ class NeoMatlabIO(BaseIO):
         This Python code generates the same block as in the previous scenario::
 
             import neo
-            import quantities as pq
+            import neo.units
             from scipy import rand, array
 
             bl = neo.Block(name='my block with neo')
@@ -164,14 +164,14 @@ class NeoMatlabIO(BaseIO):
                 seg = neo.Segment(name='segment' + str(s))
                 bl.segments.append(seg)
                 for a in range(5):
-                    anasig = neo.AnalogSignal(rand(100)*pq.mV, t_start=0*pq.s,
+                    anasig = neo.AnalogSignal(rand(100)*neo.units.mV, t_start=0*neo.units.s,
                                               sampling_rate=100*pq.Hz)
                     seg.analogsignals.append(anasig)
                 for t in range(7):
-                    sptr = neo.SpikeTrain(rand(40)*pq.ms, t_start=0*pq.ms, t_stop=10*pq.ms)
+                    sptr = neo.SpikeTrain(rand(40)*pq.ms, t_start=0*neo.units.ms, t_stop=10*neo.units.ms)
                     seg.spiketrains.append(sptr)
-                ev = neo.Event([0, 10, 30]*pq.ms, labels=array(['trig0', 'trig1', 'trig2']))
-                ep = neo.Epoch([10, 20]*pq.ms, durations=[4, 10]*pq.ms, labels=array(['a0', 'a1']))
+                ev = neo.Event([0, 10, 30]*neo.units.ms, labels=array(['trig0', 'trig1', 'trig2']))
+                ep = neo.Epoch([10, 20]*neo.units.ms, durations=[4, 10]*neo.units.ms, labels=array(['a0', 'a1']))
                 seg.events.append(ev)
                 seg.epochs.append(ep)
 
diff --git a/neo/io/nestio.py b/neo/io/nestio.py
index f2ed17b65..561755a43 100644
--- a/neo/io/nestio.py
+++ b/neo/io/nestio.py
@@ -32,12 +32,13 @@ class NestIO(BaseIO):
 
     Usage:
         >>> from neo.io.nestio import NestIO
+        >>> import neo.units
 
         >>> files = ['membrane_voltages-1261-0.dat',
                  'spikes-1258-0.gdf']
         >>> r = NestIO(filenames=files)
-        >>> seg = r.read_segment(gid_list=[], t_start=400 * pq.ms,
-                             t_stop=600 * pq.ms,
+        >>> seg = r.read_segment(gid_list=[], t_start=400 * neo.units.ms,
+                             t_stop=600 * neo.units.ms,
                              id_column_gdf=0, time_column_gdf=1,
                              id_column_dat=0, time_column_dat=1,
                              value_columns_dat=2)
diff --git a/neo/io/proxyobjects.py b/neo/io/proxyobjects.py
index 7f372bfb2..46ba53685 100644
--- a/neo/io/proxyobjects.py
+++ b/neo/io/proxyobjects.py
@@ -83,7 +83,7 @@ class AnalogSignalProxy(BaseProxy):
                                          block_index=0,
                                          seg_index=0)
     >>> anasig = proxy_anasig.load()
-    >>> slice_of_anasig = proxy_anasig.load(time_slice=(1.*pq.s, 2.*pq.s))
+    >>> slice_of_anasig = proxy_anasig.load(time_slice=(1.*neo.units.s, 2.*neo.units.s))
     >>> some_channel_of_anasig = proxy_anasig.load(channel_indexes=[0,5,10])
 
     """
diff --git a/neo/io/tiffio.py b/neo/io/tiffio.py
index f9370da73..fdccbd5fc 100644
--- a/neo/io/tiffio.py
+++ b/neo/io/tiffio.py
@@ -17,9 +17,9 @@ class TiffIO(BaseIO):
 
     *Usage*:
         >>> from neo import io
-        >>> import quantities as pq
-        >>> r = io.TiffIO("dir_tiff",spatial_scale=1.0*pq.mm, units='V',
-        ...               sampling_rate=1.0*pq.Hz)
+        >>> import neo.units
+        >>> r = io.TiffIO("dir_tiff",spatial_scale=1.0*neo.units.mm, units='V',
+        ...               sampling_rate=1.0*neo.units.Hz)
         >>> block = r.read_block()
         read block
         creating segment
diff --git a/neo/units/__init__.py b/neo/units/__init__.py
new file mode 100644
index 000000000..9b9776068
--- /dev/null
+++ b/neo/units/__init__.py
@@ -0,0 +1 @@
+from quantities.units import *
\ No newline at end of file