Consider the heat equation on real line: $$ \begin{matrix} \frac{\partial u}{\partial t} = c^2 \frac{\partial^2 u}{\partial x^2} \\ u(x,0) = h(x) \\ x \in \mathbb{R}, t \geq 0 \end{matrix} $$
The solution is: $$ u(x,t) = \frac{1}{\sqrt{\pi}} \int_{-\infty}^{\infty} e^{-r^2} h(x+2cr\sqrt{t}) dr $$
On the other hand, let $Y_{x,t}$ be the Gaussian process $ Y_{x,t} \sim \mathcal{N} (x, 2c^2 t) $, then the solution is also given by the expectation: $$ u(x,t) = E[h(Y_{x,t})] $$
Why is that? Of course it is related to the Gaussian heat kernels being fundamental solutions of the heat equation, but I am looking for deeper explanations, presumably probabilistic and physical explanations. For example, does it have anything to do with the Brownian motion of particles on real line, whose temperature distribution is described by the heat equation?
Of course, asking the question
why?is always a bit dangerous in math. I can prove a theorem, but does that explain why that theorem should be true? Anyway, for your question, there indeed is a deeper answer. Read the page on the Fokker-Planck equation, https://en.wikipedia.org/wiki/Fokker%E2%80%93Planck_equation. This shows that for a very wide range of stochastic processes, the distribution can be computed by solving a PDE. Maybe diving into the derivation of the Fokker-Planck equation can answer your question of `why'.