Supp_MultiControl_2021_02

# Structural underpinnings of control in multiplex networks

Pragya Srivastava  Department of Bioengineering, University of Pennsylvania,
Pennsylvania, PA 19104, USA  Peter J. Mucha  Department of Mathematics,
University of North Carolina, Chapel Hill, NC 27599, USA  Department of
Applied Physical Sciences, University of North Carolina, Chapel Hill, NC
27599, USA  Emily Falk  Annenberg School for Communication, University of
Pennsylvania, Philadelphia, PA 19104, USA  Fabio Pasqualetti  Department of
Mechanical Engineering, University of California, Riverside, CA 92521, USA
Danielle S. Bassett  Department of Bioengineering, University of Pennsylvania,
Pennsylvania, PA 19104, USA  Department of Physics & Astronomy, University of
Pennsylvania, Pennsylvania, PA 19104, USA  Department of Electrical & Systems
Engineering, Pennsylvania, PA 19104, USA  Department of Neurology,
Pennsylvania, PA 19104, USA  Department of Psychiatry, Pennsylvania, PA 19104,
USA  Santa Fe Institute, Santa Fe, NM 87501, USA

###### Abstract

To design control strategies that predictably manipulate a system’s behaviour,
it is first necessary to understand how the system’s structure relates to its
response. Many complex systems can be represented as multilayer networks whose
response and control can be studied in the framework of network control
theory. It remains unknown how a system’s layered architecture dictates its
control properties, particularly when control signals can only access the
system through a single input layer. Here, we use the framework of linear
control theory to probe the control properties of a duplex network with
directed interlayer links. We determine the manner in which the structural
properties of layers and relative interlayer arrangement together dictate the
system’s response. For this purpose, we calculate the exact expression of
optimal control energy in terms of layer spectra and the relative alignment
between the eigenmodes of the input layer and the deeper target layer. For a
range of numerically constructed duplex networks, we then calculate the
control properties of the two layers as a function of target-layer densities
and layer topology. The alignment of layer eigenmodes emerges as an important
parameter that sets the cost and routing of the optimal energy for control. We
understand these results in a simplified limit of a single-mode approximation,
and we build metrics to characterize the routing of optimal control energy
through the eigenmodes of each layer. Our analytical and numerical results
together provide insights into the relationships between structure and control
in duplex networks. More generally, they serve as a platform for future work
designing optimal _interlayer_ control strategies.

##  I Introduction

The field of network science provides a powerful framework for studying the
structure and organization of many complex systems including socioeconomic
networks, traffic and communication networks, physical networks, and
biological networks Gosak_bioNetsreview_2018 ; bassett2018on ;
LiaGranNetRev_2018 ; jackson_social_evol ; jackson_socio_book . Markedly
useful and still an active area of research, the network representation of a
complex system for a given research problem most often takes into account
pairwise interactions of only one type LeoAnn_complexreview . However, in most
complex systems, different types of interactions among system components can
coexist and influence one another bianconimultilayer ; kivela2014multilayer .
For example, in a transportation network people can move along bus, subway, or
train routes. In a socioeconomic network, the same set of individuals or
institutions can interact socially as well as via financial transactions.
Financial networks themselves can have multiple layers of interactions amongst
their components bookstaber_multifinance ; Wojcik_2018 . In a social network,
the same set of users can interact in-person or via multiple social platforms
Dickison_multisocial ; atkisson_arxiv_multi . Finally, in neuronal networks
information can propagate along electrical and chemical synapses. The
resultant structural picture of the system is thus that of a multilayer
network where each component (node) can engage in multiple types of
interactions with another component (node) across many temporal and spatial
scales Buldu_2017 ; Buldu_comment_2018 ; Muldoon2018comment ;
PilosofMultilayerEcology_2017 ; domenico2013mathematical ;
bentley2016multilayer ; kidd2014unifying ; Schmidt2018:MultiScale ;
bassett2011understanding ; braun2018maps .

Multilayered networks in true biological, physical, social, and technological
realms do not exist in isolation, but instead impact —and are impacted by
—other systems around them through a variety of forces, influences, and causal
interactions Magnusson_1996 ; Bronfenbrenner_1979 . Yet, precisely how
multilayer networks respond to perturbative signals, either beneficent or
maleficent, remains far from understood. A promising route to progress is the
utilization of linear systems theory, which provides a standard framework to
understand the response of a complex system to both extrinsic and intrinsic
perturbations kailath1980linear ; kubo_toda_NEQ_linear . The theory of linear
control then utilizes the response of the system to determine its control
properties and to inform the development of appropriate control strategies
kailath1980linear . An important part of characterizing system response is
determining how a system’s control properties depend upon its structure. While
in monolayer networks the link between the structure and control properties
has been investigated JasonKim_2018 ; Tang_2018 , this task becomes
particularly challenging for multiplex networks.

Several recent studies have investigated the controllability of multiplex
networks from various perspectives and have provided important initial
insights. Using the minimum number of driver nodes $n_{D}$, that ensures the
full controllability of _directed_ duplex networks, Posfai _et al_. calculated
the contribution of a given layer to the full controllability of a duplex
network as a function of the relative time scales of activity update in the
two layers posfai2016congtrollability . In a complementary study, Menichetti
_et al_. investigated the role of degree correlations across layers on the
minimum number of driver nodes necessary for control, particularly in settings
where only a fixed set of nodes are control nodes in all layers
menichetti2016control . Moving beyond directed multiplex networks, Yuan _et
al_. obtained an exact expression for $n_{D}$ in terms of network size and the
rank of the adjacency matrix associated with the duplex network
Yuan_exact_control_multi . Finally, Wang _et al_. identified a trade-off
between controllability and control energy for different patterns of
interlayer connections in a duplex network Wang_multilayer_control_2017 .
Although these studies have explored different aspects of the structure-
control relationship, a generic framework that allows us to extract the
dependence of multiplex control on layered architecture is currently lacking.

Here, we investigate the relationship between structure and control in duplex
networks, where we seek to drive the state of a deep, inaccessible target
layer by injecting control signals into a superficial, accessible input layer
(Fig. 1(A)). Such control requirements could arise in real multilayered
systems where external signals can manipulate a specific type of interaction
due to constraints or convenience. For instance, in the context of social
networks an objective could be to modify the behaviour of individuals in their
professional environment (professional layer) by providing resources of
support in their informal social environment. In multilayered socioeconomic
networks, reduction of vulnerability to financial crises in a network of
financial institutions can be achieved by designing optimal intervention
strategies by the government and public. In multilayered brain networks, an
important control objective could be to manipulate inter-regional connectivity
(i.e. the structural connectivity matrix) by targeting the activity of
specific brain-regions.

Having specified such a control configuration, we ask how the structural
properties of individual layers (eigenspectra) and their relative arrangement
(alignment between layer eigenmodes) together impact control. Using the
framework of linear control theory in a layer-diagonalized representation, we
determine the optimal energy required to drive the system to a given final
state, and show analytically how it depends upon these two structural
parameters. Using numerical calculations, we then estimate the optimal control
energy for duplex networks whose layers are Erdős-Rényi (ER), Watts-Strogatz
(WS), Barabási-Albert (BA), or Random-Geometric (RG) graphs. We calculate the
average and maximum control energy of the target layer as a function of the
layer density and topology. We show that the alignment between the eigenmodes
of the two layers determines the cost of controlling specific eigenmodes of
the target layer. Moreover, the angle of alignment determines the number of
eigenmodes that are excited by the optimal control energy. More specifically,
we demonstrate the routing of optimal control energy through the eigenmodes of
the input layer in a manner that depends systematically on their alignment
with the target eigenmode. The approach and results presented in this work can
be utilized to examine the interlayer relationships in complex multilayer
systems, identify and predict the easy versus difficult control directions,
and analyse the physical paths along which control is exercised.

The remainder of this paper is structured as follows. In Section II, we
formulate the problem of controlling a duplex with linear time-invariant (LTI)
dynamics and outline our approach. In Section III, we present the numerical
results on the structural dependence of the average control energy that
emphasize the role of target-layer topology, interlayer connections, and
alignment. Then, in Section IV we deepen our understanding of these results in
the simplified limit of a one-mode approximation. In Section V, we use the
insights from our results to build metrics to identify the paths along which
optimal control signal is distributed. Finally, in Section VI, we conclude our
investigation and discuss future directions.

Figure 1: Control-theoretic framework for duplex networks. (A)-(B) Control
signals that are injected into the superficial, accessible input layer of the
duplex network influence the state of the input layer, which in turn
influences the state of the deep, inaccessible target layer via directed
interlayer connections. (B) Stylized schematics of the energy landscapes of
the input and the target layers on which the state trajectories are navigated
by control signals. (C) The subdiagonal block of the supradjacency matrix for
the duplex networks in this work has non-zero diagonal entries corresponding
to the directed interlayer connections between replica nodes. The intra-layer
networks have undirected connections resulting in a symmetric adjacency
matrix. (D) Optimal energies for both layers are calculated corresponding to
the navigation of their state along their eigendirections (denoted by
$p_{\alpha}$ for the input layer and $q_{\beta}$ for the target layer) by a
unit amount ($w_{\alpha F}=1$ for the input layer, and $v_{\beta F}=1$ for the
target layer) from an initial state at the origin. In the eigenspace of each
layer, final states ($X_{1F}$ for the input layer and $X_{2F}$ for the target
layer) fall on the surface of the unit sphere, while the corresponding
energies ($E_{1}^{\alpha}$ for the input layer and $E_{2}^{\beta}$ for the
target layer) constitute the axes of an ellipsoid Wicks1988energyapproach .

##  II Formulation

We consider a duplex network in which control signals are injected into the
first layer, also called the ‘input’ layer (Fig. 1(A)). Practically, we think
of the input layer as a superficial layer in a real system which is accessible
to external perturbations. Since different layers of a multiplex network
represent different types of interactions between system-components, this
control setup mimics the case when only one type of interaction can be
directly affected by the external perturbations. Any changes in the state of
the first layer flow to the second ‘target’ layer via directed interlayer
links (Fig. 1(A)-(B)). For simplicity, we assume that the intralayer
connectivity matrices are undirected. The supra-adjacency matrix $M$ for such
a duplex network has the following form kivela2014multilayer :

$\displaystyle M=\left(\begin{array}[]{cc}A^{(1)}&0_{N}\\\
K&A^{(2)}\end{array}\right),$ (3)

where $A^{(1)}\in\mathbb{R}^{N\times N}$, $A^{(2)}\in\mathbb{R}^{N\times N}$,
and $K\in\mathbb{R}^{N\times N}$ (Fig. 1(C)). Here, $A^{(1)}$ and $A^{(2)}$
denote the _intralayer_ adjacency matrices of the two layers, and the matrix
$K$ encodes the interlayer connections. In a duplex topology, each layer
consists of the same set of nodes. Moreover, the interayer connections exist
only between the same nodes; that is, a given node in a layer connects to its
‘replica’ node in a different layer. Thus, the interlayer connectivity matrix
$K$ assumes the form $K=\kappa I_{N}$, where $\kappa$ denotes the weight of
interlayer connections, and $I_{N}\in\mathbb{R}^{N\times N}$ is the identity
matrix in $N-$ dimensions. Finally, $0_{N}$ denotes an $N-$ dimensional square
matrix with zero entries.

The control signal $u(t)$ is injected into the first layer of this duplex via
an input matrix $B$. We model the evolution of the network’s state by a linear
time-invariant dynamics in continuous time:

$\displaystyle\frac{dX}{dt}=MX(t)+Bu(t),$ (4)

where $X(t)\in\mathbb{R}^{2N\times 1}$ is a vector that denotes the state of
nodes in all network layers at time $t$. We assume that the control signal is
injected into all of the nodes of the input layer, and thus $B$ is $B\in
R^{2N\times N}$ and $B=\left(I_{N}\,\,0_{N}\right)^{T}$. Clearly,
$u(t)\in\mathbb{R}^{N\times 1}$. Finally, the cost of control is measured by
the control energy $E$, which is defined as:

$E=\int_{0}^{t_{F}}u^{T}(t)u(t)dt,$ (5)

where $t_{F}$ is the time horizon of the control, and the superscript $T$
denotes the matrix transpose.

In order to explicitly link the structure of a duplex network to the cost of
control, we first express state equations in a layer-diagonalized
representation. Specifically, we first separate the state vector $X$ in
$X^{(1)}\in\mathbb{R}^{N\times 1}$ and $X^{(2)}\in\mathbb{R}^{N\times 1}$,
where the former represents the state of the nodes in the input layer and the
latter represents the state of the nodes in the target layer. Then, the LTI
dynamics in continuous time (Eq. (4)) can be written as:

$\displaystyle\dot{X}^{(1)}$ $\displaystyle=A^{(1)}X^{(1)}+u$ (6a)
$\displaystyle\dot{X}^{(2)}$ $\displaystyle=KA^{(1)}+A^{(2)}X^{(2)}.$ (6b)

We diagonalize the real-symmetric matrices $A^{(1)}$ and $A^{(2)}$ as
$A^{(1)}=PD_{1}P^{T}$ and $A^{(2)}=QD_{2}Q^{T}$, where
$Diag(D_{1})=\\{\xi_{\alpha}|\alpha\in\\{1,2,\ldots,N\\}\\}$ and
$Diag(D_{2})=\\{\mu_{\alpha}|\alpha\in\\{1,2,\ldots,N\\}\\}$. In addition, the
matrices $P$ and $Q$ have the form
$P=\left[\mathbf{p}_{1},\mathbf{p_{2}},\ldots,\mathbf{p}_{\alpha},\ldots,\mathbf{p}_{N}\right]$,
and
$Q=\left[\mathbf{q}_{1},\mathbf{q_{2}},\ldots,\mathbf{q}_{\alpha},\ldots,\mathbf{q}_{N}\right]$,
such that $\mathbf{p}_{\alpha}$ and $\mathbf{q}_{\alpha}$ are the
$\alpha^{\rm{th}}$ eigenvectors of $P$ and $Q$, respectively. Using the above
decomposition, we next write,
$X^{(1)}=\underset{\alpha=1}{\Sigma}w_{\alpha}\mathbf{p}_{\alpha}$ and
$X^{(2)}=\underset{\alpha=1}{\Sigma}v_{\alpha}\mathbf{q}_{\alpha}$. Rewriting
Eqs. (6a) and (6b), we obtain:

$\displaystyle\dot{w}_{\alpha}$
$\displaystyle=\xi_{\alpha}w_{\alpha}+U_{\alpha}$ (7a)
$\displaystyle\dot{v}_{\alpha}$
$\displaystyle=\mu_{\alpha}v_{\alpha}+\kappa\underset{j}{\Sigma}w_{\beta}\cos\theta_{\alpha\beta},$
(7b)

where we have used the fact that $u(t)=\Sigma U_{\alpha}\mathbf{p}_{\alpha}$,
and that $K=\underset{\beta}{\Sigma}\mathbf{q}_{\beta}\mathbf{k}_{\beta}$. For
the duplex network considered here, $K=\kappa I_{N}$ leading to
$\mathbf{k}_{\beta}=\mathbf{q}_{\beta}^{T}$. This, in turn, leads to the
identification of
$\theta_{\alpha\beta}=\arccos(\mathbf{q}_{\beta}^{T}\mathbf{p}_{\beta})$ as
the angle between $\mathbf{q}_{\alpha}$ and $\mathbf{p}_{\beta}$. These
decompositions correspond to the projection of $u(t)$ and the matrix $K$ along
the eigenvectors of the input layer and the target layer, respectively. The
resulting Eqs. (7a) and (7b) have an explicit dependence on the eigenvalues of
each layer as well as on the angles $\theta_{\alpha\beta}$ that characterize
the alignment between the eigenmodes of the two layers. The state equations
Eqs. (7a) and (7b) are obtained using the transformation which diagonalises
the adjacency matrix of each layer, which we refer to as the ‘layer-
diagonalised’ representation of the state equations.

In a controllable system, there could be multiple control inputs that navigate
the system between a pair of states. The task then becomes to find the minimum
control energy that is required to control the system from an initial state to
a desired final state. Mathematically, this minimization problem is expressed
as:

$\displaystyle\underset{U(t)}{{\rm{minimize}}}\quad E$
$\displaystyle=\int_{0}^{t_{F}}U^{T}(t)U(t)dt,$ $\displaystyle\text{subject to
}\dot{w}_{\alpha}$
$\displaystyle=\xi_{\alpha}w_{\alpha}+\mathbf{b}_{\alpha}\cdot U$
$\displaystyle\dot{v}_{\alpha}$
$\displaystyle=\mu_{\alpha}v_{\alpha}+\kappa\underset{j}{\Sigma}w_{\beta}\cos\theta_{\alpha\beta}.$
(8)

This problem can be solved analytically using the standard framework of
optimal control theory lewis_optimal_control_2012 . We apply this framework to
the problem of Eq. (8) and provide the exact solutions for the state variables
$\\{\\{w_{\alpha}\\},\\{v_{\alpha}\\}\\}$ and the optimal control input
$u_{\alpha}=\mathbf{b}_{\alpha}\cdot U$ in closed form (see the supplementary
section  S1).

We note that our solutions to the minimization problem display a complex
dependence on the two structural parameters of interest: the eigenvalues of
each layer and the alignment of the eigenmodes of the two layers. This
dependence renders our solutions rather opaque, and it remains difficult to
extract simple insights into the structural determinants of control. To make
progress, we therefore now turn to a numerical computation of the optimal
control energy for duplex networks with layers constructed from canonical
graph models. After this numerical assessment, we work in a single-mode limit
to more deeply understand our results.

##  III Structure-control relationships in duplex networks: A numerical
investigation

Figure 2: Control properties of simulated duplex networks. (A) Schematic of
the process by which we construct duplex networks. Each layer belongs to one
of the following graph families; Erdős-Rényi (ER), Watts-Strogatz (WS),
Barabási-Albert (BA), and Random Geometric (RG). An example network from each
of these graph families is shown on the left. (B) The plots of average control
energy and (C) the maximum control energy of the target layer as a function of
the target-layer density for all graph families studied. Each column
corresponds to the topology of the target layer. Throughout this paper, we use
differently-colored curves to indicate the topology of the input layer, and we
use different markers to indicate the topology of the target layer. The
average control energy is lowest for the duplex network with a WS network as
the input layer and an RG network as the target layer. The loops in the plots
corresponding to a BA network in the target layer are indicative of a non-
monotonic dependence of network density on the generating parameter (see
supplementary section and Fig.  S3) (D)-(E). To further understand our results
in graph model families, we build a set of hierarchically modular (HM)
networks as a function of two generating parameters: $\gamma$, which
determines the degree heterogeneity, and $f$, which determines the modularity.
Here we show heatmaps of degree heterogeneity and average clustering for the
HM networks. (F)-(G) The heatmaps of average and maximum optimal control
energies of the target layer in duplex networks constructed from HM networks.

In this section, we present our numerical results on the control properties of
a duplex network as a function of the properties of its layers and interlayer
alignment. In the previous section, we determined the dependence of control
energy on: (i) the two sets of eigenvalues corresponding to the two layers of
the duplex, and (ii) the angles characterizing the alignment between the
eigenvectors of the two layers. Network properties that directly influence the
eigenvalues are the density of the network and the topology of the network. To
simulate the variations in topology, we construct duplex networks where each
layer belongs to one of the following four graph families: Erdős-Rényi (ER),
Watts-Strogatz (WS), Barabási-Albert (BA), and Random-Geometric (RG). Networks
in these four families differ in their topology, exhibiting distinct degree
heterogeneity and clustering coefficient Albert-Barabasi_1999 ; Albert-
Barabasi_2002_RMP ; Penrose2004_RGG ; watts_strogatz_1998  (see the
supplementary section and Fig. S2); moreover, the RG family is spatially
embedded while the other 3 are not. In addition, several real systems can be
represented by networks from these four families, making their study germane
to a generic understanding biological, physical, and social systems
Albert_SFBiol_2005 ; Keller_2005_SFessay ; Telesford_2011_WSubiquity ;
Penrose2004_RGG . To simulate variations in network density, we changed the
density of the target layer while holding the density of the input layer fixed
(see the supplementary section and Fig. S3).

To explicitly relate the control properties of the duplex networks to their
structure, we calculate the optimal control energy required to drive the
system from the origin to a set of pre-specified final states. We chose the
final states to be those that have unit norm and that are aligned with the
eigenvectors of the target layer (Fig. 1(D)). Control energies thus calculated
define the intersections of the control energy surface with the
eigendirections of each layer. Intuitively, these control energies
characterize the difficulty of navigating the state of a layer along a given
eigendirection. We use the quantity $E_{\ell}^{\alpha}$ to denote the control
energy of navigating the state of the $\ell^{\rm{th}}$ layer along the
$\alpha^{\rm{th}}$ eigendirection by a unit amount. Then, a measure of the
average controllability of a given layer can be defined as

$\displaystyle\bar{E}_{\ell}=\underset{\alpha}{\Sigma}E_{\ell}^{\alpha},\quad\ell\in\\{1,2\\}.$
(9)

Similarly, the maximum control energy corresponds to the eigendirection along
which it is most difficult to control the system. To extract this direction,
we define:

$\displaystyle
E_{\ell}^{\rm{max}}=\rm{max}\\{E_{\ell}^{\alpha}:\alpha\in\\{1,2,\ldots,N\\},\ell\in\\{1,2\\}\\}.$
(10)

For a given duplex network, the $E_{\ell}^{\alpha}$ quantify the difficulty of
navigating the state of the $l^{\rm{th}}$ layer along its $\alpha^{\rm{th}}$
eigenvector, thus giving access to the complete spectrum of control
properties. We calculate $E_{\ell}^{\alpha}$ for all duplex networks
constructed from the same four graph families studied earlier (see the
supplementary Section and Fig. S4). Since we primarily focus our attention on
controlling a layer by applying the control in a different layer, in the
following section we discuss the control energies only for the target layer.
The following subsections report our results first for the eigenvalues of each
layer and then for the alignment between the eigenmodes of the two layers.

###  III.1 Density- and topology-dependence of the target-layer control

To calculate the control energies of the target layer, we create sixteen
different duplex networks by constructing each layer from the 4 graph families
(Fig. 2(A)). We normalize the supra-adjacency matrix of each duplex by the
largest eigenvalue $\xi_{\rm{max}}$ of the first layer, which enables us to
consistently compare the results obtained for the different duplex networks.
Since each eigenvalue has the dimensions of a rate, intuitively this
normalization sets the unit of time as the inverse of the largest eigenvalue
of the input layer. We set the strength of the interlayer connections,
$\kappa$, to be the same as the strength of the intralayer connections in the
_unnormalized_ network. Thus in normalized duplex networks, the strength of
the interlayer connections is set to $\kappa/\xi_{\rm{max}}$. We return to the
dependence of control energy on $\kappa$ later in this paper (also see the
supplementary section and Fig. S5). Finally, we calculate the quantities
$\bar{E}_{2}$ and $E^{\rm{max}}_{2}$ as a function of the density of the
target layer, for 40 realizations of each duplex network.

We find that the average control energy $\bar{E}_{2}$ depends upon the density
as well as the topology of the input and target layers (Fig. 2(B)). For all
duplex networks, the dependence of $\bar{E}_{2}$ on the target layer density
is non-monotonic; consistently, the intermediate densities are the hardest to
control. Further, $\bar{E}_{2}$ is lowest when the topology of the target
layer is that of an RG graph, for all the duplex networks and for all density
values. For a fixed topology of the target layer, $\bar{E}_{2}$ is lowest when
the topology of the input layer is that of a WS network. Thus, a duplex
network with a WS network in the input layer and an RG network in the target
layer has the lowest average control energy.

Next, we find that the maximum control energy $E^{\rm{max}}_{2}$ also depends
upon the density of the target layer (Fig. 2(C)). For all duplex networks, the
maximum control energy is highest in the intermediate range of densities. The
topology of the input layer has a weak effect on $E^{\rm{max}}_{2}$, when the
target layer has the topology of an ER or an RG network. However if the target
layer is a WS network or a BA network – and if the input layer is a WS
network–, then the $E^{\rm{max}}_{2}$ is quite high. This effect is most
striking when the target layer has the topology of a BA network. One possible
reason for this behavior is that BA networks are characterized by a power-law
degree distribution, implying the existence of many nodes that have either
much higher or much lower degree than expected in an ER graph. Indeed, we show
a similarity between the degree heterogeneity of the target layer and
$E_{2}^{\rm{max}}$ in the next subsection. Intuitively, a larger number of
sparsely connected nodes could contribute to the eigendirections that are much
harder to control in BA networks.

To evaluate the robustness of our findings, we performed several additional
sensitivity analyses. First, we note that in Fig. 2(B-C) the average and
maximum control energy are plotted for a fixed density of the input layer. We
verify that our results remain unchanged for different densities of the input
layer (see the supplementary section and Fig.  S6). Additionally, similar
plots for $\bar{E}_{1}$ and $E^{\rm{max}}$ can be made (see the supplementary
section and Fig.  S7). Note that, for the _unnormalized_ interlayer connection
strength $\kappa=1$, we find that the optimal control energies to control the
target layer are much larger than that of the input layer. Varying the
strength of $\kappa$ changes the relative values of the optimal control energy
for input and target layers: a large value of $\kappa$ leads to a faster
propagation of control signals between the two layers, thus lowering the
control energies (see the supplementary section and Fig.  S5). Note that we
explain this dependence of control energy on $\kappa$ more formally in the
single-mode limit in Section IV of this main text.

In sum, unlike the input layer in which control signals are injected directly,
the effect of control signals on the target layer is mediated by the input
layer and the interlayer connections. This makes the strength of interlayer
connections and the arrangement of the target layer relative to the input
layer important parameters in determining the cost of controlling the target
layer. In this section, we characterized global metrics of controllability,
namely, the average and the maximum control energy of the target layer, as a
function of the density of the target layer and the topologies of the two
layers. However, the trends we observe (Fig. 2(B-C)) are due to the compound
effects of network density and topology. In the next section, we separate the
role of network topology from that of the network density. In particular, we
focus on the role of heterogeneity in the target layer on its control
properties. To do this, we construct duplex networks whose target layers
continuously span a range of degree heterogeneity and clustering coefficient.
This approach allows us to explicitly link the two topological properties of
the target layer to the cost of their control.

###  III.2 Control properties of a heterogeneous target layer

Many real-world networks exhibit topological properties that differentiate
them from Erdős-Rényi networks. Examples include heterogeneity in node
degrees, high clustering coefficient, and small geodesic distances. Out of the
four graph families considered in the previous sections, BA networks have the
highest degree heterogeneity while RG networks have the highest average
clustering coefficient for a fixed network density (see the supplementary
section and Fig. S2). To gain deeper insight regarding the link between
control and topology, we now calculate the control properties of the duplex
networks whose target layers have a continuous variation in their degree
heterogeneity and in their average clustering coefficient. We fix the topology
of the input layer to be ER, while the target layers are constructed using a
numerical scheme that generates hierarchically modular (HM) networks
Chris_Nphys2020 . The network generation method has two parameters: (i) the
parameter $\gamma$ that shapes the degree distribution by the mechanism of
preferential attachment, and (ii) the parameter $f$ that specifies the
fraction of edges in a module. The resulting networks exhibit a continuous
range of degree heterogeneity and clustering coefficient (Fig. 2(D-E)).

Using these continuously varying hierarchical architectures, we calculate the
optimal control energies for the target layers of duplex networks in a manner
similar to that outlined in the previous section. In particular, we
investigate average control energy $\bar{E}_{2}$ and the maximum control
energy $E_{2}^{\rm{max}}$ as a function of the two generative parameters,
$\gamma$ and $f$ (Figs. 2(F-D)). We note a similarity between the plot of
maximum control energy and the plot of degree heterogeneity (Fig. 2(D) and
(G)), suggesting that the former is partially explained by the latter. The
similarity between the plot of average control energy and that of degree
heterogeneity is less strong; while along the $f$\- axis, we note a positive
correlation between the two variables, along the $\gamma$\- axis this trend is
inverted. We also note a negative correlation between the heatmaps of average
clustering coefficient and average control energy (Fig. 2(E) and (F)). This
observation is consistent with the results we report in the previous section,
where the average control energy of the target layer was lowest for an RG
network, a topology that displays the highest average clustering among the
four graph families studied (see the supplementary section and Fig.  S2).
Collectively, these observations serve to support the observed dependence
between control and topology in canonical graph families in a more continuous
assessment of topological variation.

Figure 3: The role of alignment in determining the control cost. (A) The
scheme of rotation. The initial eigenspaces for each layer are shown on the
left. On the right, we show how the eigenspace of the target layer is rotated
relative to the input layer by defining the rotation matrix in the plane of
eigenmodes corresponding to the two most positive eigenvalues of the input
layer. This procedure changes the alignment of eigenmodes of the target layer
only in the rotation plane, while maintaining their original alignment with
the remaining (N-2)-dimensional subspace. (B) Control energy and (C) the
number of excited modes to control the layer state along its dominant
eigenmode as a function of alignment between the dominant eigenmodes of the
two layers. The energy required to control the target layer along its dominant
eigenmode is lowest when the two layers are perfectly aligned. Simultaneously,
the number of eigenmodes that carry the optimal control energy is lowest when
the two layers are perfectly aligned. Colors identify the topology of the
input layer and markers identify the topology of the target layer. Each curve
is an average of 100 simulations. (D-F) The fraction of control energy in the
first (D) and the second (E) eigenmode of the input layer as a function of
their alignment with the target eigenmode. (F) Total fraction of energy in the
first two eigenmodes as a function of the alignment between the dominant
eigenmodes of the two layers of duplex networks.

###  III.3 The interplay between topology and alignment in determining control

From the analysis presented in the previous sections, we see that the
alignment $\cos\theta_{\alpha\beta}$, henceforth denoted as $C_{\alpha\beta}$,
emerges as an important factor that encodes the relative arrangement of the
eigenmodes of the two layers. We now investigate how this parameter affects
the cost of controlling a specific eigenmode in the target layer. For a clear
demonstration of the role of $C_{\alpha\beta}$, we set the control task to be
the navigation of the target layer to a final state that has a unit norm and
that is aligned to the eigenvector corresponding to the largest eigenvalue
(the dominant eigenmode). In order to explicitly control the alignment of the
dominant eigenmode relative to the eigenmodes of the input layer, we construct
a rotation matrix $\mathcal{R}(\phi)$ in the plane spanned by the eigenmodes
corresponding to the two most positive eigenvalues of the input layer (Fig.
3(A)) (see the supplementary section and Fig.  S8) ChrisDoran_GeomAlgeb .
Then, the transformation
$A_{2}^{\prime}=\mathcal{R}(\phi)^{T}A_{2}\mathcal{R}(\phi)$, changes the
alignment of the eigenmodes of the target layer relative to the first two
eigenmodes of the input layer. In particular, the alignment between the most
dominant eigenmodes of the two layers is predictably altered by changing the
continuous parameter $\theta$ (Fig. 3(A) and the supplementary section  S8).

How does the alignment of the dominant eigenmode in the target layer with that
of the input layer affect the optimal cost of control? We hypothesize that
when the target eigenmode is well-aligned with the dominant eigenmode of the
input layer, the corresponding control energy will be small. As a corollary,
the optimal control signal will be channeled through the eigenmode of the
input layer that is aligned to the dominant eigenmode of the target layer.
Thus, as the dominant eigenmode of the target layer changes its alignment with
the dominant eigenmode of the input layer, from a perfect alignment to an
orthogonal orientation, the optimal control signal should excite an increasing
number of eigenmodes in the input layer. To quantify this effect, we define
$\mathcal{N}(\phi)$, which counts the number of eigenmodes
$\mathbf{p}_{\alpha}$ of the input layer such that $|u_{\alpha}|^{2}>10^{-3}$.
Formally,

$\displaystyle\mathcal{N}(\phi)=|\\{\mathbf{p}_{\alpha}:\quad|u_{\alpha}|^{2}>10^{-3}\\}|.$
(11)

We say that $\mathcal{N}(\phi)$ is the number of eigenmodes
$\mathbf{p}_{\alpha}$ of the input layer through which the optimal control
signal is ‘routed’.

For each duplex, we calculate the optimal energy to control the dominant
eigenmode and $\mathcal{N}(\phi)$, as a function of the angle $\phi$ (see the
supplementary section and Fig.  S8). Note that, since the rotation is defined
in the plane of the first two eigenvectors of the input layer, the alignment
of the target layer relative to the remaining $N-2$ eigenvectors that form a ‘
_rotation axis_ ’, remains unchanged. The optimal control energy as well as
the number of channels that are excited in the input layer to control the
dominant eigenvector of the target layer decrease as its alignment with the
dominant eigenvector of the input layer decreases (Fig. 3(B-C)). These results
imply that the alignment of a target eigenvector relative to the eigenvectors
of the input layer plays an important role in determining the optimal control
cost, and the eigenvectors that are excited in the process of control.

Note that, even as the rotation of the target eigenmode decreases the
alignment with the first axis of the rotation plane, and enhances the
alignment with the second axis of the rotation plane, we do not observe a
decrease in energy or the number of paths. To explain this point further, as
the rotating target eigenmode becomes more aligned with the second axis in the
rotation plane, the optimal energy and the number of channels must decrease as
the optimal energy must now be routed via the second axis. However, this
effect is not observed for two reasons. First, in addition to the alignment,
control energy also depends on the corresponding eigenvalue. In the analysis
of the next section, we find that a more unstable eigenmode is associated with
a lower energy. Thus, the second axis that corresponds to the second most
unstable eigenmode has higher cost associated with it. Second, the plots of
the fraction of energy carried by the first two eigenmodes show that even when
the target eigenmode is aligned with the second axis, the fraction of energy
routed along this axis is less than one (Fig. 3(D-F)). This effect is
particularly prominent in WS and RG networks. These observations indicate that
the distribution of optimal control energy, while strongly influenced by the
alignment, is also affected by the topology and the eigenvalues of the
corresponding eigenmodes.

In an LTI system with the dynamics in continuous time as considered here,
eigenvalues of the adjacency matrix have the unit of a ‘rate’. Thus, an
eigenvalue corresponding to a given eigenmode determines the time scale at
which the perturbations in the corresponding eigenmodes decay or grow. The
results in this section indicate that the eigenvalues, and the alignment
between eigenmodes, together may determine the optimal control energy. In the
next section, we attempt to understand this point further in the simplified
limit of a one-mode approximation.

##  IV Intuitions derived from a one-mode approximation

Figure 4:  Understanding structure-control relationship in single-mode limit.
(A) Control input directly couples to the input layer and can excite an
arbitrary number of modes. In a one-mode approximation, we assume that the
response of a layer is dominated by a single mode. (B) Control energy required
for navigating the system to the $\\{0,1\\}$ state as a function both of
$\mu_{\beta}/\xi_{\alpha}$ and of the overlap between retained eigenmodes. (C)
Control energy to navigate the system to a $\\{0,1\\}$ state as a function of
$\mu_{\beta}/\xi_{\alpha}$. (D) A faster decaying mode costs more control
energy whereas a faster growing mode is easier to control.

To gain some insight regarding the dependence of optimal control energy on
layer spectra and the interlayer alignment, we consider a simplified limit
where we retain only one mode in each layer (Fig. 4(A)). Studying this limit
can have practical applications when the phenomenon of interest happens at a
specific time-scale or spatial scale. In real networks, several dynamical
processes are dominated by extremal eigenvalues, such as the stability of
dynamical systems on networks, synchronization, collective behaviour, and
epidemic spreading Sarkar_Chaos_2018 ; Yang_epidemic_2003 ; Newman_book_2010 .
Here we operationalize this simplification by revisiting Eqs. (7a) and (7b).
We retain the $\alpha^{\rm{th}}$ eigenmode in the input layer and the
$\beta^{\rm{th}}$ eigenmode in the target layer to obtain the simplified
equations:

$\displaystyle\dot{w}_{\alpha}$
$\displaystyle=\xi_{\alpha}w_{\alpha}+U_{\alpha}$ (12a)
$\displaystyle\dot{v}_{\beta}$
$\displaystyle=\mu_{\beta}v_{\beta}+\kappa\cos\theta_{\beta\alpha}w_{\alpha}.$
(12b)

Applying the framework of optimal control framework in this effectively one-
dimensional representation of each layer, we obtain the solutions for the
optimal state trajectory, the optimal co-state trajectory, and the optimal
control energy (see the supplementary section  S9). These optimal solutions
depend upon the model parameters: the eigenvalues $\xi_{\alpha},\mu_{\beta}$,
the overlap $C_{\beta\alpha}=\cos(\theta_{\beta\alpha})$ between the retained
eigenmodes, and the strength $\kappa$ of the interlayer connections. We
calculate the optimal control energy $E^{\ast}$ to reach the states
$\\{w_{\alpha F},v_{\beta F}\\}$, where $w|_{t=T}=w_{\alpha F}$ and
$v|_{t=T}=v_{F}$.

In order to separately calculate the control energy required to navigate the
states of the two layers away from the origin, we set the final states to be
$\\{w_{\alpha F},v_{\beta F}\\}=\\{1,0\\}$ and $\\{w_{\alpha F},v_{\beta
F}\\}=\\{0,1\\}$. We then consider the optimal control energy as a function of
the eigenvalues $\xi_{\alpha}$ and $\mu_{\beta}$ ( Fig. 4(C) and the
supplementary section  S9). We observe that a faster decay rate makes it
harder to control the state of a given layer, while a faster growth rate makes
it easier to control the state of that layer along the retained eigenvector.
This behavior is intuitive: the faster the decay, the harder it is to navigate
the state away from the origin (Fig. 4(D)).

In the one-mode limit, we also calculate the dependence of the optimal control
energy on the alignment parameter $C_{\beta\alpha}$. The dependence of control
energy on $C_{\beta\alpha}$ has the form
$E^{\ast}_{\alpha}=S_{-2}C_{\beta\alpha}^{-2}+S_{-1}C_{\beta\alpha}^{-1}+S_{0}$
(see the supplementary section  S9 for detailed expressions). Here, the
coefficients $S_{-2},S_{-1}$ and $S_{0}$ depend on the eigenvalues
corresponding to the retained eigenmodes of the two layers, and the final
states $w_{F}$ and $v_{F}$. It is intuitive to note the dependence of the
above coefficients on the final states $\\{w_{F},v_{F}\\}$ given by
$S_{-2}\approx v_{\beta F}^{2},\quad S_{-1}\approx v_{\beta F}w_{\alpha F}$
and $S_{0}\approx w_{\alpha F}^{2}$. Note that when $v_{\beta F}\neq 0$,
$E\approx C_{\beta\alpha}^{-2}$. Thus, the energy to control the state of the
target layer when the control signal is injected into the input layer diverges
when $C_{\beta\alpha}\rightarrow 0$. Intuitively, it is therefore impossible
to control a mode by injecting control into a mode that is orthogonal to it.
Also we note that, when $v_{\beta F}=0$ but $w_{\alpha F}\neq 0$,
$E^{\ast}_{\alpha}$ becomes independent of $C_{\beta\alpha}$; this
independence is expected because the control is directly injected into the
mode that is also being controlled. Additionally, $E^{\ast}_{\alpha}$
inevitably depends on the alignment parameter $C_{\beta\alpha}$ when $v_{\beta
F}\neq 0$; that is, the control of specific eigenmodes in the target layer is
always affected by their alignment to those in the input layer. We summarize
these results in Fig. 4(B) by plotting the optimal control energy as a
function of the alignment $C_{\beta\alpha}$ and the eigenvalues $\xi_{\alpha}$
and $\mu_{\beta}$.

The decrease in optimal control energy with increasing alignment is a result
that we also obtained numerically in Section III (Fig. 3(B)). We note one
important difference. Unlike these new observations in the one-mode limit, our
numerical experiments did not show a divergence of the optimal control energy
when the alignment approaches zero. This discrepancy arises from the fact that
in numerically constructed networks, the eigenspaces are generally misaligned.
This fact implies that even when the target eigenmode is rotated in the plane
of the first two eigenvectors of the input layer, it has non-zero overlap with
the remaining $N-2$ eigenvectors of the input layer, which remain unaltered by
the rotation. As a result, the optimal control energy is ‘routed’ (in the
sense that is defined in Section III) through those eigenmodes that have non-
zero overlap with the target eigenmode.

Finally, we ask whether the optimal control energy for the numerically
constructed duplex networks shows association with the eigenvalues of the
target eigenmodes in the second layer that are qualitatively similar to that
seen in the one-mode approximation, as depicted in Fig. 4(C). For all the
duplex topologies and layer densities, we find an association between the
control energy of the target layer and the corresponding eigenvalues similar
to that obtained in the one-mode limit (see the supplementary section and Fig.
S10). This behavior is also generically observed for the control energy of the
input layers, except in some specific cases when input layer has a WS
topology. Collectively these findings indicate that we can gain significant
understanding of the broad relation between control and topology by performing
an analytic investigation of the simplified one-mode approximation.

Figure 5: Distribution of control energy in the eigenmodes of the input layer.
(A) The distribution of optimal control energy into the eigenvectors of the
input layer for duplex networks with identical layers. Here we construct the
final state by linearly combining randomly selected target-layer eigenmodes.
We then ask whether the modes that are excited in the input layer by the
optimal control signals are identical to those that constitute the target
state. (B) The distribution of optimal control energy into the eigenvectors of
the input layer for duplex networks with non-identical layers. Here, the final
state is selected to be along the eigenvectors of the target layer, and the
time-averaged energy along each eigenvector of the input layer is calculated.
(C) In the duplex network with identical layers, the eigenmodes excited in the
input layer are identical to those that constitute the final state of the
target layer (indicated by the green dots on $x$-axis). The instantaneous
energy along the $\beta^{\rm{th}}$ eigenmode of the input layer corresponding
to the $\alpha^{\rm{th}}$ eigenmode of the target layer is defined as
$|\omega_{21}^{\alpha\beta}(t)|^{2}$. The figure shows that
$|\omega_{21}^{\alpha\beta}(t)|^{2}$ (gray data points) is non-vanishing only
for eigenmodes identical to those that constitute the final state. (D) For all
of the duplex networks, the mean energy along an eigenvector of the input
layer to control the state of the target layer along its eigenvectors is
plotted as a function of the alignment between input eigenvectors and target
eigenvectors. The mean energy is normalized by its maximum value across all of
the input eigenvectors for a given target eigenvector. This normalization is
performed in order to highlight the role of alignment. A strong correlation
between the two quantities emerges, indicating the important role that
alignment plays in control.

##  V Path-dependence of the control energy

The results obtained in the previous sections for the numerically constructed
duplex networks and for the one-mode approximation indicate an important role
of the alignment – between the modes of the target layer and the modes of the
input layer – in determining the cost of control. In particular, the number of
modes that are excited in the input layer increase as the alignment between
the target eigenmode and the eigenmodes of the input layer decrease. In this
section, we formalize the notion of ‘routing’ the optimal control signal
further by building metrics to identify the eigenmodes that are excited in the
input layer by the control signal optimized to move the state of the target
layer along a given eigenvector.

We begin by denoting the optimal input for controlling the state of the
$\alpha^{\rm{th}}$ eigenvector of the $l^{\rm{th}}$ layer by
$u_{l}^{\alpha}(t)$. To measure how $u_{l}^{\alpha}$ is distributed in the
eigenmodes of the $m^{\rm{th}}$ layer, we project $u_{l}^{\alpha}$ along the
eigenvectors (labelled by $\beta$) of the $m^{\rm{th}}$ layer:

$\displaystyle\mathbf{u}_{l}^{\alpha}(t)=\Sigma\omega_{lm}^{\alpha\beta}(t)\mathbf{V}_{m}^{\beta},$
(13)

where $\omega_{\beta}(t)$ is the projection of the optimal input
$\mathbf{u}^{\alpha}(t)$. For a duplex network, $l,m\in\\{1,2\\}$. Then,
$l=1,m=1$ corresponds to a case where the $\alpha^{\rm{th}}$ eigenmode in the
input layer is the target eigenmode, and the corresponding optimal input
$u^{\alpha}_{1}(t)$ is projected along the eigenvectors of that eigenmode.
Thus, the parameters $\omega_{11}^{\beta}(t),\omega_{22}^{\beta}(t)$ and
$\omega_{21}^{\beta}(t)$ provide access to the directional information of the
optimal $\mathbf{u}^{\alpha}(t)$. Further, the quantity
$|\omega_{lm}^{\beta}(t)|^{2}$ can be identified as the energy channeled into
the $\beta^{\rm{th}}$ eigenmode at time $t$.

The quantities $\omega^{\beta}_{lm}(t)$ can be utilized to identify the modes
that are excited by the optimal control signal corresponding to a generic
state in the state-space (i.e. not necessarily along the eigenmodes). From the
insights gained in the previous sections on the role of alignment, we
hypothesize that in a duplex network with identical layers, and hence
perfectly aligned eigenmodes, the optimal control for a given eigenmode in the
target layer will excite the identical eigenmode in the input layer. We
further hypothesize that in duplex networks with non-identical layers, the
fraction of the control energy along a given eigenmode of the input layer will
be proportional to its alignment with the target eigenmode in the second
layer. The above will be true if mode alignment is the only factor that
determines the direction of optimal control. In the following we show that
while the first hypothesis holds true, a test of the latter hypothesis reveals
an additional (and admittedly expected) dependence on the topology of the
duplex.

We first construct duplex networks with identical layers from all 4 graph
families. We then construct a generic final state using a finite number of
randomly selected eigenmodes of the target layer (Fig. 5(A)). We observe that
the optimal control signal excites only those eigenmodes in the input layer
that are identical to those forming the target state (Fig. 5(C)). This
observation verifies the first hypothesis. To test the second hypothesis and
to determine the fraction of control energy along a given eigendirection of
the input layer, we calculate a time-averaged
$\bar{E}_{\alpha\beta}=|\omega^{\alpha\beta}_{21}(t)|^{2}$, for all
$\alpha\in\\{1,N\\}$ (eigenmodes of the target layer) and for all
$\beta\in\\{1,N\\}$ (eigenmodes of the input layer). The time-averaged
$\bar{E}_{\alpha\beta}$ then represents the mean energy along the $\alpha$
eigenvector corresponding to the control of the target layer along the
$\beta^{\rm{th}}$ eigenvector. We assess the correlation between
$\bar{E}_{\alpha\beta}$ and the alignment of the target eigenmode with the
eigenmodes of the input layer for all the constructed duplex networks and all
densities (Fig. 5(C)). We find that the alignment indeed emerges as a strong
factor in determining the fraction of control energy along an eigenvector of
the input layer. This fact is evidenced by the presence of a strong trend of
$\bar{E}_{\alpha\beta}\propto f(|C_{\alpha\beta}|)$ where $f$ is an increasing
function of the alignment $|C_{\alpha\beta}|$. However, deviations from this
behavior appear especially for higher densities (Fig. 5(C)), and when the
target layer has the topology of a BA or an RG network.

##  VI Conclusion and Future Directions

In our investigation, we have sought to better understand the relationship
between structure and the cost of control, in a principled and theoretically
motivated way. Our efforts, which combine analytical results and numerical
experiments, allow us to build a generic framework to extract the dependence
of multiplex control on layered architecture. We find that the response rates
as characterized by the eigenspectra of each layer, and the alignment between
the eigenmodes of the input and target layers emerge as key parameters that
determine the cost of optimal control as well as the distribution of optimal
control energy in different channels. While results from a one-mode
approximation provide some useful insight regarding the correlation between
eigenvalues and control cost of individual eigenmodes, the global
characteristics of overall network topology (such as density, degree
distribution, and clustering coefficient) determine the specific manner in
which optimal control energy is distributed in the eigenmodes of the two
layers. The effect of the global topology is particularly prominent in RG and
WS networks, both of which show high values of average clustering coefficient.
Our work indicates that a high value of average clustering coefficient is
associated with lower average control energy while a high value of degree
heterogeneity is associated with the maximum control energy. In addition, we
also studied the distribution of optimal control energy into the eigenmodes of
the input layer in interlayer control. Our approach can be utilized in
identifying and predicting the paths of optimal control in specific
applications.

In summary, we investigate the dependence of control properties of layers in
duplex networks as they are dictated by layer spectra and the alignment
between different eigenmodes of the two layers. We establish the key role of
mode alignment in controlling the state of the target layer along specific
eigendirections and further highlight the interplay between alignment and
layer topologies in determining the cost and routing of optimal control
energies. Our approach of determining the optimal energies for state
transitions in a _layer-diagonalised_ representation can be easily generalized
to multiplex networks with more than two layers, and to more generic
interlayer connections. The specific manner in which a more generic pattern of
interlayer connections interplays with the relative alignment between the
input and the target layers is an exciting direction of future research in the
theory of multiplex control. Our results will serve as a platform for future
work in narrowing the space of feasible intervention strategies in
applications where cross-layer control is required. Indeed, recent multilayer
representations of networks in biology have already indicated that different
connectivity layers can be misaligned Becker_SciRep2018 , impacting the
system’s potential response to perturbations. Using our framework, future work
can explore the consequences of such layered architectures on designing
optimal interlayer control strategies.

##  VII Citation Diversity Statement

Recent work in several fields has identified a bias in citation practices such
that papers from women and other minorities are under-cited relative to the
number of such papers in the field (Dworkin2020.01.03.894378, ;
maliniak2013gender, ; caplar2017quantitative, ;
chakravartty2018communicationsowhite, ;
YannikThiemKrisF.SealeyAmyE.FerrerAdrielM.Trott2018, ; dion2018gendered, ).
Here we sought to proactively consider choosing references that reflect the
diversity of the field in thought, gender, race, ethnicity, and other factors.
We obtained the predicted gender of the first and last author of each
reference by using databases that store the probability of a first name being
carried by a woman (Dworkin2020.01.03.894378, ). By this measure (and
excluding self-citations to the first and last authors of our current paper),
our references contain $60\%$ man/man, $20\%$ man/woman, $4.4\%$ woman/man,
$15.6\%$ woman/woman, and $0\%$ unknown categorization. Open-source code for
this estimation process is publicly available (github_dale, ). This method is
limited in that a) names, pronouns, and social media profiles used to
construct the databases may not, in every case, be indicative of gender
identity and b) it cannot account for intersex, non-binary, or transgender
people. We look forward to future work that could help us to better understand
how to support equitable practices in science.

###### Acknowledgements.

This work was primarily supported by the Schwartz Foundation, the Army
Research Office (Falk-W911NF-18-1-0244), and the Paul G. Allen Family
Foundation. We would also like to acknowledge additional support from the John
D. and Catherine T. MacArthur Foundation, the Alfred P. Sloan Foundation, the
ISI Foundation, the Army Research Office (Grafton-W911NF-16-1-0474,
DCIST-W911NF-17-2-0181), the National Institute of Mental Health
(2-R01-DC-009209-11, R01-MH112847, R01-MH107235, R21-M MH-106799), and the
National Science Foundation (NSF PHY-1554488, BCS-1631550, and IIS-1926757).
The content is solely the responsibility of the authors and does not
necessarily represent the official views of any of the funding agencies.

## References

  * (1) M. Gosak, R. Markovič, J. Dolenšek, M. S. Rupnik, M. Marhl, A. Stožer, and M. Perc, “Network science of biological systems at different scales: A review,” Physics of Life Reviews, vol. 24, pp. 118–135, 03 2018.
  * (2) D. S. Bassett, P. Zurn, and J. I. Gold, “On the nature and use of models in network neuroscience,” Nat Rev Neurosci, vol. 19, no. 9, pp. 566–578, 2018\.
  * (3) L. Papadopoulos, M. A. Porter, K. E. Daniels, and D. S. Bassett, “Network analysis of particles and grains,” Journal of Complex Networks, vol. 6, no. 4, p. 485–565, 2018.
  * (4) M. O. Jackson and A. Watts, “The evolution of social and economic networks,” Journal of Economic Theory, vol. 106, no. 2, pp. 265 – 295, 2002.
  * (5) M. O. Jackson, Social and Economic Networks.  Princeton University Press, 2008.
  * (6) L. Torres, A. S. Blevins, D. S. Bassett, and T. Eliassi-Rad, “The why, how, and when of representations for complex system,” arXiv, vol. 2006, p. 02870, 2020.
  * (7) G. Bianconi, Multilayer Networks: Structure and Function.  Oxford University Press, 2018.
  * (8) M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, “Multilayer networks,” Journal of Complex Networks, vol. 2, no. 3, pp. 203–271, 2014.
  * (9) R. Bookstaber and D. Y. Kenett, “Looking deeper, seeing more: A multilayer map of the financial system,” Office of Financial Research Brief Series, 2016\.
  * (10) D. Wójcik, The global financial network.  2018\.
  * (11) M. E. Dickison, M. Magnani, and L. Rossi, Multilayer Social Networks.  Cambridge University Press, 2016.
  * (12) C. Atkisson, P. J. Górski, M. O. Jackson, J. A. Holyst, and R. M. D’Souza, “Why understanding multiplex social network structuring processes will help us better understand the evolution of human behavior,” arXiv, vol. 11183, 2019.
  * (13) J. M. Buldú and M. A. Porter, “Frequency-based brain networks: From a multiplex framework to a full multilayer description,” Network Neuroscience, vol. 2, no. 4, pp. 418–441, 2017.
  * (14) J. M. Buldú and D. Papo, “Can multilayer brain networks be a real step forward?,” Physics of Life Reviews, vol. 24, pp. 153–155, 2018.
  * (15) S. F. Muldoon, “Multilayer network modeling creates opportunities for novel network statistics: Comment on ”network science of biological sysetms at different scales: A review” by gosak et al.,” Phys Life Rev., vol. 24, pp. 143–145, 2018.
  * (16) S. Pilosof, M. A. Porter, M. Pascual, and S. Kéfi, “The multilayer nature of ecological networks,” Nature Ecology and Evolution, 03 2017.
  * (17) M. De Domenico, A. Solé-Ribalta, E. Cozzo, M. Kivelä, Y. Moreno, M. A. Porter, S. Gómez, and A. Arenas, “Mathematical formulation of multilayer networks,” Physical Review X, vol. 3, no. 4, p. 041022, 2013.
  * (18) B. Bentley, R. Branicky, C. L. Barnes, Y. L. Chew, E. Yemini, E. T. Bullmore, P. E. Vértes, and W. R. Schafer, “The multilayer connectome of Caenorhabditis elegans,” PLoS Comput Biol, vol. 12, no. 12, p. e1005283, 2016.
  * (19) B. A. Kidd, L. A. Peters, E. E. Schadt, and J. T. Dudley, “Unifying immunology with informatics and multiscale biology,” Nat Immunol, vol. 15, no. 2, pp. 118–127, 2014.
  * (20) M. Schmidt, R. Bakker, C. C. Hilgetag, M. Diesmann, and S. J. van Albada, “Multi-scale account of the network structure of macaque visual cortex,” Brain Structure and Function, vol. 223, no. 3, pp. 1409–1435, 2018.
  * (21) D. S. Bassett and M. S. Gazzaniga, “Understanding complexity in the human brain,” Trends Cogn Sci, vol. 15, no. 5, pp. 200–209, 2011.
  * (22) U. Braun, A. Schaefer, R. F. Betzel, H. Tost, A. Meyer-Lindenberg, and D. S. Bassett, “From maps to multi-dimensional network mechanisms of mental disorders,” Neuron, vol. 97, no. 1, pp. 14–31, 2018.
  * (23) D. Magnusson and R. B. Cairns, Developmental science: Toward a unified framework, pp. 7–30.  Cambridge University Press, 1996.
  * (24) U. Bronfenbrenner, The Ecology of Human Development.  Cambridge, MA: Harvard University Press, 1979.
  * (25) T. Kailath, Linear Systems.  Prentice-Hall, 1980.
  * (26) R. Kubo, M. Toda, and N. Hashitsume, “Statistical mechanics of linear response,” in Statistical Physics II. Springer Series in Solid-State Sciences, vol. 31, Springer, Berlin, Heidelberg.
  * (27) J. Z. Kim, J. M. Soffer, A. E. Kahn, J. M. Vettel, F. Pasqualetti, and D. S. Bassett, “Role of graph architecture in controlling dynamical networks with applications to neural systems,” Nature Physics, vol. 14, pp. 91–98, 2018\.
  * (28) E. Tang and D. S. Bassett, “Colloquium: Control of dynamics in brain networks,” Review of Modern Physics, vol. 90, p. 031003, 08 2018.
  * (29) M. Pósfai, J. Gao, S. P. Cornelius, A.-L. Barabasi, and R. M. D’Souza, “Controllability of multiplex, multi-time-scale networks,” Phys Rev E, vol. 94, no. 3, p. 032316, 2016.
  * (30) G. Menichetti, L. Dall’Asta, and G. Bianconi, “Control of multilayer networks,” Scientific Reports, vol. 6, no. 20706, 2016.
  * (31) Z. Yuan, C. Zhao, W.-X. Wang, Z. Di, and Y.-C. Lai, “Exact controllability of multiplex networks,” New Journal of Physics, vol. 16, p. 103036, 2014.
  * (32) D. Wang and X. Zou, “Control energy and controllability of multilayer networks,” Advances in Complex Systems, vol. 20, no. 0, p. 1750008, 2017\.
  * (33) M. A. Wicks and R. A. DeCarlo, “An energy approach to controllability,” Proceedings of the 27th conference on Decision and Control, pp. 2072–2077, 1988\.
  * (34) F. L. Lewis, D. L. Vrabie, and V. L. Syrmos, Optimal Control.  John Wiley & Sons, 2012.
  * (35) A.-L. Barbási and R. Albert, “Emergence of scaling in random networks,” Science, vol. 286, pp. 509–512, 1999.
  * (36) R. Albert and A.-L. Barbási, “Statistical mechanics of complex networks,” Reviews of Modern Physics, vol. 74, pp. 47–97, 01 2002.
  * (37) M. Penrose, Random Geoemtric Graphs.  Oxford University Press, 2004.
  * (38) D. J. Watts and S. H. Strogatz, “Collective dynamics of small-world networks,” Nature, vol. 393, pp. 440–442, 1998.
  * (39) R. Albert, “Scale-free networks in cell biology,” Journal of Cell Science, vol. 118, pp. 4947–4957, 2005.
  * (40) E. F. Keller, “Revisiting “scale-free” networks,” BioEssays, vol. 27, 2005\.
  * (41) Q. Telesford, K. E. Joyce, S. Hayasaka, J. H. Burdette, and P. Laurienti, “The ubiquity of small-world networks,” Brain Connectivity, vol. 1, no. 5, 2011\.
  * (42) C. W. Lynn, L. Papadopoulos, A. E. Kahn, and D. S. Bassett, “Human information processing in complex networks,” Nature Physics, vol. 16, pp. 965–973, 2020.
  * (43) C. Doran and A. Lasenby, Geometric Algebra for Physicists.  Cambridge University Press, 2003.
  * (44) C. Sarkar and S. Jalan, “Spectral properties of complex networks,” Chaos, p. 102101, 2018.
  * (45) Y. Wang, D. Chakrabarti, C. Wang, and C. Faloutsos, “Epidemic spreading in real networks: an eigenvalue viewpoint,” Proceedings of 22nd International Symposium on Reliable Distributed Systems, pp. 25–34, 2003.
  * (46) M. Newman, Networks: an introduction.  Oxford University Press, 2010.
  * (47) C. O. Becker, S. Pequito, G. J. Pappas, M. B. Miller, S. T. Grafton, D. S. Bassett, and V. M. Preciado, “Spectral mapping of brain functional connectivity from diffusion imaging,” Scientific Reports, no. 1411, 2018\.
  * (48) J. D. Dworkin, K. A. Linn, E. G. Teich, P. Zurn, R. T. Shinohara, and D. S. Bassett, “The extent and drivers of gender imbalance in neuroscience reference lists,” bioRxiv, 2020.
  * (49) D. Maliniak, R. Powers, and B. F. Walter, “The gender citation gap in international relations,” International Organization, vol. 67, no. 4, pp. 889–922, 2013.
  * (50) N. Caplar, S. Tacchella, and S. Birrer, “Quantitative evaluation of gender bias in astronomical publications from citation counts,” Nature Astronomy, vol. 1, no. 6, p. 0141, 2017.
  * (51) P. Chakravartty, R. Kuo, V. Grubbs, and C. McIlwain, “# communicationsowhite,” Journal of Communication, vol. 68, no. 2, pp. 254–266, 2018.
  * (52) Y. Thiem, K. F. Sealey, A. E. Ferrer, A. M. Trott, and R. Kennison, “Just Ideas? The Status and Future of Publication Ethics in Philosophy: A White Paper,” tech. rep., 2018.
  * (53) M. L. Dion, J. L. Sumner, and S. M. Mitchell, “Gendered citation patterns across political science and social science methodology fields,” Political Analysis, vol. 26, no. 3, pp. 312–327, 2018.
  * (54) D. Zhou, E. J. Cornblath, J. Stiso, E. G. Teich, J. D. Dworkin, A. S. Blevins, and D. S. Bassett, “Gender diversity statement and code notebook v1.0,” Feb. 2020.