I wonder if this can be solved using a Feynman-Kac type argument.
Let $X_{t,s}(a)$ be a path labeled at the time $t$ (i.e. $X_{t,t}(a)=a$) going back to $s$.  Let $X_{t,s}(a)$ satisfy the following stochastic equation:
$$
{\rm{d}}X_{t,s}(a) =\sqrt{2\sigma(u(x,s))}\ \hat{{\rm{d}}} W_s
$$
where $W_s$ is a 1-dimensional Wiener process and $\hat{{\rm{d}}}$ indicates that this is a backwards Ito differential $\sigma$ is some arbitrary smooth function of $u$.  Then, the following is true:

$\textbf{Claim:}$  A function $u(x,t)$ is a solution to the equation
$$
\partial_t u = \sigma(u) \triangle u
$$
if and only if the pair $(u,X)$ satisfies the following stochastic system:
$$
{\rm{d}}X_{t,s}(a) =\sqrt{2\sigma(u(x,s))}\ \hat{{\rm{d}}} W_s
$$
$$
u(x,t)= \mathbb{E}\left[u_0(X_{t,0}(a))\right]
$$
where the expectation is taken over Brownian motions.

Here I am assuming we are solving the heat equation on the real line with no boundary conditions but it is straight-forward to incorporate boundaries. Of course, you can choose $\sigma$ to be the function you desire (for example the one in your question). I hope this helps.