An integral involving many exponential terms with quadratic exponents in the denominator

Question

Given $k$ points $\{p_1,\cdots, p_k\}$ in $\mathbb{R}^n$ and positive constants $r_1, ..., r_k$ and another positive constant $\alpha>0$. Is there a way to compute/approximate the following (Riemann) integral in terms of $\{p_1, ..., p_k\}$, $\{r_1, ..., r_k\}$ and $\alpha$? $$ \int_{\mathbb{R}^n} \frac{dx}{1+\alpha(e^{r_1\lVert x - p_1\rVert^2}+\cdots +e^{r_k\lVert x-p_k\rVert^2})} $$ where $\lVert v\rVert$ denotes the standard $2$-norm in $\mathbb{R}^n$. For example, in $\mathbb{R}^2$, I am considering integrals like $$ \int_{\mathbb{R}^2}\frac{dxdy}{1+2(e^{2\lVert (x,y)-(3,4)\rVert^2}+e^{5\lVert(x,y)-(6,7)\rVert^2})}. $$ By "approximating the integral", I mean a sequence of expressions $\{f_i(p_1, ..., p_k, r_1, ..., r_k, \alpha)\}_{i = 1}^\infty$ in terms of $\{p_1, ..., p_k\}$, $\{r_1, ..., r_k\}$ and $\alpha$, such that $f_i$ converges to the value of this integral as $i\to\infty$. Any idea or useful reference to related topics is welcome. Thanks