Lemma 1. The inequality holds when one of variables is $0$.
We will prove this lemma later, with the similar argument of the inequality itself.
Lemma 2. The inequality holds when one of variables is at least $3$.
Without loss of generality, let $x\ge3$. It is easy to see $\left(3y+\frac16\right)^2>y^2+y$, therefore$$\sqrt{x^2+x}+\sqrt{y^2+y}+\sqrt{z^2+z}+\sqrt{w^2+w}< x+\frac12+3y+\frac16+3z+\frac16+3w+\frac16\\\le x+xy+xz+xw+1\le\left(1+x\right)\left(1+y\right)\left(1+z\right)\left(1+w\right)$$
Now, let $f\left(x,y,z,w\right)=\left(1+x\right)\left(1+y\right)\left(1+z\right)\left(1+w\right)-\sum_{cyc}\sqrt{x^2+x}$. When we want to find global minimum of multivatiable function, we need to find critical points and boundary points. By lemma 1 and 2, it is enough to show for critical points. Now, notice that$$\frac{d}{dx}f\left(x,y,z,w\right)=0\iff\left(1+y\right)\left(1+z\right)\left(1+w\right)=\frac{2x+1}{2\sqrt{x^2+x}}\\\iff\left(1+x\right)\left(1+y\right)\left(1+z\right)\left(1+w\right)=\frac{\left(x+1\right)\left(2x+1\right)}{2\sqrt{x^2+x}}$$
Let $g\left(x\right)=\frac{\left(x+1\right)\left(2x+1\right)}{2\sqrt{x^2+x}}$. At critical points of $f$, $g\left(x\right)=g\left(y\right)=g\left(z\right)=g\left(w\right).$ Since $g''\left(x\right)=\frac{\left(x + 1\right) \left(2 x + 3\right)}{8 \left(x^2+x\right)^{5/2}}>0
$, there are at most two solutions for $g\left(x\right)=c$ for any positive $c$. Therefore, we may assume $x=z=w$ or $x=z\land y=w$.
Case 1. $x=z=w$, prove $3\sqrt{x^2+x}+\sqrt{y^2+y}\le\left(x+1\right)^3\left(y+1\right)$
If $y\ge\frac12$, then $3\sqrt{x^2+x}+\sqrt{y^2+y}\le3\sqrt{x^2+x}+y+\frac12$, therefore we need to show $\left(x+1\right)^3\left(y+1\right)-3\sqrt{x^2+x}-y-\frac12\ge0$. Since it is increasing linear function of $y$, it is enough to show for $y=\frac12$, in other words, $\frac32\left(x+1\right)^3-3\sqrt{x^2+x}-1\ge0$. Notice that from $$9 x^6 + 54 x^5 + 135 x^4 + 168 x^3 + 63 x^2 - 18 x + 1\ge0$$ for all positive $x$, $\left(3\left(x+1\right)^3-2\right)^2\ge36\left(x^2+x\right)$, therefore our claim is proved.
If $y\le\frac12$, then notice that $\left(x^2+x+\frac14\right)^2-\left(x^2+x\right)=\left(x^2+x-\frac14\right)^2\ge0$. We need to show $\left(x+1\right)^3\left(y+1\right)-3\sqrt{x^2+x}-y^2-y-\frac14\ge0$. Since it is concave function of $y$, it is enough to prove for extreme values of $y$, which is $0$ and $\frac12$. For $y=\frac12$, it is $\frac32\left(x+1\right)^3-3\sqrt{x^2+x}-1\ge0$ which we already proved. For $y=0$, it is $\left(x+1\right)^3-3\sqrt{x^2+x}-\frac14\ge0$. Again from $x^2+x+\frac14\ge\sqrt{x^2+x}$, $$\left(x+1\right)^3-3\sqrt{x^2+x}-\frac14\ge\left(x+1\right)^3-3\left(x^2+x+\frac14\right)-\frac14=x^3\ge0$$
Case 2. $x=z$ and $y=w$, prove $2\sqrt{x^2+x}+2\sqrt{y^2+y}\le\left(x+1\right)^2\left(y+1\right)^2$
Notice that$$\left(2x+1\right)^2\left(2y+1\right)^2=\left(2\sqrt{x^2+x}+2\sqrt{y^2+y}\right)^2+\left(4\sqrt{x^2+x}\sqrt{y^2+y}-1\right)^2$$Therefore $2\sqrt{x^2+x}+2\sqrt{y^2+y}\le\left(2x+1\right)\left(2y+1\right)<\left(x+1\right)^2\left(y+1\right)^2$.
Now we checked all critical points of $f$, so we finish the proof.
Proof of lemma 1. We want to find global minimum of $f_1\left(x,y,z\right)=f\left(x,y,z,0\right)$. We already checked boundary points where $x,y,z$ is sufficiently large. For one of variable is $0$, WLOG assume $z=0$. $$\sqrt{x^2+x}+\sqrt{y^2+y}\le x+\frac12+y+\frac12\le\left(x+1\right)\left(y+1\right)$$Therefore, we only need to check critical points of $f_1$. By the same argument as main proof, at critical points of $f_1$, $g\left(x\right)=g\left(y\right)=g\left(z\right)$ and at least two variables are same. We can assume $x=z$.
Now we need to show $2\sqrt{x^2+x}+\sqrt{y^2+y}\le\left(x+1\right)^2\left(y+1\right)$. Indeed,$$\left(2x+1\right)^2\left(y+1\right)^2=\left(2\sqrt{x^2+x}+\sqrt{y^2+y}\right)^2+\left(2\sqrt{x^2+x}\sqrt{y^2+y}-1\right)^2+4xy\left(x+1\right)+y+1$$ therefore $2\sqrt{x^2+x}+\sqrt{y^2+y}\le\left(2x+1\right)\left(y+1\right)\le\left(x+1\right)^2\left(y+1\right)$.