Multiplex decomposition of non-Markovian dynamics and the hidden layer reconstruction problem

TL;DR

This paper introduces a method to determine whether complex system interactions are best modeled as a multiplex network or a single layer, using local information from a random walk, and extends this to reconstruct non-Markovian dynamics.

Contribution

It provides a robust technique to identify hidden multiplex structures and decompose complex dynamics into Markov switching diffusive modes, with proven theorems and real-world applications.

Findings

Successfully reconstructs hidden multiplex network architecture.

Decomposes non-Markovian dynamics into Markov switching diffusive modes.

Proves the possibility of exact dynamics reconstruction for finite-order processes.

Abstract

Elements composing complex systems usually interact in several different ways and as such the interaction architecture is well modelled by a multiplex network. However often this architecture is hidden, as one usually only has experimental access to an aggregated projection. A fundamental challenge is thus to determine whether the hidden underlying architecture of complex systems is better modelled as a single interaction layer or results from the aggregation and interplay of multiple layers. Here we show that using local information provided by a random walker navigating the aggregated network one can decide in a robust way if the underlying structure is a multiplex or not and, in the former case, to determine the most probable number of hidden layers. As a byproduct, we show that the mathematical formalism also provides a principled solution for the optimal decomposition and projection of complex, non-Markovian dynamics into a Markov switching combination of diffusive modes. We validate the proposed methodology with numerical simulations of both (i) random walks navigating hidden multiplex networks (thereby reconstructing the true hidden architecture) and (ii) Markovian and non-Markovian continuous stochastic processes (thereby reconstructing an effective multiplex decomposition where each layer accounts for a different diffusive mode). We also state and prove two existence theorems guaranteeing that an exact reconstruction of the dynamics in terms of these hidden jump-Markov models is always possible for arbitrary finite-order Markovian and fully non-Markovian processes. Finally, we showcase the applicability of the method to experimental recordings from (i) the mobility dynamics of human players in an online multiplayer game and (ii) the dynamics of RNA polymerases at the single-molecule level.

Figures

Equations

$$\mathcal{D}(m):=\frac{1}{m}\sum _{\mathcal{B}(m)}{P}_X(Z_1,\dots ,Z_m)\log \frac{P_X(Z_1,\dots ,Z_m)}{P_Y(Z_1,\dots ,Z_m)},$$

$$P(\{X(t){\}}_{t=1}^N|L)=\int P(\{X(t){\}}_{t=1}^N|{\mathbf{T}}_L,{\mathbf{R}}_L)\times$$

$$\times p_{0,T}({\mathbf{T}}_L)p_{0,R}({\mathbf{R}}_L)\mu (\mathsf{d}{\mathbf{T}}_L\times \mathsf{d}{\mathbf{R}}_L),$$

$${\hat{L}}^{\text{MAP}}=\arg \underset{L\in \{1,2,\dots ,L_{+}\}}{\max }P(\{X(t){\}}_{t=1}^N|L)P_0(L),$$

$$\tilde{X}(t+1)=\{\begin{array}{l}\tilde{X}(t)+1\mathrm{mod}K,\text{if}X(t+1)>X(t)\\ \tilde{X}(t)-1\mathrm{mod}K,\text{if}X(t+1)<X(t),\end{array}$$

$$\dot{x}=-x+\xi ,$$

$$\dot{x}=-x+\eta ,$$

$$\dot{x}$$

$$\dot{y}$$

$$P_{\mathbb{K}}^h(i_t|i_{t-h:t-1}):=P(X(t)=i_t|X(t-h:t-1)=i_{t-h:t-1})$$

$$X\stackrel{d}{=}\underset{h\to \infty }{\lim }{X}_h$$

$$\underset{t\to \infty }{\lim }\underset{i_{0:t}\in {\mathbb{K}}^{t+1}}{\sup }P(X(0:t)=i_{0:t})=0.$$

$$\underset{i_{0:t}\in {\mathbb{K}}^{t+1}}{\sup }$$

$$\mathcal{M}:\mathcal{B}([0,1))\times [0,1)\to [0,1],$$

$$\mathbf{T}(y)=\left[\begin{array}{ccc}{T}_{11}(y)& \dots & T_{iK}(y)\\ ⋮& \ddots & ⋮\\ T_{K1}(y)& \dots & T_{KK}(y)\end{array}\right],y\in [0,1),$$

$$P(X(t)=i_t|X(0:t-1)=i_{0:t-1})=$$

$$=P_{i_0}\int \cdots \int \prod _{k=1}^{t-1}{T}_{i_{k+1}{i}_k}(y_k)\mathcal{M}(d{y}_k|y_{k-1}){\mathcal{M}}_0(d{y}_0),$$

$${\text{T}}^{(1)}=\left[\begin{array}{cccc}0& 1/3& 1/3& 1/3\\ 1/3& 0& 1/3& 1/3\\ 1/3& 1/3& 0& 1/3\\ 1/3& 1/3& 1/3& 0\end{array}\right]$$

$${\text{T}}^{(2)}=\left[\begin{array}{cccc}0& 1/3+2ϵ& 1/3-ϵ& 1/3-ϵ\\ 1/3+2ϵ& 0& 1/3-ϵ& 1/3-ϵ\\ 1/3+2ϵ& 1/3-ϵ& 0& 1/3-ϵ\\ 1/3+2ϵ& 1/3-ϵ& 1/3-ϵ& 0\end{array}\right]$$

$${\text{T}}^{(1)}=\left[\begin{array}{cccc}0& 1/3& 0& 2/3\\ 2/3& 0& 1/3& 0\\ 0& 2/3& 0& 1/3\\ 1/3& 0& 2/3& 0\end{array}\right]$$

$${\text{T}}^{(2)}=\left[\begin{array}{cccc}0& 1/3-ϵ& 2ϵ& 2/3-ϵ\\ 2/3+ϵ& 0& 1/3-ϵ& 0\\ 2ϵ& 2/3-3ϵ& 0& 1/3+ϵ\\ 1/3+ϵ& 0& 2/3-ϵ& 0\end{array}\right]$$

$${\text{T}}^{(1)}=\left[\begin{array}{cccc}0& 1/3& 0& 2/3\\ 2/3& 0& 1/3& 0\\ 0& 2/3& 0& 1/3\\ 1/3& 0& 2/3& 0\end{array}\right]$$

$${\text{T}}^{(2)}=\left[\begin{array}{cccc}0& 1/3& 1/3& 1/3\\ 1/3& 0& 1/3& 1/3\\ 1/3& 1/3& 0& 1/3\\ 1/3& 1/3& 1/3& 0\end{array}\right]$$

$${\text{T}}^{(1)}={\text{T}}^{(2)}=\left[\begin{array}{cccc}0& 1/2& 0& 1/2\\ 1/2& 0& 1/2& 0\\ 0& 1/2& 0& 1/2\\ 1/2& 0& 1/2& 0\end{array}\right]$$

$$\mathbf{X}(t)=\{{\mathbf{R}}_L,{\mathcal{T}}_L,X(t),\ell (t)\},t=0,1,2,...$$

$$P(Y(t)=j|X(t)=i)=\left\{\begin{array}[]{ll}1,&\mbox{if $j=i$,}\\ 0,&\mbox{otherwise,}\\ \end{array}\right.$$

$$Y(1:N)=\{Y(1),Y(2),\dots ,Y(N)\}\in \{0,1,\dots ,K-1{\}}^N.$$

$$P\left(Y(1:N)|L\right)=\int P\left(Y(1:N)|{\mathcal{T}}_L,{\mathbf{R}}_L\right)p_0({\mathcal{T}}_L)p_0({\mathbf{R}}_L)\mu (\mathsf{d}{\mathcal{T}}_L\times \mathsf{d}{\mathbf{R}}_L),$$

$$P\left(L|Y(1:N)\right)\propto P\left(Y(1:N)|L\right)P_0(L)$$

$${\hat{L}}_{MAP}$$

$$P(L|\mathcal{O})\propto P(\mathcal{O}|L)P_0(L)$$

Full text

I Introduction

Network Science has emerged as a powerful unifying framework for studying the emergence of collective phenomena in real complex systems from different domains newmanbook ; barabook , and has allowed to increase the accuracy and predictive power of minimal models of complex dynamics, including epidemic spreading vespignaniRMP , synchronisation arenasphysrep , or social dynamics castellanoRMP . One of the most fascinating challenges faced in the last few years by Network Science is the need to incorporate and couple several network structures in order to correctly capture the inherently multidimensional nature of interaction patterns in real-world systems. As a result, much effort has been recently devoted to the definition and study of multilayer and multiplex networks physrep ; jcn ; battiston2016review . The ubiquity of such structures in social, biological and technological systems has required the revision of the several canonical dynamical models that were previously studied only on isolated complex networks, including percolation buldy ; Gao012 ; grass ; BiD014 ; radi , diffusion dynamics gomez13 ; radi2 , navigation Manlio14 ; Sole16 ; Battiston2016RW , epidemics mendiola ; granell13 ; Buono14 ; Sanz14 , evolutionary games GG12 ; Wang14 ; Mata15 , synchronization DelGenio16 , or opinion dynamics Diakonova2014 ; Diakonova2016 . Notably the collective behaviour of such complex systems depends strongly on whether they can described by isolated networks or by coupled networks NatPhys16 , highlighting the importance of the multilayer architecture of real-world systems.

Multilayer network models of real-world systems face two fundamental and dual challenges. The first one is the necessity to assess in a systematic way whether a multilayer network model is adequate to represent the system, and when such model gives redundant information. This challenge was first addressed in nicosia , and constitutes nowadays an intense field of research. The dual challenge aims at understanding whether an empirical network whose multilayer character is not directly observable is genuinely monolayer or is only an aggregated projection of a hidden multilayer network (see Fig. 1b for an illustration of a multiplex network with $L=2$ layers and its aggregated projection). Such scenario has received much less attention despite being, for instance, central for networks arising in natural systems whose architecture is not directly observable, as in genetic networks or in brain functional networks where pairs of nodes modelling different brain areas can interact according to an a priori unknown range of different biological pathways Zanin .

In this article we provide a method to identify the hidden multi-layer structure of a complex system from coarse-grained dynamical measurements of its state. We show that, by using only local information extracted from simple random-walk statistics, it is possible to discriminate whether the underlying structure of the system is actually a single-layer or a multi-layer network, and in the latter case, to estimate the number of interacting layers in the system. Note that methods to infer network topological properties via random walk statistics have been explored previously PNASRosvall ; tiago3 . Notably, our discrimination method exploits the breaking of Markovianity occuring in a coarse-grained multi-layer random walk, while the method to estimate the most probable number of layers is based on a maximum a-posteriori (MAP) probabilistic criterion which can be implemented via numerical integration methods, including conventional grid-based approximations Davis07 or more sophisticated Monte Carlo algorithms Chopin12SMC which we show increase the computational efficiency.

Interestingly, this paper not only deals with a particular problem of Network Science. As a matter of fact, we show that the mathematical formulation can indeed be applied to signals of arbitrary origin –not necessarily random walkers navigating a network–, and the multiplexity estimation framework reduces in the general case to a stochastic decomposition of the signal in terms of an effective multiplex network, whose layers play the role of independent dynamical modes. More concretely, non-Markovian dynamics can be thereby decomposed into a stochastic combination of diffusive modes by projecting the dynamics into an appropriate hidden jump-Markov model.

We validate the proposed methodology with numerical simulations of (i) random walks navigating hidden multiplex networks (thereby reconstructing the true architecture of the hidden multiplex networks) and (ii) both Markovian and non-Markovian continuous stochastic processes (thereby reconstructing an effective multiplex decomposition where each layer accounts for a different dynamical mode). Furthermore, we state and prove two existence theorems guaranteeing that such multiplex decomposition is always possible for any finite order Markovian and infinite order (fully non-Markovian) processes. Specifically, we show that random sequences generated by those processes can be exactly reconstructed as a random walk over an effective multiplex network. Finally, we apply our method to experimental recordings of two complex systems of different nature, and show that the method can be leveraged to decompose noisy, non-Markovian processes into alternating combinations of simpler dynamics and extract valuable information accordingly.

The article is structured as follows: in Sec. II we propose the methods for multiplexity detection and layer estimation, and we explore their performance in a few examples. In Sec. III we discuss the analogy between multiplexity unfolding and the decomposition of non-Markovian processes as multiple-layer Markovian processes and therefore extend the methodology to continuous time processes. We also state and discuss the implications of two existence theorems, whose proofs are put in two appendices for readability. In Sec. IV we address real-world scenarios where we analyse two sets of experimental recordings, namely human mobility in an online environment and traces of RNA polymerase, that further showcase the applicability of the method. In Sec. V we provide a discussion of our results. Mathematical details and additional examples can be found in the Appendices.

II Description of the methods and illustrative examples

Multiplex networks are the most ubiquitous class of multilayer networks. They are a natural model for online social networks lamb , where a given individual can communicate with others via different platforms (e.g. Facebook, Twitter, email, etc) or transportation networks barth ; Cardillo , where a set of locations can be connected in a multimodal way (e.g. bus, train, underground, etc).

A multiplex network is defined by a set of $L\geq 1$ interaction layers (networks), all of them having the same set of $K$ nodes but different topology (different edge set), with the peculiarity that each node has a replica in each layer (see Fig. 1a for an illustration). This structure is thereby fully described by a set of adjacency matrices $\{{\bf A}^{(\ell)}\}_{\ell=1}^{L}$, where $A^{(\ell)}_{\alpha\beta}=1$ if there is an edge between nodes $\alpha$ and $\beta$ at layer $\ell$ and zero otherwise. For simplicity we label the different layers of the multiplex network with Roman letters ($i,j$, etc), and the nodes of each layer with Greek letters ($\alpha,\beta$, etc).

We consider a random walker navigating a multiplex Manlio14 defined as follows: jumps between layers are governed by a Markov chain with $L\times L$ transition matrix ${\bf R}_{L}$ ($R_{ij}$ is the probability to jump from layer $i$ to layer $j$) while the dynamics within each individual layer $\ell$ is also Markovian and determined by a $K\times K$ transition matrix ${\bf T}^{(\ell)}$ (where $T^{(\ell)}_{\alpha\beta}$ is the probability to walk from node $\alpha$ to node $\beta$ at layer $\ell$). For simplicity we only consider diffusive dynamics where at each time step the walker at node $\alpha$ on layer $\ell$ (i) remains in the same layer with probability $R_{\ell\ell}=1-r$ or instantaneously jumps with uniform probability $R_{\ell\ell^{\prime}}=r/(L-1)$ to a different layer $\ell^{\prime}$, and subsequently (ii) diffuses to one of the neighbours of node $\alpha$ in the chosen layer according to the layer internal dynamics (given by ${\bf T}^{(\ell)}$ or ${\bf T}^{(\ell^{\prime})}$). Notice that this type of dynamical model can be mathematically formalised in terms of a jump-Markov affine system, as defined in the field of Control Theory (see Garulli and references therein for a review).

When $r\ll 1$, i.e. when walkers tend to remain in the same layer, this navigation model might mimick for instance human mobility in multilayered transportation networks gallotti1 ; gallotti2 , where multimodality is minimised to avoid waiting times related to connections between different modes.

In the particular case when layers have a simple cycle-graph topology, the jump-Markov model described above is also reminiscent of the so-called discrete flashing ratchet model ratchet ; ratchet2 , better known as a Parrondo game Parrondogame ; Dinis . This is a paradoxical gambling strategy that allows winning in loosing scenarios, where gamblers can alternate between two different strategies (game A, layer 1; game B, layer 2), each of them having different rules and winning probabilties. Our model for $L=2$ can be seen as a variant of a Parrondo gambler that plays with probabilities $r$ and $1-r$ with two different biased coins.

More generally, Brownian ratchets are paradigmatic models used in Nonequilibrium Physics to describe the transport of Brownian particles embedded in periodic, asymmetric energy potentials, a paradigm originally proposed by Smoluchowski Smoluchowski and popularised by Feynman Feynman in the context of thermodynamic engines, and further shown to be a minimal model system for molecular motors in biophysics ratchet2 . Again, a Brownian particle subject to a periodic asymmetric potential that is switched on and off stochastically is formally equivalent to our random walker navigating over a multiplex with $L=2$ cycle graphs with different transition matrices 111In the definition above time is considered discrete, however a similar methodology can be applied for continuous time dynamics by sampling the continous time series..

The general navigation process on a multiplex network discussed above can be expressed as a stochastic process fully described by an infinite two-dimensional time series $\{X(t),\ell(t)\}_{t=1}^{\infty}$ where $X(t)\in\{1,\dots,K\}$ and $\ell(t)\in\{1,\dots,L\}$. As in real-world scenarios the multiplex nature of the system is not always empirically accessible, the layer indicator $\ell(t)$ is usually hidden and the only observable is the sequence of node locations $X(1),X(2),\cdots$ . In such a situation, we only have experimental access to partial information of the process, described by a finite sequence of observations of the variable $X$: ${\cal O}=\{X(t)\}_{t=1}^{N}$. This is formally equivalent to observing a dynamical process on the aggregated (projected) network. Hence the question: is it possible to discern if the system is multiplex and in that case, to estimate the most probable number of layers if we only have access to $\mathcal{O}$? We now propose and test a novel method to achieve this highly non-trivial task.

II.1 Method to detect multiplexity

Consider the Markov switching walker discussed above, navigating on a (multiplex) network from which we only have access to coarse-grained information given by $\cal O$. Initially assuming $X(t)$ is a Markov process, we can estimate directly from ${\cal O}$ the (monoplex) transition matrix ${\bf Q}$ that would describe such a Markovian dynamics. Accordingly, we can define a Markov chain associated to ${\bf Q}$ and generate a a Markovian surrogate of the original process $\{Y(t)\}_{t=1}^{N}$. If the underlying network was truly monoplex, then $X(t)$ would be actually Markov and $X(t)$ and $Y(t)$ would then have asymptotically equivalent statistics, $P_{X}(Z_{1},\dots,Z_{m})=P_{Y}(Z_{1},\dots,Z_{m})$, for all possible sequences $Z_{1},\dots,Z_{m}$ of any arbitrary length $m$. For multiplex structures however, losing information of $\ell(t)$ in general breaks Markovianity and therefore $X(t)$ is typically non-Markovian. Accordingly, $X(t)$ and $Y(t)$ now share the same joint distributions only up to blocks of size $m=2$. For blocks of size $m\geq 3$ their statistics may differ $P_{X}(Z_{1},\dots,Z_{m})\neq P_{Y}(Z_{1},\dots,Z_{m})$. To quantify such difference, we make use of the $m-$th order Kullback-Leibler divergence rate roldan2010 ; roldan2012 between data blocks of size $m\geq 1$:

[TABLE]

where ${\cal B}(m)$ enumerates all the blocks of size $m$. The statistic $\mathcal{D}(m)$ is semi-positive definite and vanishes only when the joint probabilities coincide cover . Thus by construction $\mathcal{D}(1)=\mathcal{D}(2)=0$. The Markovianity-breaking criterion implies that if $\mathcal{D}(m)>0$ for $m\geq 3$ then the underlying dynamics is multiplex 222Even if $X(t)$ is Markovian, the matrix ${\bf Q}$ cannot be computed exactly, but estimated from the observations. Similarly, one can only estimate accurately the probabilities $P_{X}(Z_{1},\dots,Z_{m})$ for moderate values of $m$. Therefore, in practical scenarios the decision rule on whether the underlying network is multiplex or not requires to introduce an error threshold $\epsilon>0$ such that $\mathcal{D}(3)>\epsilon$ implies multiplexity. Such discrimination criterion may be formally described as a statistical test and error bounds may be obtained under mild assumptions..

As a proof of concept, we first consider the simple scenario where a random walker navigates over a two-layer multiplex ring (each layer is a cycle graph of $K$ nodes), a model compatible with a discrete flashing ratchet as commented before. In the first layer, we define a Markov chain with homogeneous transition probabilities $T^{(1)}_{\alpha+1,\alpha}=2/3;\ T^{(1)}_{\alpha,\alpha+1}=1/3$ and $T^{(1)}_{\alpha\beta}=0$ if $\beta\neq\alpha+1\mod K$ or $\beta\neq(\alpha-1)\,\text{mod}\,K$. A random walker diffusing in this layer will have an induced current in the direction of decreasing node indices. In the second layer, we define a different Markov chain with homogeneous transition probabilities $T^{(2)}_{\alpha+1,\alpha}=1/2;\ T^{(2)}_{\alpha,\alpha+1}=1/2$ and $T^{(2)}_{\alpha\beta}=0$ otherwise, i.e., an unbiased random walk. While we can always estimate ${{\bf Q}}$ numerically from the observed time series, in this simple case it is easy to derive it analytically: $Q_{\alpha\beta}=W_{1}T^{(1)}_{\alpha\beta}+W_{2}T^{(2)}_{\alpha\beta}$, where $W_{1}$ ($W_{2}$) is the probability of finding the walker in layer $\ell=1$ ($\ell=2$). Now, since in this case $R_{ij}=R_{ji}=r$, the system is symmetric with respect to the switching process and the walker spends on average the same amount of time in each of the two layers, $W_{1}=W_{2}=1/2$, and then $Q_{\alpha\beta}=T^{(1)}_{\alpha\beta}/2+T^{(2)}_{\alpha\beta}/2$. For this specific example, we thus find $Q_{\alpha+1,\alpha}=7/12,\ Q_{\alpha,\alpha+1}=5/12$.

The left panel of Fig. 2 shows numerical results of $\mathcal{D}(m)$ as a function of the block size $m$ for different switching rates $r$. In order to deal with finite-size effects (which increase exponentially with $m$), we systematically increase the size of the walker series under study as a function of $m$, taking series of size $N(m)=N_{0}\cdot 2^{m}$ data extracted from the original system and from the corresponding Markovian surrogates. We used $N_{0}=10^{5}$ although smaller values yield qualitatively equivalent results. As expected, $\mathcal{D}(1)=\mathcal{D}(2)=0$, meaning that $Y(t)$ is a faithful Markovian surrogate of $X(t)$. Furthermore $\mathcal{D}(m)>0$ for $m\geq 3$, meaning that $X(t)$ is non-Markovian roldan2012 and hence the underlying network is correctly identified as a multiplex. This result is robust for a quite large range of values of $r$ (see Appendix B, and note that $r=0.5$ is a trivial exception), meaning that the method works even if the walker makes fast switches between layers. A similar scenario is found if we tune the transition probabilities such that no net induced current is found (see appendix B), pointing out that multiplexity can be unraveled even in that case. In the right panel in Fig. 2 we plot $\mathcal{D}(3)$ for different series sizes to showcase how finite size effects vanish as the series gets larger. Notably with a few thousands data points we can already accurately detect multiplexity.

Furthermore we have demonstrated the flexibility and robustness of our method to detect multiplexity in a range of additional scenarios including (i) layers with increasingly different and disordered topologies controlled by both rewiring and edge addition, (ii) Erdos-Renyi graphs and (iii) similar scenarios on larger graphs. For all cases we find a correct multiplexity detection and good scalability (see Appendix B).

II.2 Method to quantify multiplexity

The Markovianity-breaking phenomenon which we have exploited only provides a means to discriminate whether hidden layers do exist in the multiplex, but not to quantify the number underlying layers. In order to bridge this gap, we now make use of statistical inference tools to define a model selection scheme PRLNewman . We assume that two models are different if they have a different number of layers. Accordingly, the number of layers $L$ is now modelled as a random variable with prior probability mass function $P_{0}(L)$. Given the value of $L$, the motion of the random walker is determined by the Markov-switching of layers, with transition matrix ${\bf R}_{L}$, and Markov walks within each layer characterised by ${\bf T}_{L}=\{{\bf T}^{(\ell)}\}_{l=1}^{L}$. Assuming prior probability density functions for these parameters $p_{0,R}({\bf R}_{L})$ and $p_{0,T}({\bf T}_{L})$, the likelihood of a given model with $L$ layers conditional on the observed data $\{X(t)\}_{t=1}^{N}$ reads

[TABLE]

where $\mu$ is a suitable reference measure for ${\bf T}_{L}$ and ${\bf R}_{L}$ (see appendices C and D for technical details). In general, this multidimensional integral cannot be computed exactly and needs to be approximated numerically. Then, the number of layers $L$ in the system that generated the data $\{X(t)\}_{t=1}^{N}$ can be detected using a maximum a posteriori (MAP) criterion, i.e.,

[TABLE]

where $\hat{L}^{\text{MAP}}$ is the estimate of the actual number of layers, $2\leq L_{+}<\infty$ is the (assumed) maximum possible value of $L$, and $P(L|\{X(t)\}_{t=1}^{N})\propto P(\{X(t)\}_{t=1}^{N}|L)P_{0}(L)$ is the a posteriori probability mass function of $L$ given the observations $\{X(t)\}_{t=1}^{N}$.

The practical computation of the MAP estimator in Eq. (3) can be addressed in different ways. The classical literature on hidden Markov models (HMMs) Rabiner89 ; Fessler94 ; Ghahramani01 suggests the use of the expectation-maximisation (EM) algorithm (in various forms) to compute approximate maximum likelihood (ML) estimators, $\hat{\bf T}_{L}$ and $\hat{\bf R}_{L}$, of the parameters and then assume that $P(\{X(t)\}_{t=1}^{N}|L)\approx P(\{X(t)\}_{t=1}^{N}|\hat{\bf T}_{L},\hat{\bf R}_{L})$ in order to compare the models (i.e., one tries to optimise the parameters instead of averaging over them). This approach has also been applied to jump Markov affine systems Ozkan and relies on standard techniques but it has several drawbacks and is therefore not adopted here. Actually, the equation of the model likelihood with the parameter likelihood easily breaks down when the parameter estimates are poor (e.g., because of overfitting). Another major disadvantage of EM is that it converges locally, and thus performs badly when the parameter likelihood is multimodal or when the parameter dimension varies significantly for different models. More sophisticated parametric schemes have been proposed (see, e.g., Siddiqi07 ) however they are still subject to these fundamental limitations. Integration in (2), which we adopt here, has been favoured theoretically but criticised practically because of the computational cost of approximating $P(\{X(t)\}_{t=1}^{N}|L)$ numerically Ghahramani01 . However, we have found that state-of-the-art variational Bayes McGrory09 or adaptive importance sampling Cappe08 methods can be applied effectively up to moderate values of $L$. See Appendices C3 and D for further discussion, including examples of using both deterministic integration and the adaptive Monte Carlo sampler from Koblents15 .

To illustrate the MAP model selection method given by Eq. (3), we consider again the discrete flashing ratchet model formed now by $L=3$ ring-shaped layers with $K=3$ nodes and homogeneous transition probabilities (see Fig. 1b for an illustration and Appendix D for other examples). In this example, the probability to stay in the same layer is $R_{\ell\ell}=1-r=0.84$ and the probabilities for the walker to move from node $\alpha\to\alpha+1$ are $T_{\alpha,\alpha+1}^{(1)}=0.16$, $T_{\alpha,\alpha+1}^{(2)}=0.76$ and $T_{\alpha,\alpha+1}^{(3)}=0.24$ respectively ($T^{(\ell)}_{\alpha+1,\alpha}=1-T_{\alpha,\alpha+1}^{(\ell)}$ for every $\ell$, and $T^{(\ell)}_{\alpha\beta}=0$ for any other $\alpha$ and $\beta$).

In order to evaluate the likelihood function $P(\{X(t)\}_{t=1}^{N}|L)$, we need to approximate the integral in (2) for each value $L=1,\dots,4$ as discussed previously. Here we assume $L_{+}=4$. A very simple strategy is to evaluate it via numerical integration over a deterministic grid of 19 points on the interval $(0,1)$ for each unknown parameter. For the case $L=1$ this reduces to a single unknown, $T^{(1)}_{\alpha,\alpha+1}$, and for $L>1$ there are $L+1$ unknowns: $r$ and $T^{(\ell)}_{\alpha,\alpha+1}$ for $l=1,\ldots,L$. Note that the particular choice of transition probabilities was taken to make the problem more challenging, as these values are not commensurate with the integration grid points. We use an unbiased prior probability density function $p_{0,R}({\bf R}_{L})$ given by a uniform probability density function on $(0,1)$ for the unknown parameter $r$, while the prior $p_{0,T}({\bf T}_{L})$ is used to penalise system configurations with two or more identical layers 333Note that a multiplex with $L=2$ and ${\bf T}^{(1)}={\bf T}^{(2)}={\bf T}$ is equivalent to a monoplex with transition matrix ${\bf T}$.. For this particular numerical experiment the penalising prior is $p_{0,T}({\bf T}_{L})\propto\min_{(\ell,\ell^{\prime})}|T_{\alpha,\alpha+1}^{(\ell)}-T_{\alpha,\alpha+1}^{(\ell^{\prime})}|$, i.e., the prior probability density function of a given configuration ${\bf T}_{L}$ is proportional to the minimum distance between any pair of matrices ${\bf T}^{(\ell)}$ and ${\bf T}^{(\ell^{\prime})}$. Since in this example each scalar ${\bf T}^{(\ell)}$ fully characterises layer $\ell$, this prior simply penalises configurations where two or more layers are very similar and thus avoids overfitting. In general, the prior $p_{0,T}({\bf T}_{L})$ is set to penalise models where pairs of layers have similar transition matrices, to avoid redundancy, so $p_{0,T}({\bf T}_{L})\propto\min_{(\ell,\ell^{\prime})}||{\bf T}^{(\ell)}-{\bf T}^{(\ell^{\prime})}||$, for some chosen norm $||\cdot||$.

Figure 3 shows that the true model is easily estimated with the proposed scheme as the posterior probability emphatically peaks at $L=3$, which implies that $\hat{L}^{\text{MAP}}=3$. Without penalisation –i.e., with a uniform prior $p_{0,T}$– we obtain multiple equivalent solutions involving layers with identical values of the estimated parameters.

Now, it is well known that direct, grid-based deterministic integration of the posterior probability is intuitive but computationally inefficient. Accordingly, we have further considered alternative approximations of the integral in (2) using a nonlinear population Monte Carlo (NPMC) algorithm Koblents15 , and effectively reduced the runtime by a factor close to $100$ on the same computer for the example of Fig. 3 (see Appendices C and D for details). Actually, our NPMC also includes an importance sampling procedure by which an efficient grid of the parameter space is obtained. This procedure meshes in a tight way the regions where the posterior probability density of the parameters is high, and in a sparse way in the regions where the posterior probability density is low. Accordingly, Monte Carlo sampling this mesh is guaranteed to sample parameters with high likelihood (this is a global optimization process, at odds with EM), and a simple inspection of the likelihoods of each sample allows us to robustly decide which is the one with higher likelihood. Accordingly, we are able to infer that the maximum of the posterior probability density $P(\{X(t)\}_{t=1}^{N}|{\bf T}_{L},r)p_{0,T}({\bf T}_{L})$ (i.e., the integrand of Eq. (2) with uniform prior for $r$) is indeed attained at the true value of ${\bf T}_{3}$ (see appendix D5 for details). We therefore conclude that our model selection scheme correctly estimates not only the number of layers, but also the transition probabilities within each layer, i.e. the full architecture.

Finally, we have also verified that the posterior probability density is smooth close to its maximum as perturbations $\tilde{\bf T}_{3}={\bf T}_{3}+\delta{\bf T}$ systematically yield a smaller posterior probability density function for sufficiently large $N$ (see figure 21).

III A multiplex decomposition of non-Markovian dynamics

As it happens for any Bayesian inference-based method CNM ; GSP ; PRLNewman , our approach in principle would provide conclusive indication of the hidden multiplexity in the case the prior –intralayer dynamics are diffusive– is a reasonable assumption. Notwithstanding, the methodology proposed here actually extends above and beyond the reconstruction of hidden multiplex architectures using walkers with partial information. As a matter of fact, a similar approach can be considered even when the architecture is truly single-layered. Suppose for instance that a given observed time series was truly the outcome of a non-Markovian dynamics running on a physical single-layered network. In that case, our multiplex reconstruction method would still provide the most probable multiplex model, with $L>1$ due to lack of walker’s Markovianity. The key difference is that now layers in the hidden multiplex would be providing the most probable effective multiplex reconstruction with Markovian intralayer dynamics that would yield such complex dynamics. Incidentally, note that this brand of effective models is used in community detection in single-layered networks, which result from finding the optimal number of effective groups of nodes which maximise a certain likelihood function PRLNewman .

In what follows, we first capitalise on this new interpretation to extend our previous analysis on random walks on graphs to continuous stochastic processes, and then we present two existence theorems which guarantee that this stochastic decomposition is universally applicable to random sequences with arbitrary memory.

III.1 Extension to continuous processes

When the original dynamics is continuous we can discretise motion and embed the original time series into a simple graph topology, for example via the transformation $\{X(t)\}\to\{\tilde{X}(t)\}$ given by $\tilde{X}(0)=0$ and

[TABLE]

and apply our multiplexity detection methods to the discretised trace $\{\tilde{X}(t)\}$. To illustrate and validate this extension, we have considered two continuous stochastic processes: (i) an Ornstein-Uhlenbeck process (Markovian) governed by the stochastic differential equation

[TABLE]

where $\xi$ is a Gaussian white noise with zero mean $\langle\xi(t)\rangle=0$, amplitude $\sigma$ and autocorrelation $\langle\xi(t)\xi(t^{\prime})\rangle=\sigma^{2}\delta(t-t^{\prime})$ and (ii) a generalisation of the preceding Langevin equation where the noise term is not white anymore but has some colour. Such process can be described by the following Langevin equation

[TABLE]

with $\eta(t)$ being itself defined by an Ornstein-Uhlenbeck process such that $\langle\eta(t)\rangle=0$ and $\langle\eta(t)\eta(t^{\prime})\rangle\propto\exp(-|t-t^{\prime}|/\tau)/\tau$ with $\tau>0$ the correlation time of the noise. Note that in Eqs. (5) and (6) the dot denotes time derivative. When the noise term $\eta$ is not white anymore, the Langevin equation (6) generates non-Markovian trajectories for the variable $x$. Interestingly, Eq. (6) is attracting considerable attention in soft matter as a minimal model for the non-Markovian dynamics of the position of a passive particle immersed in an active (e.g. bacterial) bath active1 ; active2 .

We perform numerical simulations of these processes using an Euler-Mayurama algorithm and subsequently embed the resulting trajectories of $10^{4}$ data in a cycle-graph topology (see Fig. 1b) via Eq. (4) with $K=4$, and then we have applied our multiplex detection and estimation protocol. In the top panels of figure 4 we plot ${\cal D}(m)$ and the outcome of a layer estimation. For the Ornstein-Uhlenbeck process the method does not detect multiplexity (and a layer estimation confirms that a model with $L=1$ layer is the most probable one), in good agreement with the fact that the Ornstein-Uhlenbeck process is indeed a Markovian process. On the other hand, the Langevin equation with coloured noise has a nontrivial memory kernel and generates non-Markovian dynamics, correctly captured by the fact that ${\cal D}(3)>0$. In this case the dynamics optimally decomposes into a Markov switching combination of $L=2$ layers.

To further validate these results, we now generate synthetic trajectories from the estimated multiplex models with $L$ layers. Importantly, we should recall that our layer estimation method makes use of an importance sampling algorithm to concentrate the Monte Carlo search in the regions of the parameter space with large likelihood. As a byproduct, we are are able to estimate –given the number of layers $L$– the parameters of the most likely model, and therefore we can now generate synthetic series from this fitted model. We then compare the statistics of synthetic series from the most likely model with $L$ layers with those of the original Ornstein Uhlenbeck and non-Markovian processes, respectively. In the bottom panels of Fig. 4 we plot the $m$-th order Kullback-Leibler divergence between the original (discretised) series and the synthetic series generated by the reconstructed multiplex model. We confirm that the series generated by the model with $L=1$ and $L=2$ respectively are the ones that show higher similarity with the Ornstein-Uhlenbeck (Markovian) and non-Markovian series, respectively. Finding $L=2$ layers in the non-Markovian case is indeed reasonable having in mind that Eq. (6) allows a well-known decomposition as the following 2D (Markovian) stochastic differential equation

[TABLE]

with $\xi$ a Gaussian white noise process and $\beta>0$ a positive constant.

That is to say, in even more general terms, our methodology provides a mathematically sound solution for the stochastic projection of non-Markovian dynamics onto a base of simple diffusive dynamical modes. For random sequences which were originally generated by a random walker navigating a hidden multiplex network, the method easily reconstructs such hidden architecture, whereas for general random sequences (such as the ones generated by Eqs. 5 and 6) the effective multiplex model provides a good reconstruction of the original dynamics. Hence the question: is this type of hidden jump-Markov model able to fully (i.e., exactly) reconstruct any random sequence (i.e. high order Markovian and fully-non Markovian)? These questions are responded affirmatively in the next subsection, where we state and prove two theorems which address these matters.

III.2 Exact decomposition of Markovian and non-Markovian dynamics via multiplex models

In this subsection we pose the question of whether the statistics of arbitrary random sequences $\{X(t)\}_{t\geq 0}$, possibly with long memory, can be recovered exactly using a hidden jump-Markov model of the class described in Section II (from now on we refer to it as a multiplex model). The answer is positive (in a probabilistic sense to be made precise), as summarised by two representation theorems.

In particular, we first state and prove that every Markov model of finite (but arbitrarily large) order $h$ can be recast into a multiplex model. Then we address the representation of models with infinite memory, by letting $h\rightarrow\infty$, and show that, under mild regularity assumptions, they admit a compact representation in the form of a model with an uncountable (continuous) set of hidden layers. For readability, we have relegated as many technical details as possible (including all the mathematical proofs) to Appendices F and G, and only discuss here the key results and implications.

III.2.1 Representation of Markov models of order $h$

Let us consider a discrete state space with $K$ elements, $\mathbb{K}=\{1,\ldots,K\}$ (this will be the node set of the multiplex). A Markov sequence of (integer) order $h\geq 1$ is defined by the set of $K^{h+1}$ probability masses

[TABLE]

where $X(t-h:t-1)=\{X(t-h),\ldots,X(t-1)\}$ and $i_{t-h:t}=\{i_{t-h},\ldots,i_{t}\}\in\mathbb{K}^{h+1}$ is a sequence of $h+1$ state space observations. If we fix $i_{t-h:t-1}$, then $P_{\mathbb{K}}^{h}(i_{t}|i_{t-h:t-1})$ is a probability mass function (pmf) and, as a consequence, $P_{\mathbb{K}}^{h}(i_{t}|i_{t-h:t-1})\in[0,1]$ and $\sum_{i_{t}\in\mathbb{K}}P_{\mathbb{K}}^{h}(i_{t}|{i_{t-h:t-1}})=1$.

Now, any Markov model of finite order $h$ can be transformed into an equivalent multiplex model with a sufficiently large, but finite, number of layers $L$, as formally stated below.

Theorem 1

For every Markov model of order $h<\infty$ on the state space $\mathbb{K}=\{1,\ldots,K\}$ there exists an equivalent multiplex model, with observation space $\mathbb{K}$, and $L=K^{h}$ layers.

The proof of Theorem 1 is technical and is therefore put in Appendix F. Theorem 1 guarantees that every random Markov sequence with finite memory $h$ can be represented by a multiplex model with $K^{h}$ layers (i.e., this is an existence theorem). The theorem does not state that this representation is unique, though, and it does not state that it is minimal either. According to the numerical evidence given in the previous sections, we conjecture that a suitable selection of ${\bf R}_{L}$ and ${\bf T}_{L}$ (e.g., using the estimation methods described in this paper) can yield an accurate representation of a sequence of order $h$ with considerably less than $K^{h}$ layers.

Let us also remark that there is no contradiction between Theorem 1 and our earlier claim that multiplex models can represent systems with infinite memory. Indeed, depending on the choice of its parameters ($L$, ${\bf R}_{L}$ and ${\bf T}_{L}$), a multiplex system can yield random sequences with either finite or infinite memory. For example, in the proof of Theorem 1 we explicitly construct a multiplex model that matches the transition probabilities of a Markov model of order $h$. On the other hand, multiplex models where the transitions between layers are independent of the walker’s past trajectory yield sequences with infinite memory (except for pathological cases).

III.2.2 Representation of sequences with infinite memory

In a second part, we extend the previous existence theorem to a wide class of infinite memory (i.e. fully non-Markovian) models. Let $\{X_{h}(t)\}_{t\geq 0}$ denote a random sequence on $\mathbb{K}=\{1,\ldots,K\}$ that can be represented exactly by a Markov model of order $h$. We consider here the class of sequences with infinite memory that can be obtained from Markov models as $h\rightarrow\infty$ and refer to them as “Markov-$\infty$” sequences. To be precise, we say that $\{X(t)\}_{t\geq 0}$ is the limit of $\{X_{h}(t)\}_{t\geq 0}$ as $h\rightarrow\infty$, and write

[TABLE]

when we can approximate the transition probabilities of the sequence $\{X(t)\}_{t\geq 0}$ with an arbitrarily small error using a Markov model of sufficiently large order. Specifically, we need to introduce the following technical definition:

Definition 1

The random sequence $\{X(t)\}_{t\geq 0}$ is Markov-$\infty$ if it satisfies the regularity conditions below:

• (C1)

The joint probability of any sequence of states vanishes uniformly with the length of the sequence, i.e.,

[TABLE]

• (C2)

There exists a sequence of Markov models $\{X_{h}(t)\}_{t\geq 0}$, $h=1,2,\ldots$, such that for any $\epsilon>0$, arbitrarily small, there exists $h^{\prime}\in\mathbb{N}$, sufficiently large, that guarantees

[TABLE]

for every $t>h$ and

[TABLE]

for every $t\leq h$, whenever $h>h^{\prime}$.

Let $\{X_{L}(t)\}$ represent a random sequence generated by a multiplex model with $L$ layers. Since every Markov-$\infty$ model is the limit of a sequence of Markov systems with increasing order, $h\rightarrow\infty$, they can also be obtained (via Theorem 1), as the limit of a sequence of multiplex systems as the number of layers grows, $L\rightarrow\infty$. Moreover, it turns out that the limit $\lim_{L\rightarrow\infty}\{X_{L}(t)\}_{t\geq 0}$ can be interpreted itself as a multiplex model with an uncountable set of layers and a first-order Markov system on $\mathbb{K}$ associated to each layer. This is made precise by the following theorem:

Theorem 2

Let $\{X(t)\}_{t\geq 0}$ be a Markov-$\infty$ random sequence on $\mathbb{K}$. There exists a Markov kernel

[TABLE]

where ${\mathcal{B}}([0,1))$ denotes the Borel $\sigma$-algebra of subsets of $[0,1)$, a probability measure ${\mathcal{M}}_{0}:{\mathcal{B}}([0,1))\rightarrow[0,1]$, and an uncountable family of $K\times K$ transition matrices

[TABLE]

such that

[TABLE]

for every $i_{0:t}\in\mathbb{K}^{t+1}$, where $P_{i_{0}}=P(X(0)=i_{0})$.

See Appendix G for a proof. Theorem 2 indicates that multiplex network models can be generalised to obtain probabilistic systems with an infinite and uncountable number of layers, which are flexible enough to represent (i.e., exactly recover) a broad class of random sequences with infinite memory.

IV Applications to experimental data

We now apply our methodology to two experimental recordings of completely different nature: the mobility dynamics of human players of an online videogame and the dynamics of polymerases during RNA transcription.

IV.1 Application to mobility: the Pardus universe

Human and animal mobility mobilityJesus1 ; mobilityJesus2 ; mobilityJesus3 are often described by dynamical processes on single-layer networks, and interestingly, has been found to signatures of memory mobility ; mobility2 ; Salnikov , which can be interpreted as a deviation from Markovianity. Can such lack of Markovianity be interpreted as being the result of a Markovian dynamics taking place on a hidden multiplex network? In the case where mobility takes place across (hidden) multimodal transportation systems (as when we collect GPS traces of urban mobility), layers could be physical (underground, bus, car, etc). On the other hand, animal foraging dynamics are clearly different and alternate during day and night day . In a similar vein, human mobility patterns change when switching from work/leisure styles leisure : these would be cases where layers would be effective, rather than physical. There are a large variety of problems involving the aforementioned scenarios which could be amenable to our approach. To guarantee computational efficiency, we would only require the network over which the agents move not to be too large, something that in the general case can be achieved by coarse-graining the network via community detection. Mathematically, within this context a non-Markovian process running on a monoplex vs a Markov switching process running on a multiplex are equally valid models, much in the same way a function is equivalent to its Fourier series representation. Still, we consider this new interpretation not just suggestive but also parsimonious from a cognitive point of view, and therefore might be of relevance in the study of memory in search processes.

To illustrate this type of application, we consider experimental mobility trajectories performed by players (agents) in a virtual environment: the Pardus universe (see Szell for details). The Pardus universe is a multiplayer online role-playing videogame which is used as a large scale socio-economic laboratory to study mobility in a controlled way. It consists of an (online) physical network with $K=400$ nodes, a networked universe with social and economic activities, where the players the game move around. The mobility traces of these players can then be naturally symbolised by coarse-graining the original network via community detection into a network of 20 non-overlapping communities wired by a $20\times 20$ weighted adjacency matrix with complex topology Szell .

Szell et al analysed online player’s trajectories and reported evidences of long-term memory in the diffusing patterns of agents mobility in this universe Szell . Here we apply the multiplexity detection statistic to a long time series obtained by concatenating individual traces. In Fig. 5 we show the numerical results (black dots), along with the null case of a Markovian diffusion (red squares) over the same network. We find $\mathcal{D}(3)>0$ as a footprint of multiplexity, which can be here interpreted as the hidden presence of several effective layers. Szell et. al. introduced a long-term memory model to account for the mobility dynamics, and the heterogeneity of players. An alternative interpretation is that players are performing simple diffusion dynamics, but are switching stochastically between different effective dynamical regimes. The challenge, which we leave as an open question for future work, would be to assign a social or perhaps cognitive meaning to the different layers. In this sense, our method is here closer to an unsupervised clustering paradigm than to a supervised one. Since our experimental trace is given by the concatenation of traces of different walkers, the different effective layers could also be revealing here a taxonomy of different game strategies.

IV.2 Application to biology: transcription by eukaryotic RNA polymerases

It has been shown that the ratchet mechanism plays an important role in active transport by molecular motors in living cells ratchet2 . A paradigmatic example of a molecular motor that can be well described by a ratchet mechanism is RNA polymerase (RNAP). RNAPs are macromolecular enzymes responsible for the transcription of genetic information encoded in the DNA into RNA  CramerReview . During transcription the spatial location of RNAPs along the DNA template exhibits noise due to thermal fluctuations. In addition the dynamics of polymerases exhibits switching between two different dynamical regimes: an active polymerisation state called elongation and a passive or diffusive state called backtracking galburt . While elongating, RNAP moves along the DNA template with a net velocity of the order of $\sim 20$ nucleotides per second with $0.34$nm being the distance between two nucleotides. In the backtracking regime, the polymerisation reaction is stalled and RNAPs perform passive Brownian diffusion due to several types of noise (e.g. thermal and chemical). Optical tweezers enable measuring the motion of a single RNAP during transcription of a single DNA template at a basepair resolution abbon ; galburt . Single-molecule experimental RNAP traces are however extremely noisy and the problem of identifying the possible mixture of different underlying dynamical processes, including transitions from elongation to backtracking as well as the hidden presence of other regimes, from a single experimental trace is a challenging task.

At the single-molecule level, the dynamics RNAP corresponds to a molecular motor producing a non-Markovian dynamics which can thus be modelled as a Markov switching random walk dynamics in a biochemical multiplex with $L=2$ layers (corresponding to elongation and backtracking respectively) dang . This suggests that our methodology (using cycle graphs with homogeneous transition rates as the prior topology) is applicable and under the aforementioned assumptions we should find both $\mathcal{D}(3)>0$ and the model with $L=2$ should have a clear maximum likelihood. To test such prediction, we have applied our complete methodology to five single-molecule experimental traces of the position $X$ of RNA polymerase I (Pol I) from yeast S. cerevisiae, obtained with a dual-trap optical tweezer setup in the assisting force mode lisica . We systematically choose the first $10^{4}$ data points at $1$kHz sampling rate (see Figure 6 and Appendix E) to keep the time series short and make the inference problem harder (for instance, for the series shown in Figure 6 it is not easy to visually distinguish the elongation and backtracking regimes). Since the original traces have continuous state space, we discretise these seried by embedding the experimental recordings into a cycle graph via the transformation proposed in equation 4 for $K=4$. An example is shown in the top-right panel in Fig. 6. In the left-bottom panel of the same figure we plot $\mathcal{D}(m)$ applied to $\{\tilde{X}(t)\}$ for each of the five experimental series, consistently showing $\mathcal{D}(3)>0$. For comparison and control of finite size effects, a similar measure is computed on a null model: a time series of $10^{4}$ data points generated by a (monoplex) Markov chain whose transition matrix has been estimated from $\tilde{X}(t)$ (red dashed line). We can conclude that the underlying dynamics requires at least two alternating dynamics –a stochastic alternation between two diffusive layers–, hence the projection onto an effectively multiplex model.

Subsequently, in the bottom right panel of the same figure we provide the results on the layer estimation using our nonlinear population Monte Carlo algorithm (convergence after 12 iterations, $10^{2}$ samples per iteration). The algorithm directly provides $\log(P({\cal O}|L))$, so assuming a uniform prior on the number of layers $\log(P({\cal O}|L))\propto\log(P(L|{\cal O}))$, i.e. the logarithm of the a posteriori probability of model $L$. The true probability of the model (in natural units) is subsequently extracted and plotted accordingly. We find that the probability is essentially one for model with $L=2$ layers and negligible for the rest. Our algorithm thus reveals that the optimal hidden multiplex has $L=2$ effective layers. One possibility is that the identification of $L>1$ is the consequence of the presence of coloured noise in the single-molecule optical tweezer transcription experiment. Quite intriguingly, additional evidence points to the presence of non-Markovianity also in the backtracking regime (see Appendix E), what might suggest the presence of coloured noise even in the backtracking mode, as in the example in Sec. III. Another possibility is, as previously discussed, that the two effective layers correspond to the biochemical mechanisms of elongation (active state) and backtracking (passive state) dynamical modes. Finally, we expect our approach to be also applicable to more complicated scenarios such as to identify the number of different nucleotides in copolymerisation processes of templates with strong disorder gaspard .

V Discussion and Conclusion

In this work we have introduced a method that both detects and quantifies the degree of multiplexity in the hidden underlying structure of a networked system by only having access to local and partial statistics of a random walker. Our working hypothesis (prior) is that there is a hidden multiplex where walkers diffuse, switching layers stochastically and diffusing over each layer. Under these circumstances, any random walker for which we only see a projection of such trajectory in the aggregated network is necessarily non-Markovian, if the number of layers is larger than one. Hence our algorithm for multiplexity detection exploits such breaking of Markovianity as a means to detect multiplexity.

Incidentally, here we have focused in the specific case of multiplex networks, where every layer has the same number of (replica) nodes. Actually, in a multiplex one can even have the same topology in each layer, where only transition weights differ, as in the case of the multiplex cycle graphs considered above. On the other hand, in a generic multi-layer network each layer will have in general different number of nodes and different topology. This latter situation can be reinterpreted as having a multiplex where in each layer we can have isolated nodes which are never reached by a walker, or forbidden transitions. Accordingly, we envisage that layers in a multiplex would be in general harder to distinguish via our method than layers in a generic multi-layer network, and as such we expect that the detection method to be easily generalisable to the multi-layer case, something that should be studied in the future.

In a second step, we have introduced a probabilistic scheme to estimate the most probable number of layers composing the hidden multiplex. Note that probabilistic model selection is not new in network inference, for instance in roger a similar concept is used to estimate the most probable combination of basic block models that accounts for a certain network topology (see also tiago1 ; tiago2 ; tiago4 ; tiago5 ), whereas in PRLNewman a probabilistic framework was developed to estimate the most probable number of communities in a single-layer network. Formally similar strategies to estimate model parameters based on $\epsilon$-machines or jump-Markov system identification have also been put forward in the Nonlinear Dynamics Crutchfield and Control Theory Ozkan communities, respectively. In our case, the posterior probabilities quantify the likelihood of having a hidden multiplex with $L$ layers, i.e. our approach for multiplex model estimation is purely Bayesian. We were able to demonstrate the validity and accuracy of this second part in simple synthetic networks (multiplex cycles) up to $L=5$ layers, as reported in the main text and appendices. Since the model selection protocol can be seen as a multidimensional Bayesian inference problem, these schemes –similarly to Hidden Markov Models and other methods in statistical inference– suffer from poor scalability: essentially the computation of the model posterior probability explodes with the number of unknowns. This is a limitation of the method, and its optimisation is therefore an open problem for future work. As a matter of fact, a simple and intuitive (although inefficient) way to estimate these posteriors is to use a (deterministic) grid integration scheme. In an effort to improve scalability and optimise such calculations we have proposed a nonlinear population Monte Carlo algorithm (described in full detail in appendices C and D) which reduces the computer runtime by a factor close to $10^{2}$, without performance loss, for all the examples where we had previously used deterministic integration. Interestingly enough, our model selection scheme not only selects the most probable number of hidden layers: by capitalising on an importance sampling routine that focuses in a small region of the parameter space where the likelihood is concentrated we can also provide a Bayesian estimation of the model parameters (i.e. the topology and transition rates), thereby estimating not only the most probable number of layers but, given that number of layers, the full architecture.

We have thoroughly illustrated the validity and scalability of the whole method by addressing several synthetic systems of varying complexity, as depicted in Sec. II and the Appendices, where we show that we can reconstruct the full architecture of the hidden multiplex network by analysing the statistics of random walks over the projected network.

Interestingly, our method can be applied to signals of arbitrary origin, extending the formalism to deal also with continuous processes –i.e. time series which are not necessarily random walkers navigating a network–, after a simple graph embedding (series discretisation). Under this extension, our hidden jump-Markov model provides a decomposition of a given non-Markovian dynamics into a Markov switching combination of diffusive modes, and thus enjoys larger generality: the reconstructed multiplex model embeds the originally non-Markovian signal into a random walk navigating an effective multiplex network, where each layer accounts for a different type of diffusive dynamics. We validated this extension by analysing canonical continuous processes (Ornstein-Uhlenbeck process and a Langevin equation with coloured noise). Furthermore, we have proved two existence theorems which guarantee that an exact reconstruction is always possible, for any type of (arbitrary) finite order Markovian and fully non-Markovian (i.e. infinite memory) process.

Finally, we applied this methodology in two experimental scenarios (a case of human mobility in an online universe and the analysis of the dynamics of RNA Polymerase) and hence showcased its applicability in real, experimental data.

To conclude, starting from the question on whether it is possible to disentangle the hidden multiplex architecture of a complex system if one only has experimental access to a projection of this architecture, in this work we have elaborated a mathematical and computational framework that actually deals more generally with the decomposition of non-Markovian dynamics. When these series are indeed traces of walkers navigating a network, under the premise of having intralayer diffusion our approach provides a workable solution for the unfolding of a multiplex network from its aggregated projection. We should make clear that, generally speaking, it is not possible to assert that the effective multiplex representation is the true architecture, much like one cannot typically claim that there exist true communities in an observed network –but rather, one says that observations are better reproduced by a multiplex model than by a monoplex one, in a similar spirit as models of networks with structure sometimes reproduce better the observations that models that lack such community structure.

In the general case our method provides a potentially useful approach to disentangle combination of dynamical regimes that appear intertwined in noisy dynamics, this being for instance the case of RNA polymerase moving in a noisy environment and stochastically switching between an active and a passive state. This suggests that applications of this work not only include Network Science but extend to other fields in Biophysics, Condensed Matter Physics, gambling or Mathematical Finance, where non-Markovian signals pervade.

**Acknowledgments. **We sincerely thank Michael Szell, Roberta Sinatra and Vito Latora for sharing data on Pardus universe and for fruitful discussions, and anonymous referees for useful comments. L.L. acknowledges funding from EPSRC Early Career Fellowship EP/P01660X/1. I.P.M. acknowledges the financial support of the Spanish Ministry of Economy and Competitiveness (projects TEC2015-69868-C2-1-R ADVENTURE and TEC2017-86921-C2-1-R CAIMAN). J.M. acknowledges the financial support the Spanish Ministry of Economy and Competitiveness (project TEC2015-69868-C2-1-R ADVENTURE) and the Office of Naval Research (ONR) Global (Grant Award no. N62909-15-1-2011).

LL, IPM and JM contributed equally to this work.

Appendix A A FEW NUMERICAL CONSIDERATIONS

**Bounds on switching rate. **Numerically, the problem of detecting multiplexity via statistical differences between the Markovian surrogate $Y(t)$ and $X(t)$ should, in principle, be easier when the switching rate $r$ is small enough, such that we allow trends associated to each layer to build up in the series, but large enough such that $D(1)=D(2)=0$ (that is, large enough such that the one and two step joint distributions are still equivalent). A simple theoretical lower bound is $r<1/m$ (as the average size of a trend is $1/r$, and this should be at least as large as the block size). Note that this bound is not tight in what refers to real-world scenarios (as for human mobility), where switching rate is normally low in order to avoid unnecessary delays, or has characteristic timescales which are much lower than the diffusion timescales (as for changing foraging mode between day and night, in animal mobility).

**Finite size bias in KLD. **As a technical remark, note that $\text{KLD}(p||q)$ diverges if the distributions $p$ and $q$ have different supports (i.e., if $q(m)=0,p(m)\neq 0$ or $p(m)=0,\ q(m)\neq 0$ for some value $m$). In order to take appropriately weight this possibility while maintaining the measure finite in pathological cases, a common procedure roldan2012 is to introduce a small bias that allows for the possibility of having a small uncertainty for every contribution. Here we introduce a bias of order ${\cal O}(1/n^{2})$ where $n$ is the series size (i.e., we replace all vanishing frequencies with $1/n$, and we normalise the frequency histogram appropriately).

Appendix B INFERRING MULTIPLEXITY IN SOME ADDITIONAL SCENARIOS

We start by providing in figure 7 an additional analysis similar to Figure 2 in the main manuscript, but where there is no induced current. The method correctly detects multiplexity in this arguably more complicated scenario. From now on we use ${\cal D}(m)$ and $\text{KLD}_{m}/m$ indistinctively to express the Kullback-Leibler divergence of order $m$ between a series $X$ and its Markovianised surrogate.

In the next sections we extend the initial study on inferring multiplexity to the case where layers have different complex topologies, departing from the situation where each layer is a cycle graph (ring). Note that when each layer has a different topology we expect the algorithm to detect more easily the underlying multiplex character (in this sense, extension of this formalism to the more general case of a multi-layer network is very promising). In particular, we explore the following additional scenarios:

• •

Scenario 2: two complete graphs with different transition matrices.

• •

Scenario 3: two layers with cycle graphs where in one of the layers we introduce a shortcut which is crossed with a probability $\epsilon$. This scenario allows us to explore small perturbations in the transition matrices with respect to the original scenario studied in the main text.

• •

Scenario 4: each layer has a different topology and dynamics, the first being a cycle graph (ring) with positive net current and the second layer is a complete graph with null net current.

• •

Scenario 5: initially having two identical layers (two cycle graphs), we introduce a number of additional edges (shortcuts) in the second one, to explore topological perturbations on our initial scenario. We also explore the scalability of this general scenario by considering the effect of increasing the number of nodes.

• •

Scenario 6: initially we consider two identical layers formed by Erdos-Renyi graphs, and then rewire in one of the layers a percentage of the edges. Transition matrices are obtained here by unbiasing the walker $T_{ij}=A_{ij}/k_{i}$. We also explore the scalability of this scenario by considering the effect of increasing the size of the graphs.

B.1 Scenario 2: Complete graphs

In this scenario we consider two layers with identical topology but different transition matrices. Here we build two replicas of a complete graph. In the first layer we define an unbiased random walker with transition matrix

[TABLE]

whereas in the second layer we introduce a parametric deviation

[TABLE]

The larger $\epsilon$, the more different the statistics of a Markov Chain over each layer separately. For $\epsilon=0$ both layers are identical and a walker diffusing over the multiplex (switching layers at a constant rate $r$) reduces to a Markov Chain over one layer, whereas for $\epsilon>0$ the process is non-Markovian if we only have access to the state $X(t)$. In figure 8 we plot the results for $D(m)$ which show that multiplexity can always be detected, and such detection is easier as $\epsilon$ increases.

B.2 Scenario 3: Controlled perturbation on one layer

In this scenario we explore the effect of a controlled perturbation in the topology of one of the layers. We start by defining $L=2$ identical layers (two rings with the same transition matrix, with an homogeneous probability to flow from $i\to i+1$ $T_{i,i+1}=1/3$). In the second replica we introduce a shortcut between two nodes, weighting the probability of traversing this node as $2\epsilon$, and biasing accordingly the rest of the edges. Accordingly the transition matrices of both layers read

[TABLE]

[TABLE]

When $\epsilon=0$ this extra edge has no effect and we should expect that the multiplex is effectively monoplex, whereas for $\epsilon>0$ both layers show gradually different structure and as such $X(t)$ is non-Markovian. Similarly as in the previous case, our methodology predicts that the network is multiplex when $D(m\leq 2)=0$ and $D(m>2)>0$. Intuitively, as $\epsilon$ increases the effect of the shortcut should be higher and thus detecting the multiplex nature of the network should be easier. In figure 9 we plot $D(m)$ for a Markov Chain walking over this topology with a constant switching rate $r=0.1$ and different values of $\epsilon$. Results support the accuracy of the method, even for small values of $\epsilon$ the multiplex nature of the network is detected.

B.3 Scenario 4: Ring versus complete graph

In this scenario we focus on a multiplex with $L=2$ where each layer is totally different. The first layer is a ring with a net current described by the transition matrix

[TABLE]

whereas the second layer is a complete graph with no net current and detailed balance everywhere:

[TABLE]

Again, we find that the method can clearly detect the multiplex nature of the network as $D(m>2)>0$.

B.4 Scenario 5: Sequentially introducing shortcuts

In this scenario we consider two identical replicas (a multiplex with $L=2$ layers) where in the second layer we add a certain number of additional edges. Originally, both replicas are rings with detailed balance (unbiased walker)

[TABLE]

In this initial case we expect $D(m)=0,\forall m$ (the multiplex is effectively monoplex). We then add to this benchmark a different number of edges (and accordingly we expect $D(m>2)>0$, and larger values when the number of added edges increases), effectively interpolating between scenario 2 and scenario 3. The adjacency matrices for three concrete cases (with no added edges, one added edge and two added edges - this latter case being equivalent to a complete graph) are depicted in figure 11, and in figure 12 we show the values of $D(m)$.

B.5 Scenario 5b: Introducing edges on larger graphs

Here we investigate the scalability of the scenario 5 and in particular we explore (i) the effect of increasing the number of nodes $K$ in each layer and (ii) the effect of increasing the number of rewired edges on the detectability. We start by defining two identical replicas of a ring with $K$ nodes with unbiased transition matrix $T_{ij}=1/2$ for $i\neq j$. In the second layer we introduce a number of shortcuts and we analyse two particular behaviors as it follows.

**Effect of node increase. ** We fix the number of shortcuts $p=2$ and vary the number of nodes $K$, and explore the dependence of $D(m)$ on $K$. As we keep the series size ${\cal N}(m)=10^{5}\cdot 2^{m}$ being independent from $K$, we expect that as $K$ increases the statistics are poorer as we need larger series to capture an equivalent number of transitions.

**Detectability as a function of the number of shortcuts introduced. ** Here we fix $K=10$ and explore the multiplex detectability as the number of shortcuts $p$ is increased in the second layer. Multiplexity is detected when $D(3)>0$, and the larger $D(3)$ the easier is such detection. In figure 14 we plot $D(3)$ as a function of the number of shortcut edges $p$.

B.6 Scenario 6: Rewiring edges

In this scenario we initially consider two replicas of the same Erdos-Renyi graph (where nodes $i$ and $j$ are connected with probability $p=0.8$, above the percolation threshold to have a connected graph). We consider three different situations, namely: (i) both layers are maintained identical, (ii) we rewire at random one edge, (iii) we rewire at random two edges. The adjacency matrices of each layer for these three cases are represented in figure 15, and we choose unbiased random walkers with layer transition matrices $P_{ij}=A_{ij}/k_{i}$. In figure 16 we plot the values of $D(m)$ for these three cases. As expected, when the layers are identical we find $D(m)=0\ \forall m$, whereas when we rewire edges from a layer the network converts into a multiplex one and $D(m>2)>0$. Also, $D(m>2)$ take larger values -and thus multiplex detection is easier- when the layers are increasingly different.

B.7 Scenario 6b: Rewiring edges on larger graphs

Finally, we consider Erdos-Renyi graphs (linking probability 0.65) with $K=10$ nodes per layer and explore the multiplex detectability as we rewire a percentage of nodes $p$ in the second layer. Multiplexity is detected when $D(3)>0$, and the larger $D(3)$ the easier is such detection. In figure 17 we plot $D(3)$ as a function of $p$. Dots are the result of an ensemble average over ER graphs realisations. For an ensemble we keep fixed the number of rewired edges and compute the effective average percentage of rewired edges (which fluctuates as each realisation of an ER graph will have a different total number of edges).

Appendix C MATHEMATICAL AND ALGORITHMIC FRAMEWORK FOR LAYER ESTIMATION

Here we formalise the problem and provide a detailed derivation of the model, the probabilistic model selection scheme and some possible algorithmic implementations of this scheme. Let us remark that the probabilistic framework described herein includes the case in which the path followed by the walker across the multiplex network cannot be observed exactly (i.e., there are observation errors) even if this is not addressed in the main text.

We use an argument-wise notation to denote probability mass functions (pmf’s) and probability density functions (pdf’s). If $X$ and $Y$ are discrete random variables (r.v.’s) then $P(X)$ and $P(Y)$ are the pmf of $X$ and the pmf of $Y$, respectively, possibly different. Similarly, $P(X,Y)$ and $P(X|Y)$ denote the joint pmf of the two r.v.’s and the conditional pmf of $X$ given $Y$, respectively. We use lower-case $p$ for pdf’s. If $X$ and $Y$ are continuous r.v.’s, then $p(X)$ and $p(Y)$ are the corresponding densities, possibly different, and $p(X,Y)$ and $p(X|Y)$ denote the joint and conditional pdf’s. We may have a pdf of a continuous r.v. $X$ conditional on a discrete r.v. $Y$, $p(X|Y)$, as well as a pmf of $Y$ given $X$, $P(Y|X)$. Most r.v.’s are indicated with upper-case letters, e.g., $X$. If we need to denote a specific realisation of the r.v., then we use the same letter but lower-case, e.g., $X=x$ or $Y=y$. Matrices and vectors are indicated with a bold-face font, e.g., ${\bf T}$.

C.1 The model

We assume that a walker travels through a multiplex network taking random moves between neighbouring nodes and, occasionally, between layers. Let $L$ denote the number of layers in the multiplex network and let $K$ be the number of nodes per layer. At discrete time $t$, the random variable (r.v.) $X(t)$ denotes the in-layer walker position. Therefore, $X(t)\in\{0,\ldots,K-1\}$ and $X(t)=k$ means that the particle is located at node number $k$ at time $t$, irrespective of the layer. The r.v. $\ell(t)$ indicates the layer at time $t$, i.e., $\ell(t)\in\{1,2,...,L\}$ and $\ell(t)=l$ means that the walker is found in layer $l$ at time $t$. The state of the walker, therefore, is given by the $2\times 1$ vector ${\bf Z}(t)=[X(t),\ell(t)]^{\top}$.

At each time step, the walker may jump across layers. This motion is assumed to be Markov and hence it can be characterised by an $L\times L$ stochastic transition matrix ${\bf R}_{L}$, where the entry $R_{ij}$, $(i,j)\in\{0,\ldots,L-1\}^{2}$, represents the probability of moving from layer $i$ to layer $j$. Subsequently, the particle diffuses within the new layer to one of its neighbours. Within each single layer, the motion of the walker is also assumed to be Markov. Hence, in the $l$-th layer it is governed by a $K\times K$ transition matrix ${\bf T}^{(l)}$, such that $T_{ij}^{(l)}$ is the probability of a particle lying in layer $l$ to diffuse from node $i$ to node $j$. These probabilities are constant over time $t$. The complete Markov model is characterised by the set of matrices $\{{\bf R}_{L},{\bf T}^{(0)},\ldots,{\bf T}^{(L-1)}\}$ and we denote ${\mathcal{T}}_{L}=\{{\bf T}^{(l)}\}_{l=0}^{L-1}$ in the sequel for convenience.

We can think of this model as a discrete-time state-space dynamical system, where the state variables at time $t$ are

[TABLE]

and we assume there is a sequence of observations $\{Y(t)\}_{t\geq 1}$ taking values in the set of node labels $\{0,1,\ldots,K-1\}$. In the main text we assume that the observations are exact and, therefore, $Y(t)=X(t)$. However, we can handle a more general class of models in which $Y(t)$ is a r.v. with conditional pmf $P(Y(t)|X(t))$ (independently of the current or past layers). If the observation is exact, then

[TABLE]

however the proposed model (and related numerical methods), admit the cases in which observation errors may occur and, hence, $P(Y(t)|X(t))$ is a non-degenerate pmf. We assume there are known and independent a priori pmf’s for the node and layer at time $t=0$, $P_{0}(X(0),\ell(0))=P_{0}(X(0))P_{0}(\ell(0))$. In practical problems, the parameters ${\mathcal{T}}_{L}$ and ${\bf R}_{L}$ are unknown and we also endow them with prior pdf’s $p_{0}({\mathcal{T}}_{L},{\bf R}_{L})=p_{0}({\mathcal{T}}_{L})p_{0}({\bf R}_{L})$ with respect to (wrt) a suitable reference measure $\mu({\sf d}{\mathcal{T}}_{L}\times{\sf d}{\bf R}_{L})$. Most often, and for a general scenario, $\mu$ can be the Lebesgue measure on $\mathbb{R}^{L\times K\times K}\times\mathbb{R}^{L\times L}$, but other choices may be posible if we wish to impose constraints on ${\mathcal{T}}_{L}$ and ${\bf R}_{L}$. For the case of the network with ring-shaped layers in the main text, $\mu$ can be reduced to the Lebesgue measure on $\mathbb{R}^{L+1}$.

C.2 Bayesian model selection

Assume that we have collected a sequence of $N$ observations which we now label

[TABLE]

We wish to make a decision as to what model is the best fit for that sequence. We adopt the view that two models are different when they have a different number of layers, hence if model $A$ has $L$ layers and model $B$ has $L^{\prime}$ layers, $A=B\Leftrightarrow L=L^{\prime}$. A convenient way to tackle this problem is to model the total number of layers $L$ as a r.v., in such a way that each possible value of $L$ corresponds to a different model. If we define

• •

a prior probability mass function for $L$, say $P_{0}(L)$, for $L\in\{1,2,\ldots,L_{+}\}$, where $L_{+}$ is the maximum admissible number of layers, and

• •

a likelihood function

[TABLE]

then we can aim at computing the posterior pmf of the number of layers

[TABLE]

and choose the model according to the maximum a posteriori (MAP) criterion

[TABLE]

i.e., we choose the value of $L$ that turns out more probable given the available observations.

Expression (13) yields the optimal solution to the problem of selection the number of layers in the multiplex from a probabilistic Bayesian point of view. As discussed below, one can find alternatives to this approach in the literature on hidden Markov models (HMMs) Rabiner89 ; Ghahramani01 , however the latter suffer from a number of theoretical and practical limitations and we strongly advocate the Bayesian solution (13).

C.3 Connections with hidden Markov model estimation theory

The problem of selecting the number of layers $L$ in the multiplex model can be cast as one of selecting a hidden Markov model (HMM) where the complete state is ${\bf Z}(t)=[X(t),\ell(t)]^{\top}$ and the transition from time $t$ to time $t+1$ is governed by the (unknown) stochastic matrices ${\bf R}_{L},{\bf T}^{(0)},\ldots,{\bf T}^{(L-1)}$. Let us adapt the notation a bit in order to make it closer to the classical HMM theory. We only have partial observations of the Markov chain, namely the sequence $\mathcal{O}=Y(1:N)$, with “emission probabilities” $P\left(Y(t)|X(t)\right)$, while the sequence of layer labels $\{\ell(t)\}_{t=1}^{N}$ remains unobserved. The goal is to estimate the total number of layers $L$ from the observations $\mathcal{O}$.

The theory of HMMs has received considerable attention in the literature since the 70s, due to their application in a variety of fields, including speech processing, molecular biology, data compression or artificial intelligence. The problem of fitting a HMM, i.e., estimating its unknown parameters has been thoroughly researched Rabiner89 . The classical technique is the Baum-Welch algorithm, which is actually an instance of the expectation-maximisation (EM) method Dempster77 ; Rabiner89 . Indeed, the general EM methodology, in several forms, is the standard approach to the problem of fitting HMMs, often combined with other techniques for its implementation, such as the Viterbi algorithm, the forward-backward algorithm or the Kalman smoother (see Ghahramani01 for an excellent survey). Many of these techniques adopt the general form of an space-alternating EM algorithm where the unobserved states and the unknown parameters are iteratively estimated, one at a time. The space-alternating generalised EM (SAGE) methodology was introduced in Fessler94 and provides a common framework for many current algorithms for fitting HMMs.

However, the estimation of the number of layers in the proposed multiplex scheme does not amount to HMM fitting. Modelling $L$ as a random variable, in order to solve problem (13) we aim at computing the model posterior probabilities given the available data $\mathcal{O}$, i.e.,

[TABLE]

where $P_{0}(L)$, $L=1,2,...$, are the a priori probabilities we attribute to models with different number of layers (e.g., we may deem models with many layers less probable than simpler models with a few layers) and $P(\mathcal{O}|L)$ is the model likelihood. The latter is an integral with respect to the probability distribution of the matrix-parameters ${\mathcal{T}}_{L}=\{{\bf T}^{(0)},\ldots,{\bf T}^{(L-1)}\}$ and ${\bf R}_{L}$ for an $L$-layer system, namely

[TABLE]

which is equivalent to Eq. (12). Recall that $p_{0,T}({\mathcal{T}}_{L})$ and $p_{0,R}({\bf R}_{L})$ are a priori pdf’s w.r.t. a reference measure $\mu$. These pdf’s can be chosen differently for different values of $L$. EM methods for HMM fitting are tools to address the problem of estimating ${\mathcal{T}}_{L}$ and ${\bf R}_{L}$ via the maximisation of the parameter likelihood $P(\mathcal{O}|{\mathcal{T}}_{L},{\bf R}_{L})$ that appears in the integrand of (15).

We see from (14) and (15), however, that what we need is to be able to integrate the likelihood $P(\mathcal{O}|{\mathcal{T}}_{L},{\bf R}_{L})$, rather than maximising it. Nevertheless, most methods in the HMM literature tackle the model selection problem (in our case, selection of the number of layers $L$) by computing estimates of the parameters via the EM method and then comparing the likelihoods of the optimised parameters Ghahramani01 ; Siddiqi07 . In our setup, this means that, given two choices $L_{1}$ and $L_{2}$, we would estimate $\hat{\mathcal{T}}_{L_{1}},\hat{\bf R}_{L_{1}}$ and $\hat{\mathcal{T}}_{L_{2}},\hat{\bf R}_{L_{2}}$ (using an EM scheme to maximise $P(\mathcal{O}|{\mathcal{T}}_{L_{i}},{\bf R}_{L_{i}})$ for $i=1,2$) and then compare the likelihoods $P(\mathcal{O}|\hat{\mathcal{T}}_{L_{1}},\hat{\bf R}_{L_{1}})$ and $P(\mathcal{O}|\hat{\mathcal{T}}_{L_{2}},\hat{\bf R}_{L_{2}})$. This approach has several problems:

• •

There is no guarantee $\hat{\mathcal{T}}_{L_{1}},\hat{\bf R}_{L_{1}}$ and $\hat{\mathcal{T}}_{L_{2}},\hat{\bf R}_{L_{2}}$ are accurate estimates (e.g., they may be overfitted). It may well happen that, e.g., $\hat{\mathcal{T}}_{L_{1}},\hat{\bf R}_{L_{1}}$ are poor estimates and, hence, $P(\mathcal{O}|\hat{\mathcal{T}}_{L_{1}},\hat{\bf R}_{L_{1}})<P(\mathcal{O}|\hat{\mathcal{T}}_{L_{2}},\hat{\bf R}_{L_{2}})$, while $P(\mathcal{O}|L_{1})>P(\mathcal{O}|L_{2})$.

• •

The EM framework yields local optimisation algorithms. Even if the EM scheme converges, it may yield a local maximiser of the likelihood for $L_{1}$ and, perhaps, a global maximiser for $L_{2}$. In this case, we may have, again, that $P(\mathcal{O}|\hat{\mathcal{T}}_{L_{1}},\hat{\bf R}_{L_{1}})<P(\mathcal{O}|\hat{\mathcal{T}}_{L_{2}},\hat{\bf R}_{L_{2}})$, while $P(\mathcal{O}|L_{1})>P(\mathcal{O}|L_{2})$.

• •

Even if we manage to obtain accurate maximum likelihood estimates of ${\mathcal{T}}_{L_{1}},{\bf R}_{L_{1}}$ and ${\mathcal{T}}_{L_{2}},{\bf R}_{L_{2}}$, there is no guarantee that $P(\mathcal{O}|\hat{\mathcal{T}}_{L_{1}},\hat{\bf R}_{L_{1}})<P(\mathcal{O}|\hat{\mathcal{T}}_{L_{2}},\hat{\bf R}_{L_{2}})$ must imply $P(\mathcal{O}|L_{1})<P(\mathcal{O}|L_{2})$.

Many authors have aimed at mitigating these flaws by introducing different heuristics in the way the models to be fitted are chosen (typically, heuristics for merging and splitting candidate states, in our case candidate layers) and producing sophisticated EM parameter estimation algorithms. See Siddiqi07 for examples. This approach does not attack the core of the problem, though.

Instead, Ghahramani01 advocates Bayesian model selection as a framework to address problem (14) that automatically handles overfitting (by imposing prior probability distributions on the parameters) and the comparison of models of different complexity (by integrating over the parameters as in (15)). In Ghahramani01 , the term used for the MAP model selection method of (13) is, actually, Bayesian integration, which makes reference to the need to solve, or numerically approximate, the integral in (15). Some candidate methods to tackle this computation include:

• •

The Laplace approximation Bishop06 , which consists in searching the maximum of $P(\mathcal{O}|\hat{\mathcal{T}}_{L},\hat{\bf R}_{L})$ and then approximating the integrand $P(\mathcal{O}|{\mathcal{T}}_{L},{\bf R}_{L})p_{0,T}({\mathcal{T}}_{L})p_{0,R}({\bf R}_{L})$ by a Gaussian with the adequate covariance structure. This approach ignores the fact that $P(\mathcal{O}|\hat{\mathcal{T}}_{L},\hat{\bf R}_{L})$ is, in our case, multimodal.

• •

The variational Bayes method Watanabe04 ; McGrory09 is an approximation scheme that relies on the use of surrogate probability distributions for the parameters (which need to be analytically tractable) in order to design an EM method that tackles the maximisation of the model likelihood $P(\mathcal{O}|L)$, i.e., the integral in (15). It is a relatively “inexpensive” method in terms of computational cost, comparable to classical EM-based model fitting techniques. However, as any EM scheme, it performs a local optimisation and does not guarantee an optimal solution.

• •

Deterministic integration of (15) using either deterministic regular grids on the space of the parameters $\{{\mathcal{T}}_{L},{\bf R}_{L}\}$ or specific cubature methods for some convenient family of functions Sobolev13 . While accurate, the complexity of these methods typically grows exponentially with the dimension of the parameters, hence they can be prohibitive for larger scale models. Examples for multiplex models with up to $L\geq 4$ layers are shown.

• •

Conventional Monte Carlo integration suffers from a similar complexity limitation. Classical Markov chain Monte Carlo (MCMC) samplers Gilks96 ; Robert04 could be well-suited to solve integrals with respect to the posterior pdf $p({\mathcal{T}}_{L},{\bf R}_{L}|\mathcal{O})$; however, the integral in (15) is actually the normalising constant of this posterior, which turns out to be hard to estimate via MCMC, which limits its application to model selection in general Robert04 .

The classical alternative to MCMC in Monte Carlo integration is importance sampling (IS) Robert04 . While conventional IS suffers from a problem called weight degeneracy, that translates into poor scaling with the dimension of the parameters in the integral, recently, families of much more efficient adaptive IS schemes have been introduced DelMoral06 ; Cappe08 ; Beskos14 ; Koblents15 . These techniques yield estimates of the integral in (15) in a simple way (unlike MCMC) and can potentially work in high dimensions Beskos14 .

Below, we present a detailed description of the nonlinear population Monte Carlo (PMC) scheme of Koblents15 and show and example of model selection with up to 10 layers ($L\leq 10$). While conventional (and even state-of-the-art) importance samplers are based on the computation of weights of the form $w(z)\propto\frac{p(z)}{q(z)}$, where $p(z)$ is the target pdf and $q(z)$ is a proposal density, the key feature of the nonlinear PMC scheme is to compute transformed weights $\bar{w}(z)\propto\phi\left(\frac{p(z)}{q(z)}\right)$, where $\phi(\cdot)$ is a nonlinear function, in order to reduce the variance of the weights (if $z$ is a random variable, then $U=w(Z)$ is random as well). This very simple transformation, if properly chosen, improves significantly the numerical stability of the algorithm when the dimension of $Z$ grows, while preserving the convergence properties of conventional IS. The examples presented below, for the nonlinear PMC and a deterministic scheme based on regular grids, show that this Monte Carlo integration scheme can be as effective as a deterministic integrator with just a fraction of the running time.

C.4 Computation of the posterior probabilities via Monte Carlo integration

Let us return to the original notation where $Y(1:N)$ denotes de sequence of observations. In order to select the number of layers $L$ in the multiplex according to the Bayesian criterion in Eq. (13), we need the ability to evaluate the posterior probability

[TABLE]

where the prior $P_{0}(L)$ is known (chosen by design) but the model likelihood $P\left(Y(1:N)|L\right)$ is an integral given by Eq. (12), namely

[TABLE]

Using the Bayes theorem, we realise that the integrand in (16) is proportional to the posterior density of the parameters given the observations $Y(1:N)$, i.e.,

[TABLE]

Taken together, Eqs. (16) and (17) indicate that the model likelihood $P\left(Y(1:N)|L\right)$ is the normalisation constant of the posterior pdf of the parameters, $p\left({\mathcal{T}}_{L},{\bf R}_{L}|Y(1:N)\right)$. This normalisation constant is often termed the model evidence is the Bayesian terminology.

An efficient way of computing the normalisation constant of a target pdf via Monte Carlo integration is by using the importance sampling (IS) method.

C.4.1 Importance sampling in a nutshell

Let $p(z)$ be a target pdf that we can evaluate up to a normalisation constant $c$, i.e., we have the ability to compute

[TABLE]

point-wise, but $c$ is unknown. The IS method Robert04 enables the estimation of $c$ (actually, it enables the estimation of integrals of the form $\int f(z)p(z)dz$ in general, for any integrable test function $f$) by sampling from an alternative pdf, $q(z)$, often called proposal density or importance function. We assume that $q(z)$ is chosen to satisfy that

[TABLE]

where $w(z)$ is the weight function. The inequality in (18) typically implies, at least, that $q(z)>0$ whenever $p(z)>0$.

The basic IS algorithm proceeds as follows:

Draw $M$ independent samples $z^{1},\ldots,z^{M}$ from $q(z)$. 2. 2.

Compute weights

[TABLE] 3. 3.

Normalise the weights,

[TABLE]

It is a straightforward application of the strong law of large numbers Robert04 to prove that

[TABLE]

for any square-integrable test function $f$.

However, the most relevant result for the purpose of this paper is that

[TABLE]

is an unbiased, consistent estimator of the normalisation constant $c$ since. By the strong law of large numbers again,

[TABLE]

since $\tilde{p}(z)=cp(z)$ and $p(z)$ is a pdf (hence, it integrates to 1).

In the Bayesian model selection problem at hand, the target non-normalised function is given by $\tilde{p}(z)\equiv P\left(Y(1:N)|{\mathcal{T}}_{L},{\bf R}_{L}\right)p_{0}({\mathcal{T}}_{L})p_{0}({\bf R}_{L})$, which we can evaluate (as will be shown below), and $c=\int\tilde{p}(z)dz\equiv P\left(Y(1:N)|L\right)$ is the model likelihood.

C.4.2 Nonlinear population Monte Carlo

The main drawback of standard IS is that, whenever there is a significant mismatch between $\tilde{p}(z)$ and $q(z)$, the variance of the weights becomes very large. As a consequence, estimators converge very slowly (with $M$) and they become of little use. This issue is usually referred to as weight degeneracy Kong94 . It typically happens when the dimension of the random variable $Z$ is large or when, simply, the target pdf is very narrow.

To mitigate degeneracy, a number of adaptive IS have been proposed, especially since the publication of Cappe04 . Here we resort to one such method, that introduces a nonlinear transformation of the weights to control their variability and, hence, degeneracy. The technique is called nonlinear population Monte Carlo (NPMC) and it was originally proposed in Koblents15 . It consists of $J$ iterations, each involving the computation of both conventional importance weights (IWs) and transformed importance weights (TIWs). The transformation is a clipping or truncation operation, denoted $\phi(\cdot,\cdot)$. For a set of $M$ ordered IWs, $\tilde{w}^{i_{1}}>\tilde{w}^{i_{2}}>\cdots>\tilde{w}^{i_{M}}$, we obtain a set of $M$ TIWs, with clipping of order $M_{c}<\sqrt{M}$, as

[TABLE]

This operation truncates the $M_{c}$ bigger weights. A general algorithm is outlined below, with $M$ samples per iteration and an instrumental Markov sampling kernel $K\left(\cdot,\cdot|{\mathcal{T}}_{L}^{i},{\bf R}_{L}^{i}\right)$ centred at ${\mathcal{T}}_{L}^{i},{\bf R}_{L}^{i}$.

Initialisation.

1. (a)

Draw $M$ independent samples $(\tilde{\mathcal{T}}_{L,0}^{i},\tilde{\bf R}_{L,0}^{i})$, $i=1,...,M$, from the prior pdf’s $p_{0}({\mathcal{T}}_{L})$ and $p_{0}({\bf R}_{L})$. 2. (b)

Compute non-normalised IWs, $\tilde{w}_{0}^{i}=P\left(Y(1:N)|\tilde{\mathcal{T}}_{L,0}^{i},\tilde{\bf R}_{L,0}^{i}\right)$, $i=1,...,M$. 3. (c)

Compute non-normalised TIWs, $\bar{w}_{0}^{i}=\phi\left(i,\{\tilde{w}_{0}^{m}\}_{m=1}^{M}\right)$, $i=1,...,M$. 4. (d)

Normalise the TIWs,

[TABLE] 5. (e)

Resample $M$ times the set $\{\tilde{\mathcal{T}}_{L,0}^{i},\tilde{\bf R}_{L,0}^{i}\}_{i=1}^{M}$, with replacement and using the normalised TIWs as probability masses, to yield an unweighted sample set $\{{\mathcal{T}}_{L,0}^{i},{\bf R}_{L,0}^{i}\}_{i=1}^{M}$. 2. 2.

Iteration. For $j=1:J$:

1. (a)

Draw $M$ independent samples

[TABLE] 2. (b)

Compute non-normalised IWs,

[TABLE] 3. (c)

Compute non-normalised TIWs, $\bar{w}_{j}^{i}=\phi\left(i,\{\tilde{w}_{j}^{m}\}_{m=1}^{M}\right)$, $i=1,...,M$. 4. (d)

Normalise the TIWs,

[TABLE] 5. (e)

Resample $M$ times the set $\{\tilde{\mathcal{T}}_{L,j}^{i},\tilde{\bf R}_{L,j}^{i}\}_{i=1}^{M}$, with replacement and using the normalised TIWs as probability masses, to yield an unweighted sample set $\{{\mathcal{T}}_{L,j}^{i},{\bf R}_{L,j}^{i}\}_{i=1}^{M}$.

After the $J$-th iteration, we have the IS estimator of the model likelihood

[TABLE]

This is a consistent estimator that converges with optimal Monte Carlo error rates. In particular, assuming that the IWs are bounded and bounded away from zero, (Koblents16, , Theorem 1) states that for any, arbitrarily small $\epsilon<\frac{1}{2}$ there exists an a.s. finite random variable $U^{\epsilon}$ such that

[TABLE]

and, therefore,

[TABLE]

All that remains is to show that the parameter likelihood function $P\left(Y(1:N)|{\mathcal{T}}_{L,j},{\bf R}_{L,j}\right)$, and hence the IWs and the TIWs, can be evaluated exactly.

C.4.3 Exact calculation of the parameter likelihood $P\left(Y(1:N)|{\mathcal{T}}_{L},{\bf R}_{L}\right)$

We start with the factorisation

[TABLE]

where each factor

[TABLE]

is a pmf that can be computed exactly using a recursive algorithm. We can write the function $P\left(Y(t)|Y(1:t-1),{\mathcal{T}}_{L},{\bf R}_{L}\right)$ as a marginal of the mass $P\left(Y(t),X(t),\ell(t)|Y(1:t-1),{\mathcal{T}}_{L},{\bf R}_{L}\right)$, namely

[TABLE]

where we have used the fact that $P(Y(t)|{\bf X}(t))=P(Y(t)|X(t))$. The conditional observation pmf $P(Y(t)|X(t))$ is one of the building blocks of the state space model, so it can be readily evaluated. The predictive pmf $P(X(t),\ell(t)|Y(1:t-1),{\mathcal{T}}_{L},{\bf R}_{L})$, on the other hand, can be decomposed as

[TABLE]

which depends on the filtering pmf at time $t-1$:

[TABLE]

Note that Eqs. (23), (24) and (25) are recursively related. If we start from the prior pmf $P_{0}(X(0),\ell(0))=P_{0}(X(0))P_{0}(\ell(0))$, then we can recursively compute the required sequence of pmf’s for $t=1,2,\ldots,T$ as

[TABLE]

Substituting (27) into (21) for each $t=1,...,N$ yields the parameter likelihood.

Appendix D QUANTIFYING COMPLEXITY: ADDITIONAL RESULTS

In this section we present some additional computer simulation results that complement those shown in the main text. For the simulations, we assume (the same as in the main text) a simple ring topology for the layers of the multiplex, as well as a particular structure for the matrix ${\bf R}_{L}$, in such a way that it can be parametrised by a single scalar $r\in(0,1)$. This is described in Section D.1. In Section D.2 we show how, given the ring topology of the layer and assuming exact observations $Y(t)=X(t)$, it is possible to construct a relatively simple deterministic approximation to the posterior probabilities $P(L|Y(1:N))$. An NPMC algorithm for this specific model is also detailed. Finally, our numerical results are displayed and discussed in Section D.5.

D.1 Multiplex of ring-shaped layers

From now on we assume a ring topology (cycle graph) for every layer. In this particular case the particle may jump from one node to its neighbours (left or right on the ring) and from one layer to another layer. This motion is random, and we model it by way of the following probabilities,

[TABLE]

where $R_{ij}$ represents the probability of moving from layer $i$ to layer $j$ and ${\sf T}^{(l)}$ is the probability of moving rightwards in layer $l$, i.e., of jumping from node $i$ to node $i+1\ (\text{mod}\ K)$. The probability of moving leftwards within layer $l$ is, therefore,

[TABLE]

These rightwards-jump probabilities for layers $l=0,1,\ldots,L-1$, are put together in the vector 444In the main text, we use ${\bf T}_{L}$ to refer to the set of transition matrices ${\mathcal{T}}_{L}=\{{\bf T}^{(l)}\}_{l=0}^{L-1}$. For the multiplex network being described, where each layer is a ring, the probabilities ${\sf T}^{(l)}$, $l=0,\ldots,L-1$, contain the same information as ${\mathcal{T}}_{L}$, hence ${\mathcal{T}}_{L}$ and ${\bf T}_{L}$ are equivalent.

[TABLE]

In order to simplify the structure of the stochastic matrix ${\bf R}$ we assume that, at each time $t$, the walker may stay at the same layer (as in $t-1$) with probability $1-r$ ($r\in(0,1)$) or jump to a different layer with probability $\frac{r}{L-1}$. Hence, the transition matrix ${\bf R}_{L}$ becomes

[TABLE]

and there is a single unknown parameter, the probability to jump, $r$.

D.2 Approximation of the posterior model probabilities

In the main text we have focused on scenarios where the node visited at time $t$ can be observed exactly, i.e., $Y(t)=X(t)$ and

[TABLE]

This constraint simplifies some of the general equations in Section C.2. Specifically, the recursions given by eqs. (26)–(LABEL:eqKK3) for the evaluation of the pmf $P(Y(t)|Y(1:t-1),{\bf T}_{L},r)$ (note that we have reduced ${\bf R}_{L}$ and ${\mathcal{T}}_{L}$ to the scalar $r$ and vector, ${\bf T}_{L}$) become

[TABLE]

Eqs. (32) and (LABEL:eqRec2) can be applied recursively with initial conditions given by the prior pmf’s $P_{0}(X(0))$ and $P_{0}(\ell(0))$. In particular, note that

[TABLE]

while

[TABLE]

These priors are assumed known as part of the model specification (typically, they can be uniform pmf’s, as they have been selected in the computer experiments).

Given Eq. (32) we can evaluate the likelihood of the parameters ${\bf T}_{L}$ and $r$. In particular,

[TABLE]

with $P(X(1)|{\bf T}_{L},r)$ computed as in (34). Numerically, it is convenient to work with the log-likelihood

[TABLE]

and transform it into natural units only when strictly needed.

The likelihood of the r.v. $L$ given by the observation record $X(1:N)$ is, therefore, given by the integral

[TABLE]

where we have assumed prior densities wrt the Lebesgue measure. Indeed, in the absence of any prior knowledge, it is natural to choose $p_{0}(r)$ to be uniform over the interval $(0,1)$, i.e., $p_{0}(r)=1$. The prior $p_{0}({\bf T}_{L})$ is important because it is used to penalise system configurations with two or more identical layers. Note that a multiplex with $L=2$ and identical transition matrices ${\bf T}^{(1)}={\bf T}^{(2)}=\tilde{\bf T}$ is equivalent to a monoplex with transition matrix $\tilde{\bf T}$. Similarly, any multiplex with $L$ layers out of which $L^{\prime}\leq L$ have identical transition matrices is fully equivalent to a reduced system with $L-L^{\prime}+1$ layers 555Along these lines, also note that a multiplex with $L$ ring-shaped layers and $r=\frac{1}{2}$ can be reduced to a monoplex with probability to move rightwards $T_{i,i+1}=\frac{1}{L}\sum_{l=0}^{L-1}{\sf T}^{(l)}$.. For our computer simulations, we have chosen

[TABLE]

i.e., we penalise systems with at least two layers that have similar transition probabilities.

D.3 Deterministic integration

One conceptually simple way to approximate the integral in Eq. (36) is

[TABLE]

where ${\mathcal{G}}=\{g_{1},\ldots,g_{H}\}$ is a grid of $H$ points over the interval $(0,1)$, i.e., $0<g_{1}<g_{2}<\cdots<g_{H}<1$. The proportionality constant of the prior $p_{0}({\bf T}_{L})$ can be approximated numerically as well (over the grid ${\mathcal{G}}^{L}$). We have selected a uniform prior on $L$, i.e., assuming that $L\in\{1,2,\ldots,L_{+}\}$, we choose $P_{0}(L)=\frac{1}{L_{+}}$. In this way, the MAP selection criterion reduces to choosing the value of $L\in\{1,2,\ldots,L_{+}\}$ that maximises the model likelihood $P(X(1:N)|L)$.

D.4 NPMC algorithm for Monte Carlo integration

Let us specify a practical NPMC algorithm for the multiplex composed of ring-shaped layers. The prior densities are $p_{0}(r)=U(0,1)$ and $p_{0}({\bf T}_{L})\propto\min_{l\neq l^{\prime}}|{\sf T}^{(l)}-{\sf T}^{(l^{\prime})}|$. We approximate the normalisation constant $\hat{C}$ of the latter pdf by standard Monte Carlo integration (simulating the ${\sf T}^{(l)}$’s uniformly over $(0,1)$).

The sampling kernels $K(\cdot,\cdot|{\bf T}_{L},r)$ are truncated Gaussian pdf’s. Specifically, if $TN(x|\mu,\sigma^{2},a,b)$ denotes the truncated Gaussian density

[TABLE]

then

[TABLE]

where $\sigma^{2}=0.005$.

The number of clipped weights is $M_{c}=\lfloor\sqrt{M}\rfloor$. Recall that function $\phi(i,\{\tilde{w}_{j}^{m}\}_{m=1}^{M})$ truncates the $M_{c}$ biggest weights to be equal (and identical to the $M_{c}$-th largest weight).

Initialisation.

1. (a)

Draw $M$ independent samples $(\tilde{\bf T}_{L,0}^{i},\tilde{r}_{0}^{i})$, $i=1,...,M$, from the prior pdf’s $p_{0}({\bf T}_{L})$ and $p_{0}(r)=1$ for $r\in(0,1)$. 2. (b)

Compute non-normalised IWs, $\tilde{w}_{0}^{i}=P\left(X(1:N)|\tilde{\bf T}_{L,0}^{i},\tilde{r}_{0}^{i}\right)$, $i=1,...,M$. 3. (c)

Compute non-normalised TIWs, $\bar{w}_{0}^{i}=\phi\left(i,\{\tilde{w}_{0}^{m}\}_{m=1}^{M}\right)$, $i=1,...,M$. 4. (d)

Normalise the TIWs,

[TABLE] 5. (e)

Resample $M$ times the set $\{\tilde{\bf T}_{L,0}^{i},\tilde{r}_{0}^{i}\}_{i=1}^{M}$, with replacement and using the normalised TIWs as probability masses, to yield an unweighted sample set $\{{\bf T}_{L,0}^{i},r_{0}^{i}\}_{i=1}^{M}$. 2. 2.

Iteration. For $j=1:J$:

1. (a)

Draw $M$ independent samples

[TABLE] 2. (b)

Compute non-normalised IWs,

[TABLE] 3. (c)

Compute non-normalised TIWs, $\bar{w}_{j}^{i}=\phi\left(i,\{\tilde{w}_{j}^{m}\}_{m=1}^{M}\right)$, $i=1,...,M$. 4. (d)

Normalise the TIWs,

[TABLE] 5. (e)

Resample $M$ times the set $\{\tilde{\bf T}_{L,j}^{i},\tilde{r}_{j}^{i}\}_{i=1}^{M}$, with replacement and using the normalised TIWs as probability masses, to yield an unweighted sample set $\{{\bf T}_{L,j}^{i},r_{j}^{i}\}_{i=1}^{M}$.

D.5 Numerical results

**From deterministic integration to NPMC: efficiently validating the full architecture. **We first consider a different example than the one shown in the main text, where now the true model from which the observations are generated has only $L=2$ layers, and we aim at estimating this hidden model via posterior probability maximisation. As for the case considered in the main text, the posterior probabilities are computed via the approximation in D.2, in particular using a 19-point grid ${\mathcal{G}}=\{0.05,0.10,\ldots,0.95\}$ for each unknown variable, which includes the probability to move rightwards, ${\sf T}^{(l)}$, over each layer and the additional parameter, $r$, which is the probability to jump to a different layer at each iteration. For a maximum of $L=4$ layers, this makes a grid of $19^{5}$ nodes, plus a $19^{4}$ grid for $L=3$, plus a $19^{3}$ grid for $L=2$, plus a 19-point grid for $L=1$.

In figure 18 we plot the posterior probability of four possible models with $L=1,\dots,4$ layers. It peaks out at the correct value, although less emphatically than for the case shown in the main text. The reason is that in this example, these probabilities are obtained using a record of only $N=10^{4}$ observations, while we used a larger time series of $N=2\times 10^{4}$ observations in the case reported in the main text. Nevertheless, the method correctly classifies the underlying architecture, what suggests that the estimation protocol is robust against short series size, something extremely relevant for real-world applications. Note also that the posterior probability of $L=1$ is zero up to the computer degree of accuracy. The observations $X(1:N)$, $N=10^{4}$, were generated from a model with $1-r=0.89$, ${\sf T}^{(1)}=0.14$ and ${\sf T}^{(2)}=0.56$. This is not a node of the grid ${\mathcal{G}}^{3}$ and yet the model is identified unequivocally.

In a second step, we explore this particular case with our NPMC algorithm, with $I=10$ iterations, $M=900$ Monte Carlo samples per iteration and clipping parameter $M_{c}=\sqrt{M}=30$, to approximate the posterior probabilities $P(\tilde{L}|X(1:N))$ for $\tilde{L}=1,2,3$ and $4$. We observe in the right panel of Fig. 18 that the posterior probability peaks at $\tilde{L}=L=2$, with similar values of the approximations. The runtime was $\approx 68$ minutes, while the computations for deterministic grid integration took $\approx 6,480$ minutes (a reduction by a factor of $\approx 95$) on the same computer and with the same software.

Moreover, the NPMC also allows us to efficiently estimate, given $L=2$, the most probable architecture of the multiplex, i.e. $r$, ${\sf T}^{(1)}$ and ${\sf T}^{(2)}$. As NPMC performs an importance sampling and therefore grids parameter space tighter in the regions with high probability measure, then Monte Carlo sampling is guaranteed to look at parameters with high likelihood. As a matter of fact, since the importance sampling provides a complete representation of the probability measure, no configurations with high likelihood will be left out (multimodal densities and likelihoods are approximated accurately, including all extrema, if the number of Monte Carlo samples –the points in the grid– is large enough), and in this sense this is a global optimization method, at odds with expectation-maximization, which converges locally.

In figure 19 we depict the parameters sampled by the NPMC (100 of Monte Carlo samples), showing that they are indeed very close to the ground truth. The probability densities of the parameters are shown in figure 20. Actually, since the grid itself yields a random discrete probability measure, integrals with respect to this measure converge to integrals with respect to the posterior probability distribution of the parameters. Hence, we can compute various kinds of estimates of the parameters (and their expected errors). A simple inspection of the posterior probability of each of these configurations therefore allows to estimate the full architecture, which will coincide with our ground truth.

Estimated architecture is a likelihood maximum. As discussed above, we have also confirmed that the likelihood of the parameters, $P(X(1:N)|{\bf T}_{L},r)$ (which is the key ingredient in the approximation of $P(L|X(1:N))$), is maximised at the values corresponding to the true model that generates the observations. In order to illustrate this fact further, we again consider the system under study in the main text (a model with $L=3$ layers, with $1-r=0.84$ and ${\bf T}_{L}=[0.16,0.76,0.24]^{\top}$) together with ten slightly perturbed models with $L=3$ layers and parameters $1-r^{\prime}=1-r+\Delta$ and ${\bf T}_{L}^{\prime}=[0.16+\Delta^{\prime},0.76+\Delta^{\prime\prime},0.24+\Delta^{\prime\prime\prime}]$, where $\Delta,\Delta^{\prime},\Delta^{\prime\prime},\Delta^{\prime\prime\prime}$ represent independent Gaussian perturbations of zero mean and variance $5\times 10^{-4}$. As a function of the time series size (the observation window $X(1:N)$), we have computed the likelihood of each set of parameter values, which in every case we compare with the likelihood obtained for the true parameters. In figure 21 we plot the logarithm of the ratio of these likelihoods, where a positive values indicates that the ratio is larger than one, i.e., that the likelihood of the true parameters is larger than the likelihood of the perturbed parameters. We systematically find that such is the case and, interestingly, we find that any small perturbation grows and generates a monotonically decreasing likelihood-ratio as the size of the observed series increases.

Scalability. Finally, we consider the issue of scalability and computational cost. All simulations have been carried out using Matlab R2016a running on an Apple iMac equipped with an Intel Core i7 processor, with 4 cores and 32 GB of RAM. For a first experiment, we have considered the same model and data as in Fig. 18, i.e., a system with $L=2$ layers, parameters $1-r=0.89$, ${\sf T}^{(1)}=0.14$ and ${\sf T}^{(2)}=0.56$ and a series of $N=10^{4}$ observations. As discussed previously, using NPMC instead of deterministic integration reduced by a factor of $\approx 95$ computing time. Therefore, equipped with the NPMC algorithm we can tackle more complex systems. For the last computer experiment we have generated a sequence of $N=20\times 10^{3}$ observations from a model with $L=5$ layers, and parameters ${\sf T}^{(1)}=0.1$, ${\sf T}^{(2)}=0.2$, ${\sf T}^{(3)}=0.6$, ${\sf T}^{(4)}=0.7$, ${\sf T}^{(5)}=0.9$ and $1-r=0.85$. We have applied the NPMC algorithm, with $J=10$ iterations $M=1,000$ Monte Carlo samples per iteration and $M_{c}=31$, to approximate the posterior probabilities $P(\tilde{L}|X(1:N))$ for $\tilde{L}=1,\ldots,10$. This simulation took around 370 minutes on the same computer. Fig. 22 displays the results, where we observe that the correct model is still clearly identified.

Appendix E BIOCHEMICAL NETWORKS: ADDITIONAL DETAILS FOR THE CASE OF RNA POLYMERASE

In figure 23 we plot (right panels) the five experimental traces of RNA Polymerase I analysed in this work. To make the inference problem more challenging, we systematically select the initial $10^{4}$ time steps (about 10 seconds, as the signal is sampled at 1kHz), where it is often not evident to visually distinguish the periods where elongation and backtracking are at play. Due to the inherent noise inside the cell, and as one can see in the figure, only for two cases (traces $\#2$ and $\#3$) the two mechanisms operating –elongation and backtracking– are observable, while the rest are hidden below the noise. Our methodology (both the layer detection and layer estimation protocols) are succesful at detecting and correctly inferring the presence of these two modes as we find that the more likely number of hidden layers is $L=2$. Intringuingly, if we select a period where apparently only one mode (backtracking) is at play, then our methodology again predicts multiplexity ($D(3)>0$) and $L=2$. This is a surprising result which suggest that, either (i) short elongation times are hidden among a long backtracking pause, or that (ii) while in backtracking, stochastic fluctuations have a non-Markovian character. Either interpretations are in contrast with the state of the art and deserve further investigation.

Appendix F PROOF OF THEOREM 1

It is well known that a higher-order Markov model ($h>1)$ can be converted into a first order Markov system ($h=1$) defined over a higher dimensional state space (see, e.g., Chen00 ). In the same vein, we are going to show that every higher-order Markov model can be converted into an equivalent hidden jump-Markov model of the type discussed in Sec. II (which from now on we call multiplex model) with a sufficiently large number of layers. The proof hence consists of two parts:

• •

We first describe a constructive procedure to obtain a multiplex model from the transition probabilities of a generic Markov model of order $h$.

• •

We then prove that the probability to generate any sequence of length $h+1$ is the same for the (order $h$) Markov model and the multiplex model.

F.1 Construction of the multiplex model

Recall that the Markov sequence takes values on the space $\mathbb{K}=\{1,\ldots,K\}$ (this is the node set in the multiplex network) and is defined by the transition probabilities $P_{\mathbb{K}}^{h}(i_{t}|i_{t-h:t-1})$, $i_{t-h:t}\in\mathbb{K}^{h+1}$. We aim at constructing a model with $L=K^{h-1}$ layers, which implies that the random layer sequence $\ell(t)$ takes values on the set ${\mathcal{L}}_{h-1}:={1,\ldots,K^{h-1}}$. As a first step, we establish a bijection between the set of layers ${\mathcal{L}}_{h-1}$ and the set of vectors $\mathbb{K}^{h-1}$. To be specific, let us introduce the map $b:\mathbb{K}^{h-1}\rightarrow{\mathcal{L}}_{h-1}$ which associates every vector of states

[TABLE]

with a unique layer $b({\bf v}):=1+\sum_{s=1}^{h-1}(i_{s}-1)K^{s-1}\in{\mathcal{L}}_{h-1}$. Intuitively, we interpret the vector

[TABLE]

as a sequence of digits in a base $K$ number system, and the function $b(\cdot)$ symply returns the evaluation of this number (plus one, for notational coherence). Since $b$ is bijective, we can readily construct the inverse function $b^{-1}:{\mathcal{L}}_{h-1}\rightarrow\mathbb{K}^{h-1}$ that takes a layer number $\ell\in{\mathcal{L}}_{h-1}$ and returns its base-$K$ representation (in vector form) $b^{-1}(\ell)$.

With these ingredients, we can now construct a multiplex model with the form of a Markov chain on the state space $\mathbb{K}\times{\mathcal{L}}_{h-1}$ described by the transition probabilities

[TABLE]

where $i_{t},i_{t-1}\in\mathbb{K}$, $j_{t},j_{t-1}\in{\mathcal{L}}_{h-1}$ and

[TABLE]

Note that we have simply encoded the history of states $i_{t-2:t-h}$ into the layer index $j_{t-1}$. For the system to behave “correctly” we need to check that the transition from layer $j_{t-1}$ to layer $j_{t}$ is compatible with the sequence of states $i_{t-1:t-h+1}$ that has to be “recorded” for the next element of the chain. Unfortunately, the Markov chain described by the transition probabilities $\mathbb{M}_{*}(\cdot|\cdot)$ is semi-degenerate. To be specific, the layer value $\ell(t)=j_{t}$ is deterministic given $X(t-1)=i_{t-1}$ and $\ell(t-1)=j_{t-1}$. Indeed,

[TABLE]

with probability 1. This determinism is not coherent with our discussion of multiplex models in Sections I and II. Fortunately, this difficulty can be easily removed if we further extend the state space of the model.

Let us construct a (new but similar) Markov chain on the space $\mathbb{K}\times{\mathcal{L}}_{h}$, defined by the transition probabilities

[TABLE]

where

[TABLE]

and

[TABLE]

(Note that, in this case, the bijection $b$ is a function $\mathbb{K}^{h}\rightarrow{\mathcal{L}}_{h}$, with the same interpretation as before.) Intuitively, we have now encoded the system memory (the sequence of states $i_{t-h:t-2}$) into a subset of $K$ layers, instead of a single one. From each layer in this subset there is probability $\frac{1}{K}$ to jump to a new layer belonging to the subset that encodes the sequence $i_{t-h+1:t-1}$. Comparing Eq. (42) with our notation for multiplex models in Sections I and II we can readily identify the matrix ${\bf R}_{L}$ of transition probabilities between layers, while $T_{i_{t-1}i_{t}}^{(j_{t-1})}=P_{\mathbb{K}}^{h}(i_{t}|i_{t-1:t-h})$ (note that $j_{t-1}$ uniquely identifies the sequence $i_{t-2:t-h}$).

F.2 Identity of transition probabilities

What remains to be shown is that the conditional probabilities

[TABLE]

generated by the probability function $P_{\mathbb{K}}^{h}(\cdot|\cdot)$ (for the order-$h$ Markov model) and the Markov kernel $\mathbb{M}(\cdot|\cdot)$ in (42) (for the multiplex model) coincide.

On one hand, if $\{X(t)\}_{t\geq 0}$ is an order-$h$ Markov sequence, then

[TABLE]

by definition. On the other hand, if the sequence $\{X(t)\}_{t\geq 0}$ is generated by the multiplex model on $\mathbb{K}\times{\mathcal{L}}_{h}$ defined by $\mathbb{M}$, we have

[TABLE]

where the first equality follows from the theorem of total probabilities, the second equality is obtained from the definition of conditional probability and the third identity follows from the definition of the multiplex model on $\mathbb{K}\times{\mathcal{L}}_{h}$. Moreover, for the multiplex model constructed in Appendix F.1,

[TABLE]

Hence substituting (54) and (42) into (47) yields

[TABLE]

where

[TABLE]

However, from the definition of $R_{j_{t-1}j_{t}}$ in Eq. (43) we can further reduce Eq. (55) and obtain

[TABLE]

and, for each term in the sum of Eq. (56) we have $R_{j_{t-1}j_{t}}=\frac{1}{K}$ (from the definition in Eq. (43)). Because of the construction of the sets $S_{b^{-1}}(i_{t-h+1:t-1})$ and $S_{b^{-1}}(i_{t-h:t-2})$, there are exactly $K^{2}$ terms in the sum of Eq. (56), which readily yields

[TABLE]

and completes the proof. $square$

Appendix G PROOF OF THEOREM 2

Theorem 1 guarantees that any order-$h$ Markov sequence on $\mathbb{K}$ can be represented exactly by a multiplex model with $L=K^{h}$ layers. Hence, for every random sequence $\{X_{h}(t)\}_{t\geq 0}$ we have an equivalent sequence $\{X_{L}(t)\}_{t\geq 0}$ with $L=K^{h}$ layers and transition probabilities

[TABLE]

$i,j\in\{1,\ldots,L\}$ for the layers and

[TABLE]

for $i,s\in\{1,\ldots,K\}$, $j\in\{1,\ldots,L\}$. Therefore, every Markov-$\infty$ sequence $\{X(t)\}_{t\geq 0}$ can be represented as well by a sequence of random processes $\{X_{L}(t)\}_{t\geq 0}$ described by suitably constructed multiplex models, i.e.,

[TABLE]

where $\stackrel{{\scriptstyle d}}{{=}}$ denotes equality in distribution. What we need to prove is that the sequence of multiplex models, described by the pairs $({\bf T}_{L},{\bf R}_{L})$, $L=K,K^{2},K^{3},...$, yields the continuous-layer model described by a proper probability measure ${\mathcal{M}}_{0}$, a Markov kernel ${\mathcal{M}}$ on ${\mathcal{B}}([0,1))$ and a family of transition matrices ${\bf T}(y)$, $y\in[0,1)$, when $L\rightarrow\infty$.

Let us note that every transition probability in the multiplex model can be specified by a simple function Kolmogorov75 defined over subsets of $[0,1)$. To be specific, we construct

[TABLE]

where

[TABLE]

is the set-indicator function,

[TABLE]

and

[TABLE]

The limits $\lim_{L\rightarrow\infty}\hat{T}_{is}^{L}(y)$, $\lim_{L\rightarrow\infty}\hat{\mathcal{M}}^{L}(y|y^{\prime})$ and $\lim_{L\rightarrow\infty}\hat{\mathcal{M}}_{0}^{L}(y)$ are well defined (because $X_{L}(t)$ yields the same transition probabilities as $X_{h}(t)$, for $L=K^{h}$, and $\lim_{h\rightarrow\infty}X_{h}\stackrel{{\scriptstyle d}}{{=}}X$ uniformly over $t$ and $\mathbb{K}^{h+1}$). Therefore, from Kolmogorov75 (see Theorem 5 in Chapter 8) the functions

[TABLE]

exist and they are measurable w.r.t the Lebesgue measure. In particular, $\bar{\mathcal{M}}(y|y^{\prime})$ and $\bar{\mathcal{M}}_{0}(y)$ are probability density functions w.r.t. the Lebesgue measure and we obtain ${\mathcal{M}}(dy|y^{\prime})=\bar{\mathcal{M}}(y|y^{\prime})dy$, and ${\mathcal{M}}_{0}(dy)=\bar{\cal M}_{0}(y)dy$. $\Box$