Here is a probabilistic proof. Start with the integral
\begin{align}
I=\int_{0}^{1} ... \int_{0}^{1} \frac{1}{1-x^2_1 \ ... \ x^2_4} \ dx_1 \ ... \ dx_4.
\end{align}
Convert the integrand into a geometric series, use Tonelli's Theorem to exchange summation and integration, and integrate term by term to get the sum
\begin{align}
\sum_{n=0}^{\infty} \frac{1}{(2n+1)^4}.
\end{align}
By expanding $\zeta(4)$ into
\begin{align}
\zeta(4) & =\sum_{n=1}^{\infty} \frac{1}{(2n)^4} + I \\
& = \frac{1}{16} \zeta(4) +I,
\end{align}
we see
\begin{align}
\zeta(4)= \frac{16}{15} \sum_{n=0}^{\infty} \frac{1}{(2n+1)^4}.
\end{align}
Reconsidering $I,$ make the change of variables (discovered by Beukers, Calabi, and Kolk)
\begin{align}
x_i=
\frac{\sin \left(\frac{\pi}{2} u_i \right)}{\cos \left(\frac{\pi}{2} u_{i+1} \right)}
\end{align}
for each $i \in \lbrace 1, \ ... \ 4 \rbrace$ and $u_{4+1} := u_1.$ Such a transformation has a Jacobian Determinant
\begin{align}
\left | \frac{\partial(x_1,\ ... \ x_4)}{\partial(u_1,\ ... \ u_4)} \right |=\left( \frac{\pi}{2} \right)^4 \left( 1-x^2_1 \ ... \ x^2_4 \right),
\end{align}
and the region of integration is the open polytope
\begin{align}
\Delta^{4}= \lbrace (u_1, \ ... \ ,u_4) \in \mathbb{R}^4 : u_{i}+u_{i+1} < 1 , u_i>0 , 1 \leq i \leq 4 \rbrace,
\end{align}
whose proofs can be found in
https://pdfs.semanticscholar.org/35be/01e63c0bfd32b82c97d58ccc9c35471c3617.pdf
and
https://www.maa.org/sites/default/files/pdf/news/Elkies.pdf.
Hence, we see \begin{align}
I=\left( \frac{\pi}{2}\right ) ^4 \text{Vol}(\Delta^4).
\end{align}
Consider independent random variables $U_1, \ ... \ ,U_4 \sim \text{Unif}(0,1).$ Then we see
\begin{align}
\text{Vol}(\Delta^4)= \text{Pr} \left(U_1+U_2, \ ... \ , U_4 + U_1 < 1 \right).
\end{align}
First suppose all $U_i < \frac{1}{2}.$ Then we see $U_{i} +U_{i+1}<1$ for each $i,$ which tells us that $\text{Pr} \left(U_1+U_2, \ ... \ , U_4 + U_1 < 1 \right),$ assuming each $U_i < \frac{1}{2},$ is
\begin{align}
\int_{0}^{\frac{1}{2}} ... \int_{0}^{\frac{1}{2}} \ du_1 \ ... \ du_4 = \frac{1}{16},
\end{align}
the volume of the $4$ dimensional cube $\left(0, \frac{1}{2} \right)^4.$
Next, suppose exactly one $U_j \geq \frac{1}{2}.$ Then we see that $U_{j-1},U_{j+1}<1-U_{j}$ and the remaining $U_l<\frac{1}{2}.$ Note that since we have $4$ random variables, there in all $4$ ways this can happen. This tells us that $\text{Pr} \left(U_1+U_2, \ ... \ , U_4 + U_1 < 1 \right),$ assuming exactly one $U_j \geq \frac{1}{2},$ is
\begin{align}
4 \int_{\frac{1}{2}}^{1} \int_{0}^{1-u_j} \int_{0}^{1-u_j} \int_{0}^{\frac{1}{2}} \ du_l \ du_{j-1} \ du_{j+1} \ du_j = \frac{1}{12}.
\end{align}
Lastly, suppose exactly two $U_j, U_l \geq \frac{1}{2}$ and $U_j \geq U_l.$ Then we see $U_{j-1},U_{j+1}<1-U_j.$ Note that $j$ and $l$ must not be consecutive to one another. There are $2$ ways to choose such pairs $j$ and $l,$ and for each pair, there are $2$ ways of ordering $U_j,U_l,$ so there are $4$ instances of this happening. Thus, we have
that $\text{Pr} \left(U_1+U_2, \ ... \ , U_4 + U_1 < 1 \right),$ assuming exactly two $U_j,U_l \geq \frac{1}{2},$ is
\begin{align}
4 \int_{\frac{1}{2}}^{1} \int_{0}^{1-u_j} \int_{0}^{1-u_j} \int_{\frac{1}{2}}^{u_j} \ du_l \ du_{j-1} \ du_{j+1} \ du_j = \frac{1}{48}.
\end{align}
We cannot have more than two random variables simultaneously be at least $\frac{1}{2},$ which would lead to a contradiction to our constraints. Hence, summing the individual probabilities gives
\begin{align}
\text{Vol}(\Delta^4)= \frac{1}{6},
\end{align}
so
\begin{align}
I= \left( \frac{\pi}{2} \right)^4 \frac{1}{6} = \frac{\pi^4}{96},
\end{align}
and
\begin{align}
\zeta(4)= \frac{16}{15} I = \frac{\pi^4}{90}.
\end{align}