A variety of complex systems exhibit 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. ! arXiv:1801.10144v1 [q-bio.MN] 30 Jan 2018 • GN/ER. This system is generated by combining an Erdos-Renyi layer with one generated by the Girvan-Newman (GN) benchmark [9]. A community structure is only present in the GN layer, hence this multiplex network is used to test the impact of coupling noise to a structured population.• GN/GN. This system consists of two GN networks, with tunable cross-layer community overlap, generated as following. First, we create a single-layered GN network and duplicate it to generate a multiplex network with two layers. Then, in second layer, we iteratively swap the neighbours of pairs of randomly selected nodes to change the community structure. The swapping procedure is repeated until the ratio of community overlapping across layers, defined by the fraction of nodes that belong to the same communities across layers, is reached. Two classes of networks are generated, corresponding to a different amount of overlapping across layers (50% and 10%, respectively).• GN/CGN. This system consists of two layers, one GN network and one complementary GN (CGN). In order to explain how these two layers are generated,