Heterogeneous network flow and Petri nets characterize multilayer complex networks

Interacting subsystems are commonly described by networks, where multimodal behaviour found in most natural or engineered systems found recent extension in form of multilayer networks. Since multimodal interaction is often not dictated by network topology alone and may manifest in form of cross-layer information exchange, multilayer network flow becomes of relevant further interest. Rationale can be found in most interacting subsystems, where a form of multimodal flow across layers can be observed in e.g., chemical processes, energy networks, logistics, finance, or any other form of conversion process relying on the laws of conservation. To this end, the formal notion of heterogeneous network flow is proposed, as a multilayer flow function aligned with the theory of network flow. Furthermore, dynamic equivalence is established with the framework of Petri nets, as the baseline model of concurrent event systems. Application of the resulting multilayer Laplacian flow and flow centrality is presented, along with graph learning based inference of multilayer relationships over multimodal data. On synthetic data the proposed framework demonstrates benefits of multimodal flow derivation in critical component identification. It also displays applicability in relationship inference (learning based function approximation) on multimodal time series. On real-world data the proposed framework provides, among others, multimodal flow interpretation of U.S. economic activity, uncovering underlying empirical steady state probability distribution, as well as inherent network (economic) robustness.


Background concepts
To aid further understanding of the proposed framework, a brief introduction of background concepts is presented, including basic terminology and notation. The formal notions of graphs, network flow, multilayer networks and Petri nets are introduced individually, by recalling some relevant definitions.

Graphs.
A graph G is considered an ordered tuple G = (V , E) , with a vertex set V of cardinality n (graph order) and edge set E of cardinality m (graph size). The edge set E is a subset E ⊆ [V ] 2 , where [V ] 2 denotes all two-element subsets of V defining distinct pairs of adjacent vertices or nodes, incident to their relevant edge. An adjacency matrix A m ∈ R (n×n) is a matrix encoding among all distinct vertex pairs of a graph, those which belong to an edge set. An incidence matrix A t ∈ R (n×m) is a matrix encoding all incident vertices and edges. The graph G = (V , E ′ , µ ′ ) comprising of not necessarily distinct vertex pairs is considered a multigraph 37 . The multigraph is defined with respect to a finite edge set E ′ , where association to pairs of nodes is taking place with respect to a map µ ′ : E ′ → [V ] 2 . The map µ ′ assigns to each edge two end vertices, allowing multiple edges between adjacent nodes (note that self-loops are left out in this definition). To simplify notation, a graph is always denoted as G = (V , E) , where a multigraph is always understood in the latter context. The graph can be weighted ω : E → R + ( R + = {x ∈ R | x > 0} ) or unweighted, where ω denotes the weighting coefficient. The graph is assumed connected, or otherwise decomposed into the union of its connected graph components.
A network N is a digraph or directed graph defined on an ordered tuple N = (V , A N ) . A flow network N f is a digraph with additional structural properties i.e., N f = (N, c, s, t) ≡ (V , A N , c, s, t) . The digraph consist of an arc set A N ⊆ V × V of ordered vertex pairs or node pairs, forming directed edges or arcs (u, v) ∈ A N . Node u denotes the initial vertex or arc tail and node v the end vertex or arc head. A flow network comprises furthermore of an arc capacity function c : A N → R + . It is also defined with respect to two distinguished subsets of V, denoted as S and S , where S ⊆ V defines the set of sources of N f , and S ⊆ V the set of sinks. From a physical point of view, these subsets define the sets of points where flow enters or leaves a network, in form of exogenous flow. From a mathematical point of view, S and S are arbitrary subsets of V (normally disjoint, as any overlap is cancelled out as net flow). Vertices v ∈ V \ (S ∪S) are called intermediate. Flow network N f can, without loss of generality, be reduced to a flow circulation, a single-source s and single-sink t representation by auxiliary vertices s, t ∈ V , referred to as the environment. Single-source s is incident with all source vertices S = {u | (s, u) ∈ A N } , and such that (v, s) / ∈ A N , ∀v ∈ V . Single-sink t is incident with all sink vertices S = {v | (v, t) ∈ A N } , and such that (t, u) / ∈ A N , ∀u ∈ V . The transformation converts sources u ∈ S \ {s} and sinks v ∈S \ {t} to intermediate nodes, with intermediate arc set AN ⊆ A N obtained as the set of all arcs not incident with the environment. Arc capacity from and to the environment is assumed infinite. The direction of flow in a flow network corresponds to the direction of an arc (providing physical meaning), hence bidirectional flow is represented by two oppositely directed arcs. Network flow. Network flow is a function defined with respect to the topology of a flow network. Somewhat counter-intuitively it is not primarily concerned with the mechanics of flow, but rather with the algebraic relationships encoded in the flow process. A network flow can be considered a form of algebraic tool. The following definition is derived from 28,29 , chosen for interpretability (slightly adapted to the formal content of this paper), but can equivalently be given with respect to a flow circulation 37 . www.nature.com/scientificreports/ where f is a nonnegative flow with capacity constraint c, while (j, i), (i, j) ∈ A N correspond to inflow and outflow arcs of node i ∈ V , respectively. The definition of network flow is here introduced with respect to the notion of positive or negative direction (similar to a vector), rather than positive or negative flow (providing physical meaning and some formal convenience). It is worth noting that in a classic network flow node participation is one-fold, and corresponds to conservation of flow. In the work presented in this paper this role is extended to coupling (in form of flow conversion). It is also worth pointing out that a network here is understood as an abstract mathematical object, where network structure encodes relationships describing certain systems or processes. This will be explained in more detail in the proposed framework section, under "Heterogeneous flow networks".

Multilayer networks.
A multilayer network comprises of a node set V, as any other network or graph. It furthermore consists of layers α i , i ∈ {1, . . . , b} (Fig. 1a, horizontal rectangles), where b denotes the total number of layers. Each layer α i is an array (α i 1 , . . . , α i d ) of elementary layers α i a (Fig. 1a, Cartesian axes labels, e.g., X, Y), corresponding to aspects a ∈ {1, . . . , d} , where d denotes the total number of aspects or dimensions (Fig. 1a, Cartesian axes). The layers form a layer sequence L M = {L a } d a=1 , consisting of sets of elementary layers L a , where L a = {α | ∃i ∈ {1, . . . , b}, α i a = α} (Fig. 1a, light grey highlight). Nodes and layers form a multilayer node set V M ⊆ V × L 1 × ... × L d , comprising of node-layers (u, α i ) ≡ (u, α i 1 , . . . , α i d ) , such that u denotes a node u ∈ V existing on a respective layer α i (Fig. 1a, e.g., node-layer (2, α (Fig. 1a, solid lines), while an inter-layer edge set is a complementary set E C = E M \ E A (Fig. 1a, dashed lines). A coupling comprising of intra-layer edge set E A (solid lines), inter-layer edge set E C (dashed lines), and coupling edge set EC (dashed grey lines, corresponding to e.g., node-layer (1, α i ) , i ∈ {1, 3} ). (b) Heterogeneous (layerdisjoint) network, with path between nodes in layers α 1 , α b (left), along with network schema and meta-path (c) Petri net, with place nodes p i , p j ∈ P (circles) interpretable as resources or system states, and transition nodes t k ∈ T (bars) interpretable as processes, along with Pre and Pos functions (arcs), interpretable as conversion ratios or weights (left) (e.g., 1 unit of resource p 1 , produces 1.2 units of resource p 2 , and 1.3 units of resource p 3 , through process t 1 (right)). www.nature.com/scientificreports/ edge set is a subset EC ⊆ E C , EC = {{(u, α i ), (v, α j )} ∈ E C | u = v} (Fig. 1a, dashed grey lines, corresponding to e.g., node-layer (1, α i ) , i ∈ {1, 3} ). One can analogously define the multilayer arc set is referred to as the underlying multilayer network graph. If the graph is weighted, a weighting function ω M : E M → R + is defined. The graph is otherwise referred to as unweighted. The concept of a graph or network introduced in previous subsections is what is commonly referred to as a single-layer network or monoplex in multilayer network terms. The formal definition of the presented multilayer network concept is introduced in Definition 2. A graphical illustration is presented in Fig. 1a.
Definition 2 (Multilayer network 1,2 ) Let L M be a layer sequence L M = {L a } d a=1 of sets of elementary layers L a for each aspect a ∈ {1, . . . , d} , and let V M be a multilayer node set In a heterogeneous network all nodes associate to one specific layer or type, and can be adjacent to nodes of either the same or a different type 19 , which also applies to edges 22 . The network is formally classified as layerdisjoint, which emphasizes that each node of the network associates to a single layer only 2 . In the following definitions the formal notions of a layer-disjoint 2 and heterogeneous network 19,22 are recalled. The corresponding notions of a meta-path and information network are interpreted from 22 .
Definition 4 (Information network 22 ) Let N be a digraph N = (V , A N ) , with an object type mapping function θ : V → O T and a relation type mapping function ψ : A N → R T . An information network I is a digraph defined on an ordered tuple I = (N, θ, ψ) , such that there is a surjection ψ(l) = (θ(u), θ(v)) , where each vertex u, v ∈ V belongs to one particular object type θ(u), θ(v) ∈ O T , and each arc l = (u, v) ∈ A N belongs to one particular relation type ψ(l) ∈ R T . The digraph T I = (O T , R T ) is a network schema, defined as a meta template of information network I.
An information network represents a mathematical object, establishing correspondence between a node set V and arc set A N , and a set of object types O T and relation types R T , respectively, such that reference to each element is preserved. The formal notion enables introduction of a heterogeneous network (Fig. 1b, left panel), where nodes and edges associate to one specific layer or type (Fig. 1b, right panel), such that there is more than one object type O T or relation type R T , as introduced next (Definition 5). The notion of heterogeneity refers here to semantics (i.e., physical interpretation of disjoint layers), while other types of heterogeneity (e.g., node degree) are supported implicitly, as is any other property inherent to single-layer (classic) networks, as supported by the generalized network form, the multilayer network (where a multilayer network represents a generalization over a classic definition of a network, while a heterogeneous network represents a multilayer network class).
A network schema of a heterogeneous network corresponds to a topological projection of paths between node-layers (Fig. 1b, left panel), onto a set of composite relations between elements of object types (Fig. 1b, right  panel). The formal notion of a universal path across object types (Fig. 1b, right panel) corresponds to a metapath 22 . Note that an inverse relation R −1 l on a meta-path does not necessarily always exist in a general sense. A transpose relation R T l , however, does by definition.
Petri nets. The Petri net 30,31 is a model of information flow, with specific focus on concurrent event systems and topological representation of the underlying relationship structure. A Petri net is a tuple of the form P N = (P, T, Pre, Pos) , where: P is a set of n place nodes, T is a set of m transition nodes, Pre : P × T → R + 0 is an incidence function specifying weights from places to transitions, and Pos : P × T → R + 0 is an incidence function specifying weights from transitions to places (Fig. 1c). In a physical sense, place nodes (Fig. 1c, circles) can informally be understood as resources or system states, while transition nodes (Fig. 1c, bars) can be understood as processes. The incidence matrix is derived as C = Pos − Pre . The pre and post sets are denoted as • u, u • , respectively, where u ∈ P × T . A marking is a mapping of the form Q : P → R + 0 (that can be represented by a www.nature.com/scientificreports/ vector, once places are ordered), which assigns to each place a nonnegative quantity called mark. A P-invariant is a vector of the form x ∈ R + n 0 , where x T C = 0 with support set �x� = {p i ∈ P|x i > 0} . A T-invariant is a vector of the form y ∈ R + m 0 , where Cy = 0 with support set �y� = {t j ∈ T|y j > 0} . The two relations encode fundamental laws of mass conservation (P-invariant), and flow balance (T-invariant) 38 .
The Petri net was initially introduced as a discrete model 30,31 , but was later on extended to a continuous Petri net (CPN) 32,33 P N , Q 0 , , where: Q 0 is a vector of initial markings Q 0 i = Q 0 i (p i ) , assigned to each place p i ∈ P , |P| = n , and = [ 1 , 2 , . . . , m ] T is a vector of finite transition firing rates, assigned to each transition t j ∈ T , |T| = m . The mark enters a place p i at time τ and enables all adjacent transitions t j ∈ p • i , as soon as the marking Q i becomes available. The marking evolution is governed by the state equation such that transition t j ∈ T is considered enabled, provided that one can find a rate vector ∈ R + m 0 , with t j ∈ � � , where [C ] i ≥ 0 for all i such that Q i = 0 . The role of flow in a CPN is directly related to the notion of a stationary transition firing rate . The set of admissible firing rates at equilibrium can be computed by means of a CPN linear program 35 , which for a system of the form P N , Q 0 , , with incidence matrix C and a finite firing rate , can be formulated as

Proposed framework: heterogeneous network flow (multimodal flow)
The proposed framework is introduced with respect to two novel complementary concepts, namely: the formal notion of a heterogeneous flow network, as the underlying network topology and multilayer equivalent of a Petri net, and heterogeneous network flow, as the multilayer flow function enabling derivation of multimodal flow. The framework is thereby aligned with formal requirements from complex network theory and the theory of network flow, introduced in "Background concepts".

Heterogeneous flow network. A heterogeneous flow network (HFN) is a class of layer-disjoint multilayer
flow networks, endowed with structural properties aligned with the theory of network flow (Definition 1). For tractability and with slight abuse of notation, node-layers (u, α i ) are from here on denoted as u i , given layerdisjoint nature of an HFN, which preserves the layer reference of node u. In the following definition the formal notion of an HFN is introduced.
The HFN is layer-disjoint in the sense that a node residing in multiple layers of an HFN has parts physically belonging to different semantic domains. The correspondence can be linked to the notion of a split node 29 , a network transformation converting one node u into b distinct nodes u i , connected by objects referred to as links, (Fig. 2a, illustrating a heterogeneous flow network H f , with node u split into two distinct nodes, u and ũ , forming link (u,ũ) ). In an HFN a link corresponds to an inter-layer arc (where AC = ∅ i.e., the nodes reside in distinct layers and are layer-disjoint). The link constitutes an element of the composite relationship structure formed between nodes in different layers i.e., the HFN meta-path. The definition of an HFN meta-path is taken over from Definition 6. The link is directed due to the directed nature of a relation. In topological line-graph terms (graph dual terms), a link is defined with respect to the notion of a node. The node in an HFN meets two functional properties, namely: conservation of flow, in its original network flow sense, and coupling, in its line-graph form as a link (through flow conversion). This goes beyond the classic definition of a network flow, where a flow function does not distinguish on the role of coupling or concurrency. The HFN network flow, however, does due to the following property (Definition 8): the line-graph equivalent of a heterogeneous multilayer flow network is a Petri net (Fig. 2b). (ii) Transition nodes: each corresponding to a single intra-layer arc (u, v) ∈ A A ; (iii) Incidence functions: www.nature.com/scientificreports/ Finally, we associate a function ¯ : T → R + , which assigns to each transition a firing rate capacity bound equivalent to c, and a function ω : T → R + , which assigns to each transition a weight equivalent to ω , such that  (7), (8) such observation is trivial, given one-to-one correspondence between all elements of domain and range, for each relation individually. Though the single-source and single-sink sets, s and t, are not mapped to the set of all places P, the transformation is invertible given relation (4). The relation maps the entire intra-layer arc set A A to a set of all transitions T, making s and t the complementary node sets, relating to all transitions with no incoming or outgoing places, respectively, which is uniquely determined. For Eqs. (5) and (6) the transformation is not one-to-one. However, note that A C is a set of all inter-layer couplings, uniquely determined by its node set counterpart. Inter-layer links (u,ũ) ∈ A C correspond to reference of node u to each distinct layer α the node resides on (Fig. 2b,d), as retained by labelling of node-layers in V M , and consequently places in P, which can be traced back and is uniquely determined.
The transformation in Definition 8 is invertible, within the class of Petri nets obtained by the injection (Remark 1). It is retained by the one-to-one correspondence within its image (Fig. 2b). Note that a multi-semantic split node (node residing in multiple layers) or its links, correspond to a single network element, transformed to a transition node with multiple outgoing and/or incoming arrows (as an equivalent network object). Place and transition nodes can thereby, as specified in "Background concepts", be interpreted as e.g., resources and processes, respectively.
The transformation encodes through Eq. (1), the evolution process of a positive compartmental system q(φ) 39,40 , where q denotes donor mass (equivalent to Q), and φ the flux or mass flow (equivalent to ). Incidence matrix C, derived from Eqs. (5) and (6), encodes furthermore, a proportion between adjacent domains α i , α j , φ j = kφ i41,42 , where k denotes a coupling modulus or conversion coefficient. It is worth noting that a Petri net concurrently represents both, conservation and coupling relations, structurally encoded in incidence matrix C. A Petri net can thereby be understood as a form of mathematical object, where state space relations and physical constraints are integrated into the network structure (as in the case of a classic graph or flow network) (Fig. 2c). This property, beyond computational context (transition firing rates, marking evolution process 35 ) exploited in the following subsection, provides access to a number of structural and behavioural analysis tools 38 .

Heterogeneous network flow.
A heterogeneous network flow is a flow function, which can be described by two types of flow processes, namely: transfer, a classic form of intra-layer network flow, and coupling, a form of inter-layer flow establishing exogenous relationship between layers in form of transformation. Though coupling may conveniently be referred to as inter-layer flow, in essence, no physical flow as in the sense of transfer is taking place. The flow represents rather an exchange or update, where relationship between inflow and outflow is not necessarily an identity (as in the case of intra-layer flow). Wherever inter-layer flow is encountered, a type of flow conversion is taking place, in line with respective network semantics (Definition 9). (i) Capacity constraint: Heterogeneous flow network H f and split node u, split into two distinct nodes u and ũ , residing in two layers α , α , and forming link (u,ũ) , which acts as exogenous arc with respect to the sink layer α (mapped to single-source s as flow circulation). (b) Transformation of heterogeneous flow network (left-hand side), to line-graph Petri net equivalent (right-hand side), where arcs (lines) match transition nodes (bars), while nodes (dots) match places (circles) with corresponding end node transitions. Single-source and single-sink arcs, (s, u) and (v, t) , are mapped as transition nodes only, corresponding to sources and sinks of exogenous flow. (c) Breakdown of linegraph Petri net equivalent of network in panel a, with corresponding set of relations, where highlighted objects encode: conservation of mass (left-hand side), and flow balance (right-hand side). (d) Illustration of incidence matrix transformation, corresponding to networks in panels a and c, where A t corresponds to graph incidence matrix of network in panel a (unweighted, except for link (u,ũ) ), while C corresponds to Petri net incidence (coupling) matrix of network in panel c. The relationship can equivalently be written as:  where (ṽ,ũ) , (ũ,ṽ) are inflow and outflow arcs of node ũ in layer α , respectively, (u, v) are outflow arcs of node u in layer α , and (u,ũ) ∈ A C are coupling links of weight k, participating in flow conversion Relations (9) Remark 2 For a flow network H f associated to Petri net P f by transformation T f : H f → P f and relation (12), the following statements hold: (i) f h fulfils condition (9) iff 0 ≤ ≤¯ , following from the definition of firing rate capacity bound (Eq. 7); (ii) f h fulfils condition (10) iff the corresponding vector fulfils C = 0 ; (iii) if there exists a firing rate such that conditions (9) and (10) are satisfied, then there exists a firing rate ˜ such that conditions (9)-(11) are satisfied as well, and vice versa, as fulfilled by the standard construction where q(f h ) corresponds to the left-hand side of relation (10), are solutions of Q = C given Eq. (12), where q relates to instantaneous flow balance at each compartment ũ of a positive compartmental system 39 , with exogenous flow determined by all terms relating to arc sets A C and A A \ AĀ , while remaining flow balance is determined by all residual terms. The relationship proves HFN-PN dynamic equivalence.
From Theorem 1 and Remark 2, it can be deduced that stationary firing rate of a Petri net directly corresponds to multilayer network flow f h of an HFN. Incidence matrix C is bidirectional, with direction encoded in oppositely facing transitions (Fig. 2e). The link coupling direction may thereby exhibit two distinct forms, where positive or negative flow conversion (Fig. 2f) encodes the relationship of a multimodal node with each layer (as a form of source or sink assignment, respectively). Transpose flow (Fig. 2g), on the other hand, refers to a complementary relationship between participating layers, in the sense of the notion of a transpose relation. Incidence matrix C ultimately encodes relationships between nodes in different layers, with network flow acting as a form of algebraic tool, where time complexity of relation (12) rests on matrix-vector multiplication and is at most O(nm) , such that n is the graph order and m the graph size.
Overall, the key property of the proposed framework lies in the ability to derive multimodal flows across different types of network layers, in a compact unified form. Inter-layer flow conversion takes place across network links, which act as exogenous control flows with respect to the sink layer. The multimodal flows are implicitly either coupled (across layers) or uncoupled (in intra-layer flow derivation, specifically with regard to cycle space conditions, as presented next). This property lends some new perspective to network interdependency assessment, such as e.g., with respect to flow based centrality 43-45 , system stability or reachability of states 38 . What is important to note is that the proposed framework satisfies concurrently conditions of conservation and coupling www.nature.com/scientificreports/ of flow across network semantics. This is in contrast to current state-of-the-art, as presented in e.g. 26,27,46 , where the main processes studied to date relate to epidemic spreading and single-mode information diffusion 8,47,48 . This renders the proposed framework more applicable to certain other types of multilayer network problems, as presented next.

Illustrative examples: synthetic and real-world networks
To demonstrate application potential and illustrate implementation of the proposed framework, different examples of multilayer heterogeneous network flow are presented. The varied range of examples is thereby aimed at displaying generality and breadth of framework application. Further extensions related to domain specific applications are introduced as well.
Laplacian flow. To demonstrate implementation of the proposed framework in uncovering importance of specific network components, an illustrative example of flow based random walk node betweenness centrality C RW 44 is presented. The Petri net flow relations (statement (i)-(iv), Remark 2) are here extended, to apply to a broader range of application domains, such as random walks, which correspond to a simple type of linear dynamic flow, also known as Laplacian flow 39,49 . The Petri net flow relations are thereby complemented by node potential balance (cycle space) conditions (along the lines of the notion of fundamental cycles over states 38  where q is a donor mass vector and a compartmental matrix which is Metzler i.e., It is worth pointing out that in a random walk arc weights are symmetric i.e., ω(u, v) = ω(v, u) for each arc. The proposed framework, however, supports asymmetry in arc weights as well i.e., ω(u, v) = ω(v, u) , which is a more general case. Computational complexity of relation (13) is dictated by Gaussian elimination in producing the basis of the incidence matrix null-space, which has an arithmetic complexity of O(n 3 ) , where n is the graph order, while time complexity is algorithm dependent and is polynomial (using standard algorithm forms). Overall, the proposed framework introduces a form of multilayer Laplacian flow. By comparison, state-of-the-art Laplacian flow applied to multilayer network structures 46 does not capture the notion of coupling (conversion) between different semantic domains. It is hence not applicable to many physical problems which rely on laws of flow and mass conservation (e.g., energy networks, chemical processes, financial transactions, etc.), compared to the proposed framework, as demonstrated in the following example.
The multilayer network selected for the illustrative example (Fig. 3a), may in a broad sense represent a resource shipped from source s ′ to sink t ′ of a source layer α ′ , according to an arc weight ω (interpretable also as transition probability, arc preference, or conductance). The shifting of the resource may generate further outputs in form of e.g., products, revenues, mechanical force or similar, with respect to the sink layer α ′′ , obtained based on a conversion coefficient k. For the sake of generality, one or both layers may be assumed to meet (Eq. 13). The approach is thereby general enough to represent a number of equivalent problem formulations. In this particular case the object in Fig. 3a represents an energy network, comprising of a three tank energy resource layer α ′ (node 1 ′ , 2 ′ and 3 ′ ), and a two generator power conversion layer α ′′ (node 1 ′′ and 2 ′′ ). Both layers have exogenous sources and/or sinks, referring to inflow (generation) or outflow (demand), respectively. The flow in both layers is a linear function of donor mass, where in layer α ′ the proportion refers to a constant discharge rate ω(u ′ , v ′ ) = 0.5 , Definition 7). Coefficient k refers to a conversion constant.
Node centrality and flow interpretation. The state of each node in Fig. 3a can be determined from a steady state solution of the random walk evolution process, where based on 46 , k would be considered a weight ω (a form of layer switching probability), resulting in Fig. 3b. For the presented class of problems, however, the lack of approach 46 lies in the fact that each edge of the multilayer network is given uniform treatment, and no distinction is made with respect to inter-layer coupling in the derivation of flow (compare Fig. 3c,d, relating to Fig. 3a,b,  respectively). First of all, the induced inter-layer interaction is a flow conversion (Definition 9) and does not participate in conservation of flow of the source layer, but it does participate in overall conservation of mass. Secondly, the flow conversion acts as exogenous flow with respect to the sink layer, influencing, along with network topology, the respective state probabilities or node potentials of the sink (Fig. 3c). The outcome in the first place www.nature.com/scientificreports/ www.nature.com/scientificreports/ is a different result obtained for network flow (Fig. 3a), which directly affects derivation of the corresponding flow based centrality measure C RW 44 , as where τ st (u) corresponds to flow throughput of node u, normalized by maximum st-flow. The second difference relates to centrality measure interpretation (Fig. 3a,b), given that inter-layer mapping of flow has different effects on different layers, i.e.: it affects the corresponding sink layer (as exogenous flow), however, it does not affect the source layer it "originated" from. The difference is evident in obtained node centrality (impact/importance) rankings, where removal of the highest centrality node in Fig. 3a (node 3 ′ in layer α ′ ) would result in a collapse of the entire network (given removal of the only sink or demand node in layer α ′ , reducing all flows to 0, due to flow balance and coupling constraints). In Fig. 3b, however, the removal of the highest centrality node would result in a flow rebalance only, not reflecting the true coupling nature of the underlying evolution process. The proposed framework extends therefore differently the notion of flow based centrality to multilayer flow type networks, providing therewith a new perspective to multilayer flow based assessments.
To highlight the benefits of multilayer flow derivation and demonstrate properties of flow based centrality further, an additional comparative analysis is presented, with respect to a classic type of centrality i.e., geodesic centrality. As a convenient candidate to the flow based counterpart, classic node betweenness centrality is assessed. Node betweenness centrality 44 is derived as C BC (u) = s� =u� =t∈V σ st (u)/σ st , or the share of times a node u participates in the number of shortest paths σ st between single-source s and single-sink t of a network, where a shortest path is obtained as a path of minimum weighted length, along a sequence of non-repeating nodes and edges. It is easy to see that each node in Fig. 3a,b participates equally in each shortest path of the network, resulting in a uniform node betweenness centrality C BC across nodes. Ultimately, no distinction is made between functional properties of nodes, where again the most critical node is not identified, as opposed to the ranking derived from flow based random walk node betweenness C RW . It is here worth stressing that geodesic centrality, unlike flow based centrality, is concerned with topology only and follows a shortest path intuition. In most networks, however, transition does not unfold along geodesic trajectories and is often dictated by node potentials and exogenous forces, as demonstrated in Fig. 3c. To highlight the difference, the example is expanded with an assessment of the classic node betweenness centrality (as one type of geodesic centrality), to demonstrate disadvantages of disregarding flows in centrality evaluations of flow type problems.
Graph learning-relationship inference. Application of the proposed framework can also be found in the domain of graph learning 51 , specifically in inference of multilayer relationships over multimodal data. Graph learning aims at uncovering structure in data through network inference, where relationships between data points are encoded in the network topology, while data points themselves are represented as node signals. Models can be statistical, physical or graph signal processing based 51 , where for the example described in this subsection a simple implementation of an adaptive gradient based steepest descent method is presented for illustrative purposes, to demonstrate application potential of the proposed framework. To this end, a synthetic dataset is derived based on Eq. (14) and Fig. 4a, where q corresponds to graph node signals perturbed by additive (uniformly distributed) zero-mean noise, while U refers to independent realizations of labelled outcomes, interpretable also as exogenous flow. The model of finding relationships over features q can be formulated as a least mean square (LMS) or Widrow-Hoff adaptive algorithm 52 , where the optimal weigh matrix solving min ω E|U + �q| 2 can be approximated via recursion such that −1 is an initial guess, i an iteration step corresponding to observation i, and µ a sufficiently small constant step size. Feature vector q i and outcome vector U i , are both centred and normalized by observation (ensemble) i for adaptive updates.
It is worth noting that, unlike a supra-Laplacian 2 , which is only an augmented Laplacian matrix over links and layers in flat form, the weight matrix contains conversion ratios of the multimodal relationship as offdiagonal block elements, where such that inter-layer relationships participate in the derivation of the layer structure, as off-diagonal exogenous block entries (Eq. 18), while layers are identified from diagonal block matrices (Eq. 17) (Fig. 4b). The inference algorithm is consequently layer invariant, while identification of mono or multilayer structure resides in interpretation of the obtained weight matrix (e.g., zero column and/or row block matrix sums, Fig. 4).
Ultimately, the proposed framework is capable of picking up a layered (exchange & conversion) as opposed to a flat (exchange only) relationship structure, in relationship inference performed on multimodal time series data. If, for instance, the six nodes in the example corresponded to financial time series, there would be no clear semantic feature distinction, yet a layered structure would exist (as layer 2, the bottom layer, accrues interest/ www.nature.com/scientificreports/ exogenous inflow from layer 1, the top layer). It is worth noting that the relationships within each layer are different to the relationships between the layers, which can erroneously be omitted by a classic flat relationship inference approach. The proposed framework, however, enables capturing such complex network topology and the underlying relationships, without any a priori knowledge about the underlying data sets.
Algorithms and model interpretation. Overall, the LMS algorithm corresponds to a simple implementation of an adaptive gradient based steepest descent method, where signal statistics are replaced by suitable approximations (here instantaneous values) and convergence is dictated by step size. Different approximations may lead to different algorithms and performance 52 , where step size can be iteration dependent and regularization may be employed. The presented example is only illustrative, aimed at demonstrating capabilities of the proposed framework (where the performance of the implemented LMS algorithm could be further improved, though this is outside of the scope of the present manuscript). The key message conveyed is that the proposed framework opens up a new direction for graph learning based applications, by introducing a new type of mathematical object, where other types of relationship inference and graph learning algorithms will be explored separately, as part of future work. It is here worth stressing that a form of graph learning over multilayer networks has been explored in recent work 53 , with reference to Structural Equation Models (SEMs), in form of multilayer SEMs (ml-SEMs). It is, however, worth noting that the proposed HFN framework is a generalization over 53 , where ml-SEMs represent a special case of HFNs, such that r i ,u i =r i ω(u i , r i ) = 1 . Where layers represent temporally lagged snapshots of a monoplex, the proposed framework can be interpreted as a Structural Vector Autoregressive Model (SVARM), as in 53 .
, corresponding to algorithm performance, obtained over 30 experiments R. Note that exogenous block entries of (corresponding to K) relate to conversion of centred data. When solving for e.g., multilayer network flow, this would not pose a limitation, as one could derive the node potentials q based on the obtained weight matrix (i.e., its pseudoinverse) and the given exogenous flow U. Ultimately, this would not distort the flow derivation, as potentials are fixed to within a constant (here average). Derivation of exact conversion coefficients k, however, would require additional constraints and a different approach, which is outside of scope for now.  54 is analysed, and also put into perspective with ml-SEM 53 . The dataset 54 contains cross-industry use of commodities and GDP gross output by industry, over 71 industries, 21 industry groups, and a time range 1997-2019, expressed in millions of dollars. To simplify further exposition, a simple example (Fig. 5) is used as a guide throughout this subsection, in form of a proxy equivalent to 54 . The studied dataset (Fig. 5b) consists of domestic supply (rows) and use (columns) of commodities by industry (nodes), in form of a flow matrix F, by corresponding year α . Supply relates to node potentials in this case, where rate of transfer (arc weight) ω is derived from cross-industry flow, as a share of compartmental (supply) outflow (Fig. 5a). Compartmental weight matrix is derived from the rate of transfer ω (Fig. 5c), where r i ,u i =r i ω(u i , r i ) = 1 for certain nodes (see diagonal), as some funds get reinvested over time or further exogenous investments are introduced (see F diagonal). Inter-layer coupling corresponds to transformation of net output by industry, where conversion k refers to a proportion defined with respect to e.g., rate of return. Multilayer structure comes hence from economic activity (Fig. 5a,b,d), rather than industrial clustering suggested in 53 . Each layer can therewith be interpreted as a snapshot in time, or equivalently as a Petri net state (Fig. 5a,e,f). Note that 53 relates total use, as opposed to total supply, to node potentials. (c) Compartmental weight matrix , derived from rate of transfer arc weights ω . Note that r i ,u i =r i ω(u i , r i ) = 1 (see weight matrix diagonal) for some nodes (i.e., nodes {1 ′ , 3 ′ , 1 ′′ , 2 ′′ } ), as some funds get reinvested over time or further external investments are introduced (see flow matrix F diagonal www.nature.com/scientificreports/ To derive the exact conversion coefficient k, one requires information about the capital carried over from a previous time instance and the additional investment at a given time (compare nodes 1 ′′ and 2 ′′ , Fig. 5a,b; both nodes have nonzero F diagonal entries). From there, k is derived as a proportion between exogenous outflow from the previous time instance and accrued capital in the next one (Fig. 5a,d). The value added may correspond to e.g., accrued interest, which would not be captured by a traditional approach, e.g. 46 . The proposed framework, in contrast, enables a multilayer flow interpretation of economic activity in 54 , uncovering some interesting properties.
Economic robustness and industry ranking. The derived network equivalent finds first of all interpretation in the domain of Stochastic Petri nets (SPN), where an SPN with given transition rates induces a Markov chain on its reachable set 38 . One can hence derive a corresponding steady state probability distribution, where elements refer to probability of being in a state Q i (here node potentials q on a layer α i ).
One can, furthermore, derive a flow based random walk node betweenness centrality C RW for each industry, directly from cross-industry and cross-layer flow (Figs. 5g, 6b). It is here interesting to note that all industries exhibit a consistent centrality over a 23 year time span (Fig. 6b,c), regardless of domestic output fluctuation (Fig. 6a). This means that the relative importance of each industry i.e., the relative monetary throughput with respect to other industries and even with respect to a different monetary volume and economic environment (i.e., the 2008 economic crisis), remains absolutely the same, which suggests some form of relative network (i.e., economic activity) robustness. The industry ranking based on random walk node betweenness centrality C RW is also more stratified (Fig. 6b), as opposed to domestic output ranking alone (Fig. 6a), where it should be noted that the centrality ranking takes both, domestic output and supply of commodities by industry (Fig. 5a,g), into consideration. The proposed framework is compared with the real situation, and introduces a multilayer network flow interpretation of U.S. GDP by industry. The resulting node centralities, as a network theory tool, provide insights on the relative importance of each industrial sector, revealing consistency among relative positions of industries over a time period of more than two decades, and a finer more stratified ranking as compared to domestic output ranking alone. They identify therein five key industrial sectors (Fig. 6b), their relative positions in monetary throughput, as well as relative robustness under different economic environments (consistent centrality, Fig. 6b,c, regardless of domestic output fluctuation, Fig. 6a), where exogenous flows are interpretable as, and act in form of, control inputs. The proposed framework effectively captures the role of nodes (industries) as not necessarily generators of outputs (products) alone (as measured by the classic gross domestic output ranking, Fig. 6a), but also as suppliers of commodities required by other industries (as captured by the proposed framework and the implemented flow centrality ranking, Fig. 6b,c), providing entirely new insights on an economy, its industries, their importance, and relative robustness.

Discussion and outlook
The presented examples demonstrate some of the benefits of multimodal flow derivation in complex network analysis. Inter-layer flow conversion takes place across network links, which act as exogenous control flows with respect to the sink layer. The multimodal flows are implicitly either coupled (across layers) or uncoupled (in intra-layer flow derivation, specifically with regard to cycle space conditions). This property lends some new www.nature.com/scientificreports/ perspective to network interdependency assessment. Furthermore, establishing correspondence between two apparently distinct conceptual domains brings new insights to both formalisms, as the tools and models of one can be used to develop an understand tools and models of the other, and vice versa. As previously mentioned, the heterogeneous flow network enables derivation of a layered relationship structure (corresponding to connectivity and conversion of node data), as opposed to a classic flat relationship structure (corresponding to connectivity of node data only). This enables a physical interpretation of models with complex interactions between different semantic domains (e.g., chemical processes, energy networks, logistics, finance, or any other form of conversion process relying on the laws of conservation), in the form of multilayer network flows. Disregarding such context may lead to erroneous interpretations, as demonstrated in the presented examples. The Petri net, on the other hand, enables flattening of the layered relationship structure (not to be confused with the latter, classic flat relationship representation). The construction simplifies computation, as it reduces the complex multilayer network to a form of bipartite graph (of resource and transition nodes), at the expense of losing semantics. As the proposed framework, however, establishes reversible correspondence between the two mathematical objects, interpretation is always preserved and outputs can uniquely be transformed from one mathematical object to the other. In a hypothetical 'ablation' analysis sense, Petri nets are needed to facilitate computation, while heterogeneous multilayer flow networks provide interpretation and access to well-established network theory tools. Removing flow from the formulation would leave a very limited scope for quantitative assessments (i.e., disregarding node potentials and exogenous forces), which would inhibit reasoning, as it would remove the building blocks of the underlying relationship structure (e.g.: feature vectors and labelled data, respectively; node probabilities with respect to random walks on multilayer networks; mutual relationships between nodes and exogenous control flows; etc.). It is also worth adding that the proposed framework produces outputs that are analytically obtained and are fully interpretable.
In prospect, one of the anticipated goals is to give further attention to the extraction of the exact conversion coefficient k in graph learning based relationship inference (learning based function approximation) problems. From a synthesis point of view it is clear that, at least in some cases, k is uniquely determined. Further interesting insights could be obtained through possible integration of nonlinear relationships and dynamic topology considerations within the proposed framework, drawing on rich literature from current state-of-the-art.

Conclusion
In this paper the formal notion of heterogeneous network flow is proposed, as a multilayer flow function aligned with the theory of network flow. A dynamic equivalence with the framework of Petri nets is established, as the baseline model of concurrent event systems. The Petri net flow relations are here extended, to possibly incorporate both fundamental equations of balance, namely: flow balance, which is integral to the Petri net model, and node potential balance (cycle space condition), which may arise in relation to specific application domains. Overall, the key property of the proposed framework lies in the ability to derive multimodal flows across different semantic domains, satisfying conditions of cross-layer conservation and coupling. The proposed framework provides therewith an extension to the recently introduced field of multilayer networks.
Some covered applications include multilayer Laplacian flow and multilayer flow centrality, as well as graph learning based inference of multilayer relationships over multimodal data. On synthetic data the proposed framework demonstrates benefits of multimodal flow derivation in critical component identification. It also displays applicability in relationship inference (learning based function approximation) performed on multimodal time series. On real-world data the proposed framework provides, among others, multimodal flow interpretation of U.S. economic activity, uncovering underlying empirical steady state probability distribution, as well as inherent network (economic) robustness.