Given a wavefunction, $\psi(x,0)$, the corresponding quantity, $\psi(x,t)$, is obtained from the integral:
$$
\boxed{\psi(x,t) = \int_{-\infty}^{\infty}dy\,Z[x,t;y,0]\psi(y,0)\,}
$$
where $Z[x,t;y,0]$ is a path integral, $\int_{x,0}^{y,t} \mathcal{D}x\,e^{\frac{i}{\hbar}I[x]}$, and $I[x]$ the single particle action of interest. For free non-relativistic particles, $I[x]=\int_0^td\tau\frac{m}{2}(\frac{dx}{d\tau})^2$, and the correctly normalised path integral reads:
$$
Z[x,t;y,0] = \sqrt{\frac{m}{2\pi i\hbar t}}\,\exp\Big(i\frac{m(x-y)^2}{2\hbar t}\Big).
$$
The initial wave packet is of the form:
$$
\psi(y,0)=Ne^{\frac{i}{\hbar}ky}e^{-\frac{1}{2}b(y-x_0)^2},
$$
where you can read off the $b,N$ of interest from your formula, e.g. $b=1/(2\sigma_x^2)$.
The objective is to compute $\psi(x,t)$, which amounts to carrying out the $y$ integral in the above boxed formula. This is a Gaussian and so can easily be carried out.
You'll primarily need the result:
$$
\int_{-\infty}^{\infty} da\,e^{iza^2} = \sqrt{\frac{\pi i}{z}},\quad {\rm if }\quad {\rm Im}\,z\geq0,\quad z\neq0,
$$
which follows from Cauchy's theorem applied to $\oint_C da\,e^{iza^2}=0$, where the contour $C$ is that in the figure:
It is convenient to write $z=\rho e^{i\phi}$ (with $0\leq \phi<\pi$), and then given any allowed $\phi$ you choose $\theta$ such that $\phi+2\theta=\pi/2$. This enables you to take $R\rightarrow \infty$ while ensuring the angular contour contributions vanish, and you are left with $\int_{-\infty}^{\infty}da\,e^{iza^2}=\sqrt{i}e^{-i\phi/2}\int_{-\infty}^{\infty}da\,e^{-\rho a^2}$. Evaluating the remaining Gaussian integral and rewriting the result in terms of $z=\rho e^{i\theta}$ leads to the displayed result.
The next step is to redefine your integration variable, $a\rightarrow a'=a-b$ (with $b\in \mathbb{R}$),
$$
\int_{-\infty}^{\infty} da\,e^{iz(a+b)^2} = \sqrt{\frac{\pi i}{z}},
$$
which implies:
$$
\int_{-\infty}^{\infty} da\,e^{iza^2+2izba} = \sqrt{\frac{\pi i}{z}}\,e^{-izb^2},\quad {\rm if }\quad {\rm Im}\,z\geq0,\quad z\neq0.
$$
We are almost there. The only thing remaining to do is to notice that both the left- and right-hand sides of this last expression are analytic in $b$, so we can extend $b$ into the complex plane, $b\rightarrow w=b+ib'$. Finally, we redefine $w\rightarrow zw$. Therefore, for generic complex numbers $z$ (with ${\rm Im}\,z\geq0$, $z\neq0$) and $w$,
$$
\boxed{\,\int_{-\infty}^{\infty} da\,e^{iza^2+2iwa} = \sqrt{\frac{\pi i}{z}}\,e^{-iw^2/z},\quad {\rm if }\quad z,w\in \mathbb{C},\quad {\rm and}\quad{\rm Im}\,z\geq0,\quad z\neq0\,\,\,}
$$
This is all you need in order to evaluate the above integral that gives you $\psi(x,t)$. You just read off $z,w$ and $a$ as appropriate for your case of interest and apply this result. (Of course the same formula applies for general $x_0$.) The result is:
$$
\boxed{\psi(x,t) = \frac{1}{\sqrt{1+i\frac{\hbar b }{m}t}}\exp\Bigg\{\!\!-\frac{i}{\hbar}\frac{k^2t}{2m}\frac{\big[1+i\frac{\hbar b}{k}(x-x_0)\big]^2}{1+i\frac{\hbar b}{m}\,t}\Bigg\}\psi(x,0)\,\,}
$$
Recall that $
\psi(x,0)=Ne^{\frac{i}{\hbar}kx}e^{-\frac{1}{2}b(x-x_0)^2}.
$
