I am trying to solve this double integration, but not getting it exactly.
$$P_{int} = \text{Pr}(X_{is}<X_{ie})$$
$$P_{int} = \int \int_{x_{is}<{x_{ie}}}f_{X_{ie}}(x_{ie})f_{X_{is}}(x_{is})dx_{ie}dx_{is}\tag1$$
$$P_{int} = \int_{0}^{\infty}\Gamma\left(m_i,\frac{m_i}{\sigma^2_{is}}x_{ie}\right)f_{X_{ie}}(x_{ie})dx_{ie}\tag2$$
where: $$f_{X_{is}}(x_{is}) = \frac{1}{\Gamma(m_i)}\left(\frac{m_i}{\sigma^2_{is}}\right)^{m_i}x_{is}^{m_i-1}\exp\left(-\frac{m_ix_{is}}{\sigma^2_{is}}\right)$$ and $$f_{X_{ie}}(x_{ie}) = \frac{1}{\Gamma(k_i)}(\frac{k_i}{\sigma^2_{ie}})^{k_i}x_{ie}^{k_i-1}\exp\left(-\frac{k_ix_{ie}}{\sigma^2_{ie}}\right)$$ are PDF of $X_{is},X_{ie}$ respectively. And $\Gamma(m_i,\frac{m_i}{\sigma^2_{is}}x_{ie})$ is the lower incomplete Gamma function and is defined as $$\Gamma \left(m_i,\frac{m_i}{\sigma ^2_{is}}x_{ie}\right) = \int_0^{\frac{m_ix_{ie}}{\sigma^2_{is}}}\frac{1}{\Gamma(m_i)}x^{m_i-1}\text{exp}(-x)dx$$
My doubt is how eq. $(2)$ is obtained from eq. $(1)$.
Any help in this regard will be highly appreciated.
They are known densities: $X\sim \text{Gamma}\left(m;\frac{m}{\sigma_X^2} \right)$ and $Y\sim \text{Gamma}\left(k;\frac{k}{\sigma_Y^2} \right)$ where $\frac{m}{\sigma_X^2}$ is the rate parameter.
To simplify the notation I used $X,Y$ instead of your rv and eliminated some indexes not needed to understand the question.
To understand the passage simply observe that, assuming independence,
$$\mathbb{P}[X<Y]=\int_0^{\infty}f_Y(y)\left[\underbrace{\int_0^y f_X(x)dx }_{F_X(y)} \right]dy=\int_0^{\infty}f_Y(y)F_X(y)dy$$
and $F_X(y)$, the CDF of a gamma function can be expressed in terms of the lower incomplete gamma function
Further details about your comment
the integral region is the following
Thus, you have to integrate the joint density $f_{XY}(x,y)=f_X(x)f_Y(y)$ in that regiom. This can be done in two equivalent ways
$$\mathbb{P}[X<Y]=\int_0^{\infty}f_X(x)\left[\int_x^{\infty}f_Y(y)dy \right]dx=\int_0^{\infty}f_Y(y)\left[\int_0^{y}f_X(x)dx \right]dy$$
They chose the former.