The Systems Biology field presents several complex combinatorial problems that can be in part reduced to an instance of the Hitting Set Problem (HSP), which is NP-Hard. These reduced problems often come with a large amount of data that needs to be processed, such as gene expression profiles, resulting in prohibitive computational costs for finding the exact solutions. There are some proposals to obtain exact solutions for HSP, including an approach which uses GPUs. However, such an approach is not scalable for real input sizes (thousands of variables). We propose a novel algorithm for solving HSP instances with thousands of variables by using: (i) clause sorting, which enables the efficient discarding of non-solution candidates, (ii) parallel generation and evaluation of candidate solutions through the use of GPUs, and (iii) support for multiple GPUs. To permit the execution on heterogeneous clusters, we determine the minimum kernel size that does not incur extra overhead and distribute tasks among available GPUs on demand. Our experimental results show that the combination of these techniques results in a speedup of 118.5, when using eight NVIDIA Tesla K20c in comparison with a ten-core Intel Xeon E5-2690 processor. Consequently, our algorithm can enable the usage of exact algorithms for solving the Hitting Set problem and applying it to real world problems.