Which unfoldings of the hypercube tile 3-space: How to check for isometric space-fillers? Recently Mark McClure constructed and displayed 
the 261 unfoldings of the hypercube (tesseract)
in response to the question,
"3D models of the unfoldings of the hypercube?":



The first 9 unfoldings in Mark McClure's display

Each of the 11 unfoldings of the cube form monohedral tilings of the plane,
as so well illustrated in the "Etudes" video
to which Igor Pak pointed:

         


A polyhedron that is the prototile of a monohedral tiling is called
an isometric space-filler: $\mathbb{R}^3$ can be tiled by congruent copies
of that one shape (rotated and translated but not reflected).

Now that we have the unfoldings of the hypercube, analogy with the cube raises the
question:

Q1. Which (if any) of the 261 unfoldings of the hypercube are isometric space-fillers?

Asking this question raises another:

Q2. How can one determine if a given shape, in this case a polycube / 3D polyomino, is an isometric space-filler?


Update (7Dec2015).
Aside from the two hypercube unfoldings that Steven Stadnicki showed
tile space (below), with a student I found two more that tile $\mathbb{R}^3$,
including the Dali hypercube cross unfolding,
confirming Steven's intuition ("I don't know if the 'Dali unfolding' tiles space, though I'd be surprised if it didn't.")
We posted an arXiv note on the topic.
 A: Here is an illustration of Steven Stadnicki's first tiling
of a (portion of a) "2-cell thick infinite plane":

 


A: For your first question, I haven't gone through all of the unfoldings, but the first one there clearly tiles space: it even tiles a 2-cell thick infinite plane (which can then be stacked to tile space), something like this:
$$\begin{array}{cccc}
a & b & c & d & e & f \\
a & g & c & d & e & h \\
a & g & i & d & e & h \\
a & g & i & j & e & h \\
k & g & i & j & l & h \\
k & m & i & j & l & n \\
\\
\\
a & a & a & a & e & e \\
o & g & g & g & g & h \\
p & p & i & i & i & i \\
q & q & q & j & j & j \\
k & k & k & k & l & l \\
r & m & m & m & m & n \\ 
\end{array}$$
As for your second question, the traditional answer is 'find a region that the polyomino tiles, which in turn tiles the space by translation'.  A good example of this is the second unfolding shown here.  We can couple the given shape with a copy that's been rotated 180 degrees about two axes to get a more 'solid' shape with a rotational symmetry:
$$\begin{array}
. & . & a & a & . & . \\
a & a & a & b & b & b \\
\\
\\
b & b & b & a & a & a \\
. & . & b & b & . & . \\
\end{array}$$
Now, a copy of this shape can be paired with another, rotated copy to form a shape that can tile a 2x2xinfinity 'rod' by translation:
$$\begin{array}
. & . & c & c & d & d & d & d & d & d \\
c & c & c & c & c & c & d & d & . & . \\
\\
\\
c & c & c & c & c & c & d & d & . & . \\
. & . & c & c & d & d & d & d & d & d \\
\end{array}$$
Something similar (essentially, building shapes that fit an extended Conway criterion) can certainly be done with some of the other 'biplanar' unfoldings, but I haven't gone through them to try and figure out which ones are amenable to it.  I don't know if the 'Dali unfolding' tiles space, though I'd be surprised if it didn't.
Meanwhile, it's at least conceivable that one of the unfoldings tiles space only aperiodically, but that would be a shock; AFAIK, no single polyomino has been shown to tile only aperiodically (though sets of polyominoes can model Wang tilings, and so tile only aperiodically) in either two or three dimensions.
A: Answer to Q1: All of the 261. 
I looked at this question because of a video of Matt Parker and wrote an algorithm to find solutions. See here for an example of how a solution would look like. I dumped all solutions on github. For some cases I list multiple solutions. The files start with a number followed by an underscore and the number is between 0 and 260 corresponding to one of the unfoldings in this list.
129   <--- Number of unfolding (this time 0-indexed)
[[-2, -2, 0], [1, 1, 1]]  <---- the box, we put it in
Boolean Program (maximization, 1152 variables, 428 constraints)  <---
some information to ignore
32.0  <--- How many of the voxels in the cube are covered (here it
happens to be all of them)
[(-2, -2, 0), (-1, -2, 0), (-1, -2, 1), (-1, -1, 0), (-1, 0, 0), (-1,
1, 0), (0, -2, 1), (1, -2, 1)]
[(-2, -2, 1), (-2, -1, 0), (-2, -1, 1), (-2, 0, 0), (-2, 1, 0), (-1,
-1, 1), (0, -1, 1), (1, -1, 1)]
[(-2, 0, 1), (-1, 0, 1), (0, 0, 1), (1, -2, 0), (1, -1, 0), (1, 0, 0),
(1, 0, 1), (1, 1, 1)]
[(-2, 1, 1), (-1, 1, 1), (0, -2, 0), (0, -1, 0), (0, 0, 0), (0, 1, 0),
(0, 1, 1), (1, 1, 0)]
 

(the last rows are then the coordinates  of the copies of the
unfolding, whatever sticks out of the box fits in another one..)
Here is an example of one of the two one that fit in a 4x4x2 box ( in an exploded view, you can shift them together and they neatly fit into a 4x4x2 box, which can then be stacked to obtain a tiling.:

Answer to Q2: One way is to use integer programming.
The basic idea is to take a box with volume divisible by 8 as a fundamental domain and to cover it as much as possible with non-overlapping copies of an unfolding and make sure that the stuff that sticks out of the bottom actually fits into gaps at the top and so on.
Here's a two-dimensional example:
The tile:

A 5x6 box filled with copies of that tile:

A tiling produced from this:

This can be formulated a as an integer program and it turns out that those 261 unfoldings all have feasible solution for some smallish box (for most 4x4x4 is enough).
Setting up the integer program is just a few lines of sage:
from sage.combinat.tiling import Polyomino
import itertools
def get_mod_points(point, diff):
    for coeffs in itertools.product([-1, 0, 1], repeat=diff.length()):
        yield vector([diff[i]*coeffs[i] + point[i] for i in range(diff.length())])
def milp_mod_box(tile, inner_box, outer_box, solver=None):
    # tile, inner_box and outer_box are `Polyomino`.
    begin, end = inner_box.bounding_box()
    diff = vector(end) - vector(begin)
    diff += vector([1]*diff.length())

    tiles = [Q for Q in tile.isometric_copies(outer_box)]
    p = MixedIntegerLinearProgram(solver=solver, maximization=True)
    x = p.new_variable(binary=True)
    for i in inner_box:
        cons = p.sum(sum([int(j in T) for j in get_mod_points(i, diff)])*x[T] for T in tiles) == 1
        p.add_constraint(cons)
    for i in outer_box:
        p.add_constraint(p.sum(x[T] for T in tiles if i in T) <= 1)
    p.set_objective(p.sum(x[T] for i in inner_box for T in tiles if i in T))
    
    return p, x  

Here we optimize for having as much voxels filled in the box as possible, but only so that the solutions look nicer, it would give valid tilings even without setting that objective.
A similar integer programming approach could be used to prove that a certain polyomino is not a space filler: We can check for larger and larger boxes that if it is not possible to cover them completely with non-overlapping copies of the original polyomino, but this wasn't necessary.
Update:
I added 3d-plots of all solutions to 3d-renderings of all 261 unfoldings page. When you click on the number of each one you get a 3d-rending of copies of that unfolding mostly in a box (drawn in grey), which then can be used to tile the plane.
Note this is a different order than the one above (and it starts indexing with $1$.

*

*#213 is the Dalí-unfolding, needs 8 pieces in the prototile.


*#3 is the only one that only needs 3 pieces (in a 4x3x2 box) (what a coincidence..)


*#72 and #159 are the ones that perfectly fit in a 4x4x2 box.


*#139 needs a long 8x2x2 box.
