This is a well known result from diffusion PDE theory. For example [1] studied the semi-linear equation in form of
$$u_t - \Delta u + u^\gamma = 0, \quad (t,x) \in (0,T) \times \mathbb{R}^n, \gamma>1$$ with initial value $u(0,x) = g(x), \ x \in \mathbb{R}^n.$
And with typical Green function representation we could write (2.1) in [1].
$$u(t)=e^{t\Delta}g+{\displaystyle \int_{0}^{t}e^{(t-s)\Delta}u^{\gamma}ds}$$ with $e^{t\Delta}g(x)={\displaystyle \int_{\mathbb{R}^{n}}Gr_{t}(x-y)dy}$ and $Gr_{t}(x)=\left(4\pi t\right)^{-\frac{n}{2}}e^{-\frac{\|x\|^{2}}{4t}}$. Therefore an intuition that such a solution $u$ exists if $u^\gamma,\gamma>1$satisfy regularities. And if $f(x,u,u_x)$ is somehow regular enough then we can expand it into analytic form in terms of powers of $u,u_t$ and the Green function representation argument applied.
A rigid proof requires more work but I believe that is what I(and possibly @ Michael Renardy) had in mind about regularity on solution to the equation above.
[1]Weissler, Fred B. "Existence and non-existence of global solutions for a semilinear heat equation." Israel Journal of Mathematics 38.1 (1981): 29-40.