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}