A variety of complex systems exhibits different types of relationships simultaneously that can be modeled by multiplex networks. A typical problem is to determine the community structure of such systems that, in general, depend on one or more parameters to be tuned. In this study we propose one measure, grounded on information theory, to find the optimal value of the relax rate characterizing Multiplex Infomap, the generalization of the Infomap algorithm to the realm of multilayer networks. We evaluate our methodology on synthetic networks, to show that the most representative community structure can be reliably identified when the most appropriate relax rate is used. Capitalizing on these results, we use this measure to identify the most reliable meso-scale functional organization in the human protein-protein interaction multiplex network and compare the observed clusters against a collection of independently annotated gene sets from the Molecular Signatures Database (MSigDB). Our analysis reveals that modules obtained with the optimal value of the relax rate are biologically significant and, remarkably, with higher functional content than the ones obtained from the aggregate representation of the human proteome. Our framework allows us to characterize the meso-scale structure of those multilayer systems whose layers are not explicitly interconnected each other-as in the case of edge-colored models-the ones describing most biological networks, from proteomes to connectomes.