Solution 1: Let $p=x y z$ and $q=(x+1)(y+1)(z+1)$. Then, we get $$p q=[x(1+y)] \cdot[y(1+z)] \cdot[z(1+x)]=(1-p)(2-p)(4-p)$$ Additionally, note that $$q-p=x y+y z+z x+x+y+z+1=(x+x y)+(y+y z)+(z+x z)+1=8-3 p$$ Therefore, we have $q=8-2 p$. Substituting this into our earlier equation gives us $$p(8-2 p)=(1-p)(2-p)(4-p)$$ We can rearrange this to get $(4-p)\left(2-5 p+p^{2}\right)=0$. Solving this gives us $p=4, \frac{5 \pm \sqrt{17}}{2}$. Thus, our maximum solution is $\frac{5+\sqrt{17}}{2}$, which yields an answer of 5272. To show that such a solution exists, see Solution 2. Solution 2: Let $r=x y z-1$. Observe that $$r x=x(y+y z+x y z)-(x+x y+x y z)=2 x-1 \Longleftrightarrow x=\frac{1}{2-r}$$ Similarly, $y=\frac{2}{4-r}$ and $z=\frac{4}{1-r}$. Therefore $8=(1+r)(1-r)(2-r)(4-r)$. This factors as $r(r-3)\left(r^{2}-3 r-2\right)=0$, so the maximum possible value for $r$ is $\frac{3+\sqrt{17}}{2}$ Now let's check that this yields a valid solution for $x, y, z$. Let $r=\frac{3+\sqrt{17}}{2}$ and let $x=\frac{1}{2-r}, y=\frac{2}{4-r}, z=$ $\frac{4}{1-r}$. Then $x y z-1=\frac{8}{(2-r)(4-r)(1-r)}-1=1+r-1=r$. Now, we may do our above computations in reverse to get $$2 x-1=2 x-(2-r) x=r x=x^{2} y z-x=x(y+y z+x y z)-(x-x y+x y z)$$ Repeating the same thing for $y$ and $z$ yields that $$\left(\begin{array}{ccc} -1 & x & 0 \\ 0 & -1 & y \\ z & 0 & -1 \end{array}\right)\left(\begin{array}{l} 1 \\ 2 \\ 4 \end{array}\right)=\left(\begin{array}{ccc} -1 & x & 0 \\ 0 & -1 & y \\ z & 0 & -1 \end{array}\right)\left(\begin{array}{l} x+x y+x y z \\ y+y z+x y z \\ z+x z+x y z \end{array}\right)$$ However, since $x y z-1 \neq 0$, the determinant of the matrix is nonzero, so we may multiply by its inverse to find that $$\left(\begin{array}{l} 1 \\ 2 \\ 4 \end{array}\right)=\left(\begin{array}{l} x+x y+x y z \\ y+y z+x y z \\ z+x z+x y z \end{array}\right)$$ Therefore this construction is valid.