**Abstract** : 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 series in powers and logarithms of the parameters involved in the problem under consideration. However, in most of the cases, several series representations exist for a given integral. They 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: (a) different sets of poles which correspond to different convergent double series representations of a given integral, (b) the regions of convergence of all these series (without an a priori full knowledge of their general term), and (c) 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 N = 4 SYM theory. Mellin-Barnes integrals of higher dimension are also considered.