Clear[s, e, n, w, dir, tab, m, \[Delta], positions, allPrimes,
   twinPrimes, sph, ulam];
 m = 20; (* m is the no. of spiral "laps", so (2m)^2 numbers will be \
 displayed *)
 \[Delta] = .1; (* \[Delta] = vertical depth added to \
 pyramid for each number *)
 
 dir = {s, e, n, w}; (* directions: south, east, north, west *)
 
tab = {};
Do[dir = RotateLeft[dir]; AppendTo[tab, Table[First[dir], {k}]];
  dir = RotateLeft[dir];
  AppendTo[tab, Table[First[dir], {k}]], {k, 1, 2 m}];
positions =
  FoldList[Plus, {0, 0, 0},
   Take[Flatten[tab], 4 m^2 - 1] /. {e -> {\[Delta], -1, 0},
     n -> {\[Delta], 0, 1}, w -> {\[Delta], 1, 0},
     s -> {\[Delta], 0, -1}}];
allPrimes = PrimeQ[Range[4 m^2]];
twinPrimes =
  MapThread[
   Or, {MapThread[
     And, {allPrimes, Join[{False, False}, Drop[allPrimes, -2]]}],
    MapThread[
     And, {allPrimes, Join[Drop[allPrimes, 2], {False, False}]}]}];
sph[k_, pos_, highlight_] :=
  Translate[
   SphericalPlot3D[.3, {\[Theta], 0, Pi}, {\[Phi], 0, 2 Pi},
     Mesh -> None,
     PlotStyle ->
      Texture[Graphics[
        Style[Text[k, {0, 0}, {0, 0}, {0, 1}], 128,
         FontFamily -> "Optima"], AspectRatio -> 2,
        Background ->
         If[highlight, RGBColor[1.`, 1.`, 0.`],
          RGBColor[.5`, .5`, .5`]], ImageSize -> 400]],
     TextureCoordinateFunction -> ({#4, #5} &)][[1]], pos];
ulam = Graphics3D[{
    {Gray, Tube[positions, .05]},
    MapThread[sph, {Range[4 m^2], positions, twinPrimes}]
    }, ViewPoint -> {-1.8, 0, 0}, PlotRangePadding -> None,
   Boxed -> False, Lighting -> "Neutral", ImageSize -> 5200,
   Background -> Black];