Separation implies summation
When we perform a separation of variables of this form,$$\rho(x, t)= X(x) T(t),
\\
X\dot T = DX''T,\\
\dot T/T = D X''/X = -Dk^2 \text{ for some constant } k\\
\rho(x,t) = A \cos(kx)~e^{-Dk^2 t} + B \sin(kx)~e^{-Dk^2 t} \text{ for arbitrary } A,B,$$it is important to remember that we started from a really narrow starting place with $\rho$ and that this is not a fully general solution. The fully general solution requires, as you have said, taking a sum of these. In this case you can use the linearity of the differential equation $\dot \rho = D \rho''$ to argue that if you can determine that $$\rho(x,0) = \int_{-\infty}^{\infty}dk~\cos(kx)~f(k)$$ for some $f(k)$ then that is a sum of coefficients times the functions you got above, and the solution is then $\int dk~e^{-Dk^2t}~\cos(kx)~f(k).$
In other words: the boundary condition is the full function $\rho(x, 0)$. It's a 2D space consisting of half of $\mathbb R^2$ and your boundary condition fixes the function on that boundary $t=0.$
Fourier transforms to the rescue
When the boundary conditions are not either finite-range bounded or else periodic, the cosine-and-sine choice above proves somewhat awkward and we will find it easier to just complexify:
$$ \rho(x, t) = \alpha e^{ikx}~e^{-D k^2 t} + \beta e^{-ikx} e^{-D k^2 t}.$$ We can then see that the two cases are really overlapping and once we take a sum over all $k$ we do not need both of them because we will get both positive and negative $k$, the separated solution can be simplified to just $\alpha e^{ikx} e^{-Dk^2 t}.$ Which is good because the equivalent integral there is known as the inverse Fourier transform, $$\rho(x, 0) = \frac{1}{2\pi}\int_{-\infty}^{\infty} e^{ikx}~dk~\int_{-\infty}^{\infty} e^{-ikx'}dx'~\rho(x', 0),$$
and therefore we have that by this prescription $$
\rho(x,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty} e^{ikx} dk~e^{-Dk^2t}\int_{-\infty}^{\infty} e^{-ikx'} dx' \rho(x', 0).$$
But the boundary condition is a $\delta$-function
You have had some trouble specifying $\rho(x', 0)$ to be a probability density expressing certainty that a particle is at the origin. This is because to do that properly we need a function that is impossible (but surprisingly well-behaved!) called the Dirac $\delta$-function. It is a ‘function’ $\delta(x)$ such that $$\int_a^b dx~\delta(x) f(x) = \begin{cases}f(0)&\text{if } a < 0 < b,\\
-f(0)&\text{if } a > 0 > b,\\
0&\text{otherwise}.
\end{cases}$$
One typically thinks of it as an infinitely thin, infinitely large spike normalized so that the area under the curve is $1$, but other constructions also involve functions that oscillate wildly at $x\ne 0$ so that they cannot create nonzero integrals when multiplied by a standard smooth function. Part of this idea does mean that there is a great ambiguity about what happens if $a=0$ or $b=0$ in the integral above, too. But it can be realized as the limit of smooth functions e.g. a Gaussian with finite area 1 pushed to be thinner and thinner, so one can imagine that it itself is a smooth function too by just "backing off" from that limit by the tiniest bit.
Anyway in probability densities, the area under a probability density curve is the actual probability, and the Dirac $\delta$-function lets us define a probability density for "I know with probability $p$ that it is at point $x_0$" as the density contribution $p~\delta(x - x_0).$ In this case therefore you want to consider $\rho(x, 0) = \delta(x),$ and using the definition we immediately find that the Fourier transform of $\rho(x, 0)$ was just $e^{-ik0} = 1$,
$$
\rho(x,t)=\int_{-\infty}^{\infty}\frac{dk}{2\pi}~ e^{ikx-Dk^2t}.$$
Complete the square, complete the derivation
So now we want to just complete the square, because we know the Gaussian integral that $\int_{-\infty}^\infty du~e^{-u^2} = \sqrt{\pi}.$
$$D~t~k^2 - ix~k =D~t~\left(k - \frac{ix}{2Dt}\right)^2 + \frac{x^2}{4Dt}.$$
Performing the Gaussian integral just gives a constant and the solution then comes out as $$\rho(x, t) = \frac{e^{-x^2/(4Dt)}}{\sqrt{4\pi Dt}}, $$a normal distribution with standard deviation $\sqrt{2Dt}.$
Now, since I have just shown you $\delta$-functions, I would like to stretch you one stage further, to tell you that actually linearity lets you use this result, instead of the separation-of-variables result, to see the evolution of an arbitrary system and alleviate your concerns about negative probabilities.
This solution that we came up with is special precisely because its boundary condition is a $\delta$-function. We call this result a “Green's function solution” of the equation for that reason. Every other function $\rho_0(x)$ could be expressed as $$\rho_0(x) = \int_{-\infty}^\infty dx'~\rho_0(x')~\delta(x - x'),$$ per the definition of $\delta$ above. As a consequence under the heat equation that function evolves over time $t$ to $$\rho(x, t) = \int_{-\infty}^{\infty} dx'~\rho_0(x')~\frac{e^{-(x-x')^2/(4Dt)}}{\sqrt{4\pi D t}},$$ from which it is easier to see that if $\rho_0(x)$ is nonnegative and integrates to $1$ then $\rho(x, t)$ for any fixed $t$ will also be nonnegative and integrate to $1$.