Python 3.8 (pre-release), 286 bytes

A lazy solution.

print("".join([['  ','.'*16+'\n','.','   .','###.',' ','\n','.   '*3+'.'][int(y)]for y in oct(int('6B2GOZIS2BMT2MWC2GYF3PKVX707G4YOJQN0CVHDAN6PY30YREP0CFNWVQSJBNEUUQI4L81LMK1FATS9FEZM03WF80ZGCXD9RENRA79O7AYCTGFAW0PC8BXIM2ZB9S0U6B7W1RS5EQ82ZV3TZQQ34H4LDA1SQVJ3IQ8GDPD4K62O1B',36))[2:]]))

Try it online!

Explanation

Hardcoding the output in the form of an octal number, where the digits 0-7 correspond to the string at the corresponding index in the list [' ','.'*16+'\n','.',' .','###.',' ','\n','. '*3+'.']. The octal number is stored in base 36, and then converted to octal at runtime using oct().

Ruby, 120 103 102 100 bytes

68.times{|j|puts j%4<1??0*18:("DHN@@JLD@HJL@DXL@"[j/4].ord.digits(2).map{|i|' 1'[i]*3}*?0)[j/17,18]}

Try it online!

Uses space for white, 0 for grey and 1 for black. The latest version does not use arithmetic to generate the output characters so now any characters may be used.

We print a grid of characters 18 wide x 17 high for each stage, total 68 lines. The fact that the grids are 17 high means the first and last line of each grid are the same, and every fourth line of the whole output is all zeros, with the data shifted one cell right every 17 lines.

The black/white cells to plot comprise 17 rows and 5 columns of data, the leftmost of which is always 0. Each of the 17 characters in the magic string represents 1 row. In previous versions, bits were extracted and multiplied by 14, giving the binary representation " 0" for 0 and "1110" for 1, with the 0 forming the right wall of the cell. The latest version generates " " or "111" by selecting a character from an array and concatenates with 0.

For each line we compile a string of characters from 7 bits (including the 1 from the 64's bit which is always set for characters in the selected range of the magic string but is not used) and pick 18 characters to print by [j/17,18]

It's a little difficult to see the boundary between the first and second grids, but the boundary above 3rd and 4th grid is clear.

Using 68.times{|j|puts...} is shorter overall than returning from a function with ->{(0..67).map{|j|...}}. Changing puts to p would save 3 bytes but add " at the beginning and end of each line.

Commented code

68.times{|j|                               #Iterate through 68 lines.
 puts j%4<1?                               #If line number mod 4 0
  ?0*18:                                   #output a row of 18 zeros, else
  ("DHN@@JLD@HJL@DXL"[j/4].ord.digits(2).  #select row data from magic string and iterate through bits.
   map{|i|' 1'[i]*3}*?0)[j/17,18]          #Select ' ' or '1' according to bit value and triplicate. Concatenate with '0'
}                                          #select 18 characters from the string starting at j/17 and output.

JavaScript (ES6), 116 bytes

-2 thanks to @l4m2

Returns an array of 4 matrices of integers of size 18×18. Uses 0, 1, 2 for white, black, grey.

_=>[1858,9808,25920,27680].map((v,n)=>[..._+[]+_].map((_,y,a)=>a.map((_,x)=>x+n&3&&y+n&3?v>>(y+n&~3|x+n-4>>2)&1:2)))

Try it online! (raw output)

Try it online! (beautified)

Animated version

f=
_=>[1858,9808,25920,27680].map((v,n)=>[..._+[]+_].map((_,y,a)=>a.map((_,x)=>x+n&3&&y+n&3?v>>(y+n&~3|x+n-4>>2)&1:2)))
a = f(); n = 0; (tick = _ => { document.getElementById("p").innerHTML = a[n++ & 3].map(r => r.map(v => " X:"[v]).join(" ")).join("\n"); setTimeout(tick, 500 - document.getElementById("r").value); })();

<pre id="p" style="font-size:8px;float:left"></pre><input type="range" id="r" min=0 max=450 value=250 style="float:left;margin-left:10px" />

Commented

_ =>                     // main function taking no input
[                        // glider patterns as 16-bit masks:
  1858,                  //   (0) 0000 0111 0100 0010
  9808,                  //   (1) 0010 0110 0101 0000
  25920,                 //   (2) 0110 0101 0100 0000
  27680                  //   (3) 0110 1100 0010 0000
]                        //
.map((v, n) =>           // for each mask v at index n:
  [..._ + [] + _]        //   build an array of length 18
  .map((_, y, a) =>      //   for y = 0 to 17:
    a.map((_, x) =>      //     for x = 0 to 17:
      x + n & 3 &&       //       if x + n is not a multiple of 4
      y + n & 3 ?        //       and y + n is not a multiple of 4:
        v >> (           //         right-shift the mask to
          y + n & ~ 3 |  //         the row at position y + n and
          x + n - 4 >> 2 //         the column at position x + n - 4
        ) & 1            //         isolate the relevant bit
      :                  //       else (this is a border):
        2                //         return 2
    )                    //     end of map()
  )                      //   end of map()
)                        // end of map()

Charcoal, 53 bytes

ＥＩ⪪”←＆⊞δ¡｜8Ｙ⧴)”⁴Ｅ¹⁸⭆¹⁶⎇∧﹪⁺λκ⁴﹪⁺νκ⁴§ #÷ιＸ²↨⁴⊖÷⁺⟦λν⟧κ⁴.

Try it online! Link is to verbose version of code. Explanation: Hardcodes the frames.

ＥＩ⪪”←＆⊞δ¡｜8Ｙ⧴)”⁴

Mapping over the four values 810, 865, 1858, 9808...

Ｅ¹⁸⭆¹⁶

... drawing an 18×16 grid...

⎇∧﹪⁺λκ⁴﹪⁺νκ⁴§ #÷ιＸ²↨⁴⊖÷⁺⟦λν⟧κ⁴.

... if the current character, adjusting for the outer index, is a cell border then output a . otherwise index into the binary of the current value to determine whether to output a # or a space.

I tried to write something a bit more Charcoal-y but that weighed in at a massive 61 bytes:

ＵＢ“∨j⁸ＳＯ”Ｆ⁴«Ｍι↖Ｂ¹⁸ψＦ⪪§⪪”)“∧dδ◧Ｐ⊘“*σΣ ”χ鲫Ｊ⊗⊗§κ⁰⊗⊗§κ¹ＵＯ³#»Ｄ⎚Ｄ

Attempt This Online! Link is to verbose version of code. Explanation:

ＵＢ“∨j⁸ＳＯ”

Set the background to be spaces and .s for the cell borders.

Ｆ⁴«Ｍι↖Ｂ¹⁸ψ

Each pass through the loop, draw a background box one square further up and left. The background characters are taken relative to their position, so that e.g. -1, -1 is always a cell border, because the last character in the background string was a .. This causes the background to shift once square each time.

Ｆ⪪§⪪”)“∧dδ◧Ｐ⊘“*σΣ ”χι²

Loop over one of the four glider shapes.

«Ｊ⊗⊗§κ⁰⊗⊗§κ¹ＵＯ³#»

Draw the shape.

Ｄ⎚Ｄ

Output the canvas and clear it ready for the next frame.

Jelly,  41 38  37 bytes

“¡¢1ʠḣ7~ȷḶ’B’Ẓ;$4ÐƤ⁺x4ZỊJm4Ɗ¦ṙʋ⁺ṛɗ"J$

A niladic Link that yields a list of four lists of 32 lists of 32 integers, where dead (white) is 0, border (grey) is 1, and alive (black) is -1.

Try it online! Footer pretty prints with as dead, . as borders, and X as alive.

Or see the smaller area of the previous 38 byte version.

How?

“...’B’Ẓ;$4ÐƤ⁺x4ZỊJm4Ɗ¦ṙʋ⁺ṛɗ"J$ - Link: no arguments
“...’                           - 15393229446195444429
     B                          - to binary
      ’                         - subtract one from each
          4ÐƤ                   - for each 4-slice:
         $                      -   last two links as a monad:
       Ẓ                        -     is prime? (vectorises) -> all zeros
        ;                       -     concatenate
             ⁺                  - repeat that
                              $ - last two links as a monad:
                             J  -   indices -> I = [1,2,3,4]
                           ɗ"   -   zip with last three links as a dyad - f(Frame, i):
                        ʋ       -     last four links as a dyad:
              x4                -       repeat each element of each four times
                Z               -       transpose
                      ¦         -       apply to...
                     Ɗ          -       ...indices: last three links as a monad:
                  J             -                     indices
                   m            -                     mod slice...
                    4           -                     ...four
                 Ị              -       ...action: insignificant? (always 1 here)
                       ṙ        -       rotate left by {i}
                          ṛ     -     with the same right argument, i...
                         ⁺      -     ...repeat -> f(AlteredFrame, i)