Multiple Mellin-Barnes integrals are often used for perturbative calculations in particle physics. In this context, the evaluation of such objects may be performed through residues calculations which lead to their expression as multiple power series and logarithms of the parameters involved in the problem under consideration. However, in most of the cases, several convergent series representations exist for a given integral. These series converge in different regions of values of the parameters, and it is not obvious to obtain them. For twofold integrals, we present a method which allows to derive straightforwardly and systematically: First, different sets of poles which correspond to different convergent double series representations of a given integral, second the regions of convergence of all these series (without an a priori full knowledge of their general term), and third the general term of each series (this may be performed, if necessary, once the relevant domain of convergence has been found). This systematic procedure is illustrated with some integrals which appear, among others, in the calculation of the two-loop hexagon Wilson loop in $\mathcal {N} = 4$N=4 SYM theory. Mellin-Barnes integrals of higher dimension are also considered.