As I understand, we take for granted that the solution of the heat equation is given in terms of the appropriate heat kernel. This heat kernel is given in terms of the Gauss–Weierstrass kernel $k_t(x) = (2 \pi t)^{-1/2} \exp(-x^2 / (2 t))$; namely, we have
$$ u(t, x) = \int_{-a}^a p(t, x, y) f(y) dy, $$
where
$$ p(t, x, y) = \sum_{n \in 2 \mathbb{Z}} k_t(y - x + 2 n a) - \sum_{n \in 2 \mathbb{Z} + 1} k_t(y + x + 2 n a) $$
for the Dirichlet boundary condition,
$$ p(t, x, y) = \sum_{n \in 2 \mathbb{Z}} k_t(y - x + 2 n a) + \sum_{n \in 2 \mathbb{Z} + 1} k_t(y + x + 2 n a) $$
for the Neumann boundary condition, and
$$ p(t, x, y) = \sum_{n \in \mathbb{Z}} k_t(y - x + 2 n a) $$
for the periodic boundary condition. The question is: where these kernels are shown to be transition probabilities of the killed, reflected and periodized Brownian motion in $(-a, a)$.
Killed Brownian motion in an interval is discussed, for example, by Karatzas and Shreve, Brownian Motion and Stochastic Calculus. In Section 4.3 killed process in half-line is studied in detail, and killing in an interval is left as an exercise (with an extensive hint); see p. 267.
For extensions to arbitrary domains, see, for example, Theorem 2.13 in Chung and Zhao, From Brownian Motion to Schrodinger's Equation.
The periodic Brownian motion is typically defined as $f(B_t)$, where $B_t$ is the usual Brownian motion on $\mathbb{R}$ and $f(x) = x - 2 a \lfloor(x - a) / (2 a)\rfloor$ is a piecewise linear function. In this case the expression for the heat kernel is a straightforward consequence of the definition. It is quite likely written explicitly in some reference, but I am not aware of any.
The reflected Brownian motion can be defined in a number of equivalent ways:
The easiest way is to "fold" the periodic Brownian motion. In other words, define the reflected Brownian motion as $f(B_t)$, where $B_t$ is the usual Brownian motion and $f(x) = \tfrac{2 a}{\pi} \arcsin(\sin(\tfrac{\pi x}{2 a}))$. Then, again, the expression for the heat kernel is a direct consequence of the definition.
Another way is to define the reflected Brownian motion as a solution to a Skorokhod problem. This is equivalent to definition 1 as a consequence of Tanaka's formula. Any reference on Tanaka's formula will discuss the case of the half-line. General convex domains in $\mathbb{R}^n$ and general diffusions were already considered by Tanaka. I am not aware of a reference that would discuss the case of the Brownian motion in an interval, but I am not an expert.
Feller's work on diffusion processes in an interval defined the reflected diffusion by an appropriate limiting procedure, and linked it directly to an appropriate differential operator. Again, this is not "explicit" and likely too general.
Finally, one can define the reflected Brownian motion as a Markov process linked to the usual energy form $\int_D |\nabla u(x)|^2 dx$ for $u$ in the Sobolev space $W^{1,2}(D)$. In this case the process corresponds to the Laplace operator with Neumann boundary conditions by the very definition.
Finally, let me mention that Pitman and Yor provide an extensive list of references. All of them seem to deal with rather general questions, but perhaps some discuss the interval as an example.
