# Impact of shocks to economies on the efficiency and robustness of the
international pesticide trade networks

Jian-An Li  Li Wang  Wen-Jie Xie<EMAIL_ADDRESS>Wei-Xing Zhou  School of
Business, East China University of Science and Technology, Shanghai 200237,
China  School of Physics, East China University of Science and Technology,
Shanghai 200237, China  Research Center for Econophysics, East China
University of Science and Technology, Shanghai 200237, China  School of
Mathematics, East China University of Science and Technology, Shanghai 200237,
China

###### Abstract

Pesticides are important agricultural inputs to increase agricultural
productivity and improve food security. The availability of pesticides is
partially achieved through international trade. However, economies involved in
the international trade of pesticides are impacted by internal and external
shocks from time to time, which influence the redistribution efficiency of
pesticides all over the world. In this work, we adopt simulations to quantify
the efficiency and robustness of the international pesticide trade networks
under shocks to economies. Shocks are simulated based on nine node metrics,
and three strategies are utilized based on descending, random, and ascending
node removal. It is found that the efficiency and robustness of the
international trade networks of pesticides increased for all the node metrics
except the clustering coefficient. Moreover, the international pesticide trade
networks are more fragile when import-oriented economies are affected by
shocks.

††journal: Eur. Phys. J. B

##  1 Introduction

Pesticides include insecticides, fungicides, herbicides, disinfectants, and
rodenticides and other similar products, which are invented to protect
agricultural crops from harms caused by insects, fungi, weeds, viruses, and
rats. Therefore, the main functions of pesticides are to reduce yield losses,
regulate plant growth, and increase agricultural productivity, which is
essential to improving food security. The availability of pesticides in most
economies is partially supported by international trade. However, there are
internal and external shocks from time to time that affect the involved
economies and their ability to trade pesticides internationally. Such shocks
to economies influence the efficiency of pesticides redistribution all over
the world. It is, thus, important to quantify the efficiency and robustness of
the international pesticide trade networks (iPTNs).

The structural properties and their evolutionary behavior of the iPTNs have
been studied [1, 2]. Moreover, the structural robustness and efficiency-based
robustness of the iPTNs have also been investigated when the trade
relationships are affected by internal and external shocks [3]. When the
economies involved in the international pesticide trade are influenced by
shocks, the structural robustness of the iPTNs has been quantified [4]. In
this work, we aim to complete the unfinished puzzle by quantifying the
efficiency and efficiency-based robustness of the iPTNs under shocks to
economies. We note that the efficiency-based robustness of complex networks
are less studied [5, 6].

Concerning the structural robustness of most complex networks, accumulating
evidence shows that complex networks are robust to internal shocks (or random
failures) but may be fragile to external shocks (or intentional attacks) [7,
8, 9, 10]. Researchers have performed extensive analysis of the robustness and
fragility of complex networks in various fields, such as international oil
trade networks [11, 12, 13, 14, 15, 16], infrastructure networks [17, 18, 19,
20, 21, 22, 23, 24], and the networks of Cosa Nostra affiliates [25], to list
a few. To carry out simulations for quantifying the efficiency and robustness
of the iPTNs under shocks to economies, we utilize nine node metrics
(clustering coefficient; betweenness; in-degree, PageRank, authority, and in-
closeness; out-degree, hub, and out-closeness) and three node removal
strategies (descending, random, and ascending).

The remainder of this work is organized as follows. We briefly describe the
data in Section 2. We quantify the network efficiency of iPTNs in Section 3
and the efficiency-based robustness of iPTNs in Section 4. We summarize our
results in Section 5.

##  2 Data description

The data sets for the international pesticide trade were retrieved from the UN
Comtrade database (publicly available at https://comtrade.un.org), which
covers the period from 2007 to 2018 and contains the import and export
economies and trading values of five categories, including insecticides
(380891), fungicides (380892), herbicides (380893), disinfectants (380894),
and rodenticides and other similar products (380899).

Based on the data sets, for each year, we construct international pesticide
trade networks for each category of pesticides, where involved economies are
nodes and a directed link is drawn when two economies trade [13, 3].

##  3 Network efficiency

###  3.1 Shortest path length between economies

Figure 1:  Distance distribution among economies in the international
pesticide trade network. $d$ is the network distance between economies in the
directed pesticide trade network. $P(d)$ represents the proportion of the
shortest path with a distance of $d$. The figure shows the time evolution of
the proportion of the distance $d=1$, 2, 3, 4, 5, and $d>5$. The six plots
correspond to the trade networks of pesticide products, including (a)
aggregated, (b) insecticides, (c) fungicides, (d) herbicides, (e)
disinfectants, and (f) rodenticides and other similar products.

In the international pesticide trade networks, pesticide products flow between
economies. There are a lot of pesticide-producing and pesticide-consuming
economies in the trade network, and there are also some economies with a lot
of pesticide imports and exports. The network efficiency analysis in this
paper is based on directed networks, and the relevant metric indicators are
for asymmetric matrices, except for the clustering coefficient. Network
efficiency is also based on directed shortest paths. In the international
pesticide trade system, different economies play different roles and
contribute their own strengths to the global pesticide trade and food
production. It is difficult to distinguish between exporting and importing
economies because the international pesticide trade networks include a large
number of economies that both export and import. Network efficiency refers to
the flow efficiency of information, energy, and matter in the network, such as
the information transmission efficiency of the information network and the
traffic efficiency of the traffic network [26, 27]. Generally speaking, the
trade network efficiency can be used as the path length between the producer
and the consumer. The larger the path length, the lower the efficiency;
otherwise, the higher the efficiency.

We calculated the shortest path length between economies in the iPTNs. Figure
1 shows the evolution of occurrence frequencies $P(d)$ for $d=1$, $d=2$,
$d=3$, $d=4$, $d=5$, and $d>5$, where $d$ is the distance between the
economies in the pesticide trade network.

Overall, the occurrences of shortest path length in the six plots are
relatively similar, especially in plot (c) for fungicides, plot (d) for
herbicides, plot (e) for disinfectants, and plot (f) for rodenticides and
other similar products. In each of these four temporal iPTNs, the proportion
$P(d)$ with $d>5$ is basically the largest, which is close to 0.5 in 2007, but
most of them have an infinite shortest path length. It is due to the fact that
the iPTNs we are studying are directed and thus there is no reachable path
between some economies. With the development of globalization, there are more
and more trade relations, and the proportion of economy pairs with infinite
shortest path lengths is getting smaller and smaller until around 2015. The
main reason is that the sales in the international pesticide market plunged
8.5% year on year, the steepest decline in more than a decade. In most years,
we have

$P(2)>P(3)>P(1)>P(4)>P(5).$ (1)

This relation does not hold for the aggregated iPTN in Fig. 1(a) and the iPTN
for insecticides in Fig 1(b). The aggregated iPTN integrates the five
pesticide products and has the largest number of links. Moreover, compared
with the other four pesticide trade products, the insecticides trade network
showed the largest trade volume and the largest trade relationships.

The proportion $P(1)$ of distance $d=1$ in Fig. 1 measures the density of the
pesticide trade network. It can be seen from the figure that, with the
evolution of time $t$, $P(1)$ becomes larger and larger, indicating that there
are more and more pesticide trade relations, but after 2015, there is also a
downward trend.

The proportion $P(2)$ of distance $d=2$ measures that the economies in the
pesticide trade network do not have direct trade relations, but have indirect
trade relations through intermediate economies. It can be seen from Fig. 1
that, with the evolution of time $t$, $P(2)$ becomes larger and larger,
indicating that there are more and more pesticide trade relations, but after
2015, there is also a downward trend. In particular, the $P(2)$ curves for the
aggregated iPTN in Fig. 1(a), the iPTN of insecticides in Fig. 1(b), and the
iPTN of fungicides in Fig. 1(c) showed obvious fluctuations in 2015.

Generally speaking, the shorter the trade distance between economies in the
international trade network, the higher the trade efficiency. However, the
international trade network is affected by geographical location, political
factors, resource endowment distribution, etc. It is impossible that the
shortest path length of the trade network is too small. Therefore, we
introduced the index of network efficiency to measure the trade flow
efficiency and the change rule of the trade structure of pesticide products in
the international pesticide trade network.

###  3.2 Network efficiency

Network efficiency is defined as the average reciprocal path length between
nodes in the network. The path length between nodes is defined differently
between directed networks and undirected networks, as well as between weighted
networks and unweighted networks. In an unweighted, directed network, the
network efficiency $E$ is defined as follows [26, 27],

$E=\frac{1}{N(N-1)}\sum_{i\neq j}e_{ij},$ (2)

where $e_{ij}$ represents the path efficiency between economy $i$ and economy
$j$, measuring the transmission efficiency of material flow, information flow,
or energy flow between two nodes, specifically defined as

$e_{ij}=\frac{1}{d_{ij}},$ (3)

where $d_{ij}$ represents the shortest path length between economy $i$ and
economy $j$.

Figure 2: Evolution of the efficiency of the international pesticide trade
networks. Each solid line corresponds to an international pesticide trade
network, including aggregated, insecticides, fungicides, herbicides,
disinfectants and rodenticides and other similar products.

Fig. 2 shows the evolution of the efficiency of the international pesticide
trade networks. Each solid line corresponds to an international pesticide
trade network, including aggregated, insecticides, fungicides, herbicides,
disinfectants and rodenticides and other similar products. The five types of
pesticide trade network efficiency have similar evolutionary patterns, which
can be basically divided into three stages. The first stage is from 2007 to
2010, which can be called the rapid growth period of network efficiency. The
period from 2010 to 2015 can be regarded as a stable period, and the
efficiency of the international pesticide trade network has not changed
significantly. After 2015, it can be seen as a period of declining network
efficiency. The efficiency of the international pesticide trade network has
declined.

Table 1: The nine node metrics used in this work. Variable |  Definition |  Description |  References
---|---|---|---
In-degree |  $\sum_{i=1}^{N}a_{ij}$ |  $N$ is the number of nodes in the trade network. $a_{ij}$ is the element of adjacent matrix of trade network. |  [16]
Out-degree |  $\sum_{j=1}^{N}a_{ij}$ |  It is the count of importing partners of economy $i$. |  [16]
In-closeness |  $\frac{N-1}{\sum_{j\in{{\mathbf{n}}_{i}^{\mathrm{in}}}}l_{ij}}$ |  It is defined as the reciprocal of the mean distance of $i$ to all other nodes. |  [16]
Out-closeness |  $\frac{N-1}{\sum_{j\in{{\mathbf{n}}_{i}^{\mathrm{out}}}}l_{ij}}$ |  It is defined as the reciprocal of the mean distance of all other nodes to $i$. |  [16]
Authorities/Hubs |  - |  The authority score and hub score of a node can be calculated in an iterative way. |  [28]
PageRank |  - |  It is an important variant of eigenvector centrality. |  [29]
Betweenness |  $\sum_{st}\frac{n_{st}^{i}}{g_{st}}$ |  It is the ratio of the number of shortest paths go through the investigated node to the number of shortest paths between all pairs of nodes in the network. |  [30]
Clustering |  $\frac{2n_{i}}{k_{i}(k_{i}-1)}$ |  It characterizes the connectance of its trade partners. $n_{i}$ is the number of triangles adjacent to $i$. |  [31]

##  4 Network robustness against shocks to economies

###  4.1 Node metrics

There are many node metrics capturing the local or global characteristics of
nodes that can be used to measure the importance of nodes in different aspects
[32, 33]. We consider here nine node metrics, including clustering
coefficient, betweenness, PageRank, in-degree, out-degree, in-closeness, out-
closeness, authorities, and hubs, as shown in Table 1. These node metrics have
different traits. The clustering coefficient is defined for undirected
networks, while the other eight metrics are extracted from directed networks.
Concerning the other eight node metrics, betweenness does not distinguish
between importing and exporting economies; in-degree, PageRank, authority, and
in-closeness are calculated from exporting (thus the node under consideration
is a target node); and out-degree, out-closeness, and hub are obtained from
importing economies (thus the node under consideration is a source node). In
this section, we quantify the importance of the economy based on these nine
node metrics and further investigate their mutual correlations.

###  4.2 Network efficiency under shocks to economies

Complex network structures and network functions are closely related, and
different network structures have different roles and functions for specific
problems and backgrounds. In the analysis of network stability, indicators
such as in-degree, authorities, and betweenness have a great impact on the
network structure, but not all indicators are equally important for all
network functions. Therefore, we analyzed the differences and impacts of
important indicators in the study of network efficiency. Network efficiency
describes the transmission efficiency of material flow, information flow, or
energy flow in a complex network, which has important research and application
value.

To compare the changes in network efficiency after removing network nodes
according to the indicators of different nodes, we take the network efficiency
of the original network as the benchmark, recalculate the efficiency of the
residual network after removing network nodes, and calculate the ratio between
the two so that we can horizontally compare the changes in network efficiency
of different pesticide trade networks. The formula of the ratio $\beta^{I}$ is

$\beta^{I}(p)=\frac{E_{p}}{E},$ (4)

where $E$ represents the network efficiency of the original network, and
$E_{p}$ represents the network efficiency after removing the nodes with a
proportion of $p$ according to the given indicators.

We consider three node removal rules. Based on the node metric index $I$, we
remove the economies from the pesticide trade network and recalculate the
network efficiency for the remaining networks. We then compare the results of
the three node removal strategies. The “descend” removal rule is to remove a
proportion of $p$ of the largest trading economies based on indicator $I$ and
recalculate the network efficiency $E_{p}$ of the remaining network. The
“ascend” removal rule is to remove a proportion of $p$ of the economies with
the smallest $I$. The “random” removal method is to randomly remove a
proportion of $p$ of the trading economies and then calculate the network
efficiency $E_{p}$ of the remaining network. In the process of removing
network nodes, in addition to randomly removing nodes, node removal is mostly
performed in descending order of indicators. We introduce the removal of
indicators in ascending order, the main purpose of which is to compare with
the other two removal methods. Not all indicators meet the rule that the
higher the value, the greater the impact on network efficiency. Considering
both ascending and descending removal methods, it is not easy to miss some key
indicators. For example, the presence of weak ties during edge removal plays a
key role in network connectivity.

Figure 3: Analysis of network efficiency under shocks to the economies for the
aggregated iPTN in 2018. The lines are the $\beta^{I}$ curves calculated from
the largest connected subnetworks after removing nodes based on indicator $I$.
The words “descend”, “ascend”, and “random” in the legend correspond to the
cases where the nodes with largest $I$, smallest $I$, and random $I$ are
removed preferentially. For the random node removal stategy, we repeated 20
runs and presented the averages.

For the aggregated iPTN in 2018, we illustrate in Fig. 3 the $\beta^{I}(p)$
curves for the nine indicators $I$. The curves associated with “descend”,
“ascend”, and “random” in the figure correspond to the priorities of removing
nodes with the largest $I$ values, the smallest $I$ values, and random $I$
values, respectively. We find that the results are similar for different node
metrics, except for the clustering coefficient. For the clustering
coefficient, the two $\beta^{I}$ curves corresponding to the descending and
ascending node removal strategies almost overlap, while the $\beta^{I}$ curve
corresponding to the random node removal strategy decreases faster than the
other two curves. For the other eight node metrics, the corresponding
$\beta^{I}$ curves have similar patterns. Specifically, the curves obtained
from the descending node removal strategy decrease the fastest, while those
from the ascending node removal strategy decrease the lowest. For the
descending node removal strategy, the $\beta^{I}(p)$ value is almost nil when
the economies are affected by around 40%. In comparison, the giant component
still contains more than 40% nodes [4].

Figure 4: Network efficiency based robustness of the aggregated iPTN (a), the
iPTN of insecticides (b), the iPTN of fungicides (c), the iPTN of herbicides
(d), the iPTN of disinfectants (e), and the iPTN of rodenticides and other
similar products (f) in 2018. The words “descend”, “ascend”, and “random” in
the legend correspond to the cases where the nodes with largest $I$, smallest
$I$, and random $I$ are removed preferentially. For the random node removal
strategy, we repeated 20 runs and presented the averages.

###  4.3 Network-efficiency-based robustness

To compare the differences in network efficiency under different node removal
strategies in a more detailed way, a more quantitative method is also used to
analyze the network efficiency. There are many similar studies on the
robustness and fragility of complex networks, such as Ref. [8]. The approach
in Ref. [8] is based on percolation theory and analyzes the general condition
for the critical fraction of nodes. Our approach is based on efficiency and
robustness as presented in Ref. [34]. The area under a curve in Fig. 3
represents the impact of the corresponding indicator on the network
efficiency. The more the area is, the more obvious the network efficiency will
be. This efficiency-based robustness can be expressed as follows [34],

$R_{\beta}^{I}=\frac{1}{n}\sum_{k=0}^{n}\beta_{\rm{node}}^{I}\left(p_{k}\right),$
(5)

where $p_{k}=k/n$ indicates the proportion of economies in the pesticide trade
network to be deleted.

Figure 4 shows the extent to which the efficiency of the international
pesticide trade network is affected under the impact of shocks to the
economies, where Fig. 4(a) shows the combined trade network of all five global
trade networks of pesticide products, and Fig. 4 (b-f) correspond to the trade
networks of five pesticide products. Different markers correspond to different
node removal strategies.

As can be seen from Fig. 4, the strategies that remove economies from the
iPTNs in descending order according to the clustering coefficient show
different behaviors. For the six networks, we have

$R_{\beta,\mathrm{ascend}}^{\mathrm{Clustering}}>R_{\beta,\mathrm{descend}}^{\mathrm{Clustering}}>R_{\beta,\mathrm{random}}^{\mathrm{Clustering}},$
(6)

showing that the random node removal strategy has the largest impact on
network efficiency. The distinct behavior is mainly reflected in the
descending node removal strategy. In contrast, for the other eight node
metrics, we have

$R_{\beta,\mathrm{ascend}}^{I}>R_{\beta,\mathrm{random}}^{I}>R_{\beta,\mathrm{descend}}^{I},$
(7)

showing that the descending node removal strategy has the largest impact on
network efficiency and the ascending node removal strategy has the least
impact on network efficiency.

We also observe that the strategies based node metrics associated with import
(authority, in-closeness, and in-degree) have a larger impact on network
efficiency than those with export (hub, out-closeness, and out-degree). The
importing economies have a greater impact on network efficiency. This
conclusion is consistent with many studies, such as the study of oil network
efficiency and robustness, in which the impact of removing importing economies
is greater than that of exporting economies [13, 16].

Figure 5: Evolution of the network-efficiency-based robustness of the
aggregated iPTN and the iPTNs of insecticides, fungicides, herbicides,
disinfectants, and rodenticides and other similar products for the descending
node removal strategies based on the nine node metrics.

###  4.4 Evolution of network-efficiency-based robustness

We now turn to investigate the evolution of the network-efficiency-based
robustness of the aggregated iPTN and the five iPTNs of insecticides,
fungicides, herbicides, disinfectants, and rodenticides and other similar
products from 2007 to 2018 for the descending node removal strategies based on
the nine node metrics. The results are presented in Fig. 5. Each plot
corresponds to a node metric. Based on this metric, we remove the nodes and
then analyze the changes in network efficiency. By comparing the robustness of
different networks, we can find that the aggregated network is the most
robust. From the definition of network efficiency, it is known that the higher
the density of connected edges of the network, the smaller the average path
length between nodes, and the higher the network efficiency. Removing the
edges of the network has less impact on the connectivity and the average
shortest path for a network with a high density of edges, so the corresponding
network robustness is higher.

Once more, we can see that while the results for the other eight node metrics
are qualitatively more similar, the results for the clustering coefficient are
noticeably different. The six robustness curves exhibit an overall decreasing
trend. More precisely, the robustness curves decreased in the early years and
remained relatively stable with local fluctuations.

In contrast, for the other eight node metrics, the robustness curves show an
overall upward trend, which increased in the early years and remained stable.
Among the six networks, the aggregated iPTN is the most robust. Among the
other five networks, the iPTN of insecticides is the most robust. We also
observed a sharp decrease in robustness in 2015, especially for the two iPTNs
of herbicides and fungicides. We contend that international pesticide trade
networks are becoming more robust in terms of network efficiency in the face
of shocks to economies.

##  5 Summary

We have investigated the efficiency of the international pesticide networks of
insecticides (380891), fungicides (380892), herbicides (380893), disinfectants
(380894), and rodenticides and other similar products (380899) from 2007 to
2018, as well as the corresponding aggregated networks of all the five
categories of pesticide. We found that the network efficiency increased in the
first four years and decreased in the last four years. In addition, for each
year, the aggregated iPTN had the highest efficiency and the iPTN of
insecticides had the second highest efficiency. These observations are by and
large consistent with the time-varying pattern of the number of links [1].
There are, of course, other factors to explore.

We further investigated the robustness of iPTN efficiency, or the efficiency-
based robustness of the iPTNs, as adopted in Refs. [13, 3]. To simulate
different types of shocks to economies, we utilized three strategies by
removing nodes with descending, random, and ascending orders of nine node
metrics. We found that the efficiency-based robustness of the international
pesticide trade networks increased for all the node metrics except the
clustering coefficient. Moreover, the international pesticide trade networks
are more vulnerable when shocks hit import-oriented economies than export-
oriented economies. We also found that the aggregated iPTN is the most robust
against shocks, and the iPTN of insecticides is the second most robust. The
robustness curves are mostly $\cup$-shaped, decorated by a sharp drop in 2015
in some robustness curves, which was caused by an 8 percent plunge in market
sales in 2015.

## Acknowledgements

This work was supported by the National Natural Science Foundation of China
(72171083), the Shanghai Outstanding Academic Leaders Plan, and the
Fundamental Research Funds for the Central Universities.

## Author contributions

Funding acquisition: W-XZ; investigation: J-AL, LW, W-JX and W-XZ;
methodology: W-JX and W-XZ; supervision: W-JX and W-XZ; writing—original
draft, J-AL and W-JX; writing—review and editing: W-XZ.

Data availability statement This manuscript has no associated data or the data
will not be deposited. [Authors’ comment: The associated data in this
manuscript can be retrieved from the UN Comtrade database at
https://comtrade.un.org.]

## References

  * [1] J.A. Li, W.J. Xie, W.X. Zhou, Front. in Phys. 9, 681788 (2021).  DOI 10.3389/fphy.2021.681788
  * [2] J.A. Li, L. Wang, W.J. Xie, W.X. Zhou, Evolving community structure in the international pesticide trade networks (2023).  ArXiv:2207.08344
  * [3] W.J. Xie, J.A. Li, N. Wei, L. Wang, W.X. Zhou, Sci. Rep. 12, 19641 (2022).  DOI 10.1038/s41598-022-21777-1
  * [4] J.A. Li, L. Wang, W.J. Xie, W.X. Zhou, Journal of Management Science and Engineering p. submitted (2022)
  * [5] M. Oehlers, B. Fabian, Mathematics 9(8), 895 (2021).  DOI 10.3390/math9080895
  * [6] S.E. Schaeffer, V. Valdes, J. Figols, I. Bachmann, F. Morales, J. Bustos-Jimenez, J. Complex Netw. 9(2), cnab018 (2021).  DOI 10.1093/comnet/cnab018
  * [7] R. Albert, H. Jeong, A.L. Barabási, Nature 406(6794), 378 (2000).  DOI 10.1038/35019019
  * [8] R. Cohen, K. Erez, D. ben Avraham, S. Havlin, Phys. Rev. Lett. 85(21), 4626 (2000).  DOI 10.1103/PhysRevLett.85.4626
  * [9] D. Callaway, M. Newman, S. Strogatz, D. Watts, Phys. Rev. Lett. 85(25), 5468 (2000).  DOI 10.1103/PhysRevLett.85.5468
  * [10] R. Cohen, K. Erez, D. ben Avraham, S. Havlin, Phys. Rev. Lett. 86(16), 3682 (2001).  DOI 10.1103/PhysRevLett.86.3682
  * [11] N.J. Foti, S. Pauls, D.N. Rockmore, J. Econ. Dyn. Control 37(9), 1889 (2013).  DOI 10.1016/j.jedc.2013.04.009
  * [12] K.R. Fair, C.T. Bauch, M. Anand, Sci. Rep. 7, 7177 (2017).  DOI 10.1038/s41598-017-07202-y
  * [13] W.J. Xie, N. Wei, W.X. Zhou, J. Stat. Mech. 10, 103401 (2021).  DOI 10.1088/1742-5468/ac21da
  * [14] N. Wei, W.J. Xie, W.X. Zhou, J. Complex Netw. 10(1), cnab053 (2022).  DOI 10.1093/comnet/cnab053
  * [15] S. Chen, Y. Ding, Y. Zhang, M. Zhang, R. Nie, Energy 239, 122139 (2022).  DOI 10.1016/j.energy.2021.122139
  * [16] N. Wei, W.J. Xie, W.X. Zhou, Energy 251, 123939 (2022).  DOI 10.1016/j.energy.2022.123939
  * [17] R. Albert, I. Albert, G. Nakarado, Phys. Rev. E 69(2), 025103 (2004).  DOI 10.1103/PhysRevE.69.025103
  * [18] V. Latora, M. Marchiori, Phys. Rev. E 71(1), 015103 (2005).  DOI 10.1103/PhysRevE.71.015103
  * [19] S. Arianos, E. Bompard, A. Carbone, F. Xue, Chaos 19(1), 013119 (2009).  DOI 10.1063/1.3077229
  * [20] G.A. Pagani, M. Aiello, Physica A 392(11), 2688 (2013).  DOI 10.1016/j.physa.2013.01.023
  * [21] X. Wang, S. Guan, C.H. Lai, New J. Phys. 11, 033006 (2009).  DOI 10.1088/1367-2630/11/3/033006
  * [22] M. Zanin, F. Lillo, Eur. Phys. J.-Spec. Top. 215(1), 5 (2013).  DOI 10.1140/epjst/e2013-01711-9
  * [23] K. Diao, C. Sweetapple, R. Farmani, G. Fu, S. Ward, D. Butler, Water Res. 106, 383 (2016).  DOI 10.1016/j.watres.2016.10.011
  * [24] S. Jo, L. Gao, F. Liu, M. Li, Z. Shen, L. Xu, Z.Y. Gao, Int. J. Mod. Phys. C 32(08), 2150103 (2021).  DOI 10.1142/S0129183121501035
  * [25] F. Musciotto, S. Micciche, EPJ Data Sci. 11(1), 11 (2022).  DOI 10.1140/epjds/s13688-022-00323-z
  * [26] V. Latora, M. Marchiori, Phys. Rev. Lett. 87(19), 198701 (2001).  DOI 10.1103/PhysRevLett.87.198701
  * [27] P. Holme, B. Kim, C. Yoon, S. Han, Phys. Rev. E 65(5), 056109 (2002).  DOI 10.1103/PhysRevE.65.056109
  * [28] J. Kleinberg, J. ACM 46(5), 604 (1999).  DOI 10.1145/324133.324140
  * [29] S. Brin, L. Page, Comput. Netw. ISDN Syst. 30(1-7), 107 (1998).  DOI 10.1016/S0169-7552(98)00110-X
  * [30] L.C. Freeman, Sociometry 40(1), 35 (1977).  DOI 10.2307/3033543
  * [31] D.J. Watts, S.H. Strogatz, Nature 393, 440 (1998).  DOI 10.1038/30918
  * [32] D. Chen, L. Lu, M.S. Shang, Y.C. Zhang, T. Zhou, Physica A 391(4), 1777 (2012).  DOI 10.1016/j.physa.2011.09.017
  * [33] L. Lu, D. Chen, X.L. Ren, Q.M. Zhang, Y.C. Zhang, T. Zhou, Phys. Rep. 650, 1 (2016).  DOI 10.1016/j.physrep.2016.06.007
  * [34] C.M. Schneider, A.A. Moreira, J.S. Andrade, S. Havlin, H.J. Herrmann, Proc. Natl. Acad. Sci. U. S. A. 108(10), 3838 (2011).  DOI 10.1073/pnas.1009440108


# Learnt dynamics generalizes across tasks, datasets, and populations

U. Mahmood1, M. M. Rahman1, A. Fedorov2 , Z. Fu1, V. D. Calhoun1, 2, 3, S. M.
Plis1
Tri-institutional Center for Translational Research in Neuroimaging and Data
Science:
1Georgia State University, 2Georgia Institute of Technology, 3Emory University
Atlanta, GA, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Differentiating multivariate dynamic signals is a difficult learning problem
as the feature space may be large yet often only a few training examples are
available. Traditional approaches to this problem either proceed from
handcrafted features or require large datasets to combat the $m\gg n$ problem.
In this paper, we show that the source of the problem—signal dynamics—can be
used to our advantage and noticeably improve classification performance on a
range of discrimination tasks when training data is scarce. We demonstrate
that self-supervised pre-training guided by signal dynamics produces embedding
that generalizes across tasks, datasets, data collection sites, and data
distributions. We perform an extensive evaluation of this approach on a range
of tasks including simulated data, keyword detection problem, and a range of
functional neuroimaging data, where we show that a single embedding learnt on
healthy subjects generalizes across a number of disorders, age groups, and
datasets.

##  1 Introduction

Many interesting domains and systems are best characterized by their
mutivariate dynamics. Short and long term weather fluctuations, speech and
ambient sound, and, importantly, function of human brain to name a few. Often
in these applications we need to differentiate one condition from another
based on these dynamics. Although handcrafting features for a classifier is a
common and widespread approach, discarding potentially valuable information,
when we know that dynamics is the key to the system’s behavior, is suboptimal.
Especially when one of the goals is to learn what determines a condition of
the system.

Consider a simple problem of training devices to recognize custom keywords
with minimal training. This is a setting where only a few training samples can
be expected. Arguably more important domain to deal with is the area of mental
health.

Mental disorders manifest in behavior that is driven by disruptions in brain
dynamics [1, 2]. Functional MRI captures the nuances of spatio-temporal
dynamics that could potentially provide clues to the causes of mental
disorders and enable early diagnosis. However, the obtained data for a single
subject is of high dimensionality $m$ and to be useful for learning, and
statistical analysis, one needs to collect datasets with a large number of
subjects $n$. Yet, for any kind of a disorder, demographics or other types of
conditions, a single study is rarely able to amass datasets large enough to go
out of the $m\gg n$ mode. Traditionally this is approached by handcrafting
features [3] of much smaller dimension, effectively reducing $m$ via
dimensionality reduction. Often, the dynamics of brain function in these
representations vanishes into proxy features such as correlation matrices of
functional network connectivity (FNC) [4]. Efforts that pull together data
from various studies and increase $n$ do exist, but it is difficult to
generalize to study of smaller and more specific disease populations that
cannot be shared to become a part of these pools or are too different from the
data in them.

Our goal is to enable the direct study of systems dynamics in smaller
datasets. In the case of brain data it, in turn, can enable an analysis of
brain function. In this paper, we show how one can achieve significant
improvement in classification directly from dynamical data on small datasets
by taking advantage of publicly available large but unrelated datasets. We
demonstrate that it is possible to train a model in a self-supervised manner
on dynamics of healthy control subjects from the Human Connectome Project
(HCP) [5] and apply the pre-trained encoder to a completely different data
collected across multiple sites from healthy controls and patients. We show
that pre-training on dynamics allows the encoder to generalize across a number
of datasets and a wide range of disorders: schizophrenia, autism, and
Alzheimer’s disease. Importantly, we show that learnt dynamics generalizes
across different data distributions, as our model pre-trained on healthy
adults shows improvements in children and elderly. The generality of the
approach is also demonstrated in an application to the keyword detection
problem.

##  2 Related Work

Unsupervised pre-training is a well-known technique to get a head start for
the deep neural network. It may be considered as a regularizer which compares
to classical regularizers (i.e. L1/L2) may not vanish even with more data and
could find a robust local minima for better generalization [6]. Classical
methods are Deep Beliefs Networks (DBMs) [7] and stacked denoising
autoencoders (SDAE) [8]. Unsupervised pre-training has broad implications in
fields such as computer vision [9], natural language processing (NLP) (GPT2
[10], BERT [11], Word2Vec [12]) and automatic speech recognition (ASR) (with
SDAE [13], with DBN-HMMs [14]). However, this unsupervised pre-training is
considered to be less popular in fields other than NLP [15]. Specifically, in
computer vision, researchers usually use the model which is pre-trained in
supervised fashion on Imagenet as a starting point for downstream tasks.
Furthermore, given enough data and technical strategies, it is possible to
achieve better results on COCO object detection without supervised pre-
training on Imagenet [16].

Recent advances in unsupervised learning using self-supervised methods with
mutual information objectives have reduced the gap between supervised and
unsupervised learning on standard computer vision classification datasets [17,
18, 19, 20, 21] and scaled pre-training to very deep convolutional networks
(e.g., 50-layer ResNet). Furthermore, it influences the neuroimaging field for
classification of progression to Alzheimer’s disease from sMRI [22], learning
useful representation of the states from the frames in Atari games [23] and
also from the speech chunks for speaker identification [24]. Specifically,
authors [19] have shown that contrastive based self-supervised pre-training
can outperform supervised methods by a large margin in case of small data
(e.g., $13$ images per class in ImageNet [25]).

In most cases, due to practical reasons, researchers in brain imaging are
constrained to work with small datasets. In addition, earlier work [26, 27] in
brain imaging have been based on unsupervised methods to learn the dynamics
and structure of the brain while supervised approaches are used to perform
predictions at individual level. Such unsupervised methods include models as
linear ICA [28], HMM framework [29]. Moreover, some other nonlinear approaches
are also proposed to capture the dynamics. Examples include using Restricted
Boltzman Machines (RBMs) [30], RNN modification of ICA [31], and
reconstructions by recurrent U-Net architecture [32]. In some cases, where
dataset is very small, transfer learning as observed in some neuroimaging
application [33, 34, 35] is considered as a way to enable learning from data
and thus improve results in downstream classification. To achieve improved
performance, another idea is the data generating approach [36] which uses
synthetic data generator for pre-training, relieving the scarcity of data.

##  3 Methods

Our method is two fold. We first pre-train our encoder on large unrelated
dataset to learn improved representation of the latent factors, and then use
the pre-trained encoder for downstream task. We explain both steps in the
following sections.

###  3.1 Pre-training

####  3.1.1 Spatio-Temporal DeepInfoMax

For self-supervised pre-training we use Spatio-Temporal DeepInfoMax [23] to
maximize predictability between current latent state and future spatial state
and between consecutive spatial states (for example, on encoded time points of
the resting-state fMRI (rsfMRI)).

Let $D=\\{(u_{t},v_{s}):1\leq t,s\leq N,t\neq s\\}$ be a dataset of pairs of
values at time point $t$ and $s$ sampled from sequence with length $N$. Then
$D^{+}=\\{(u_{t},v_{s}):1\leq t\leq N-1,s=t+1\\}$ is called a dataset of
positive pairs and $D^{-}=\\{(u_{t},v_{s}):1\leq t,s\leq N,s\neq t+1\\}$ — of
negative pairs. The dataset $D^{+}$ refers to a joint distribution and $D^{-}$
— a marginal distribution. Eventually, the lower bound with InfoNCE estimator
[37] $\mathcal{I}_{f}(D^{+})$ is defined as:

$\mathcal{I}(D^{+})\geq\mathcal{I}_{f}(D^{+})\triangleq\sum^{N}_{t=1}\log\frac{\exp
f((u_{t},v_{t+1})^{+})}{\sum^{N}_{s=1}\exp f((u_{t},v_{s})^{-})},$ (1)

where $f$ is a critic function [38]. Specifically, we are using separable
critic $f(u_{t},v_{s})=\phi(u_{t})^{\intercal}\psi(v_{s})$, where $\phi$ and
$\psi$ are some embedding functions parameterized by neural networks. Such
embedding functions are used to calculate value of a critic function in same
dimensional space from two dimensional inputs. Critic learns an embedding
function such that critic assigns higher values for positive pairs compared to
negative pairs: $f(D^{+})\gg f(D^{-})$.

We define a latent state as an output $z_{t}$ of encoder $E$ and a spatial
state $c^{l}_{t}$ as the output of $l$th layer of the encoder for input
$x_{t}$ at time point $t$. To optimize the objective between current latent
state and future spatial state, the critic function for input pair
$(x_{t},x_{s})$ is $f_{global}=\phi(z_{t})^{\intercal}\psi(c^{l}_{s})$ and
between consecutive spatial states —
$f_{local}=\psi(c^{l}_{t})^{\intercal}\psi(c^{l}_{s})$. Finally, the loss is
the sum of the InfoNCE with $f_{global}$ and InfoNCE with $f_{local}$ as
$L={I}_{f_{global}}+{I}_{f_{local}}$. The scheme of the ST-DIM is shown in
Figure 1.

Figure 1: Spatio-Temporal Deep InfoMax (ST-DIM) and autoencoder architecture
used in pre-training. To enable comparison, the same encoder architecture is
used for both ST-DIM and autoencoder models.

####  3.1.2 Autoencoder

An autoencoder [39] is often viewed as an unsupervised learning representation
technique useful for obtaining an informative representation of the input for
downstream tasks [15]. As implied in its name, it consists of two phases,
namely, reduction (encoding) and reconstruction (decoding). In the reduction
phase, the autoencoder compresses its input preserving useful latent
representation of the input. In the reconstruction phase, it learns to create
outputs similar to inputs based on the latent representation achieved in
reduction phase. Both phases are trained together to obtain optimal results.
The reduction and reconstruction phases can be considered as two mathematical
functions $\Phi:\mathcal{X}\mapsto\mathcal{F}$ and
$\Psi:\mathcal{F}\mapsto\mathcal{X}$ respectively. We use the mean-squared
error for the loss function which is defined as:
$\mathcal{L}(x,\hat{x})=\|x-\hat{x}\|^{2}=\|x-\Psi(\Phi(x))\|^{2}$, where
$x,\hat{x}\in\mathcal{X}$ and refer to input and output of the autoencoder
respectively. Though autoencoder can be implemented using different types of
networks, we use the same network architecture as used for ST-DIM pre-training
(See Section 4.1) to make them ideally comparable. The architecture used in
pre-training is shown in Figure 1.

###  3.2 Transfer and Supervised Learning

In the downstream task, we use the representation of the encoder pre-trained
using ST-DIM or autoencoder as input to a biLSTM [40] and add a simple binary
classifier on top. We also use the same encoder architecture to learn end-to-
end supervised model (i.e., no pre-trained weights) for the same
classification tasks. The overall architecture used for downstream tasks is
shown in Figure 2. Refer to section 4.1 for further details.

Figure 2: Architecture used in downstream task. Windows (time slices) are
encoded using encoder described in Section 3.1.2 and passed through a biLSTM.
Outputs of the last forward and backward LSTM layers are concatenated and
passed through a dense layer for classification.

##  4 Experiments

In this section we study the performance of our model on both, synthetic and
real data. To compare and show the advantage of pre-training on large
unrelated dataset we use three different kind of models — 1) FPT (Frozen Pre-
Trained): The pre-trained encoder is not further trained on the dataset of
downstream task, 2) UFPT (Unfrozen Pre-Trained): The pre-trained encoder is
further trained on the dataset of downstream task and 3) NPT (Not Pre-
trained): The encoder is not pre-trained at all and only trained on the small
dataset of downstream task. The models are shown in Figure 3. In each
experiment, we compare all three models to demonstrate the effectiveness of
unsupervised pre-training.

Figure 3: Three different models are used for downstream tasks: 1) In FPT, we
pre-train our encoder but keep it frozen during downstream task. 2) In UFPT,
we pre-train our encoder but continue training during downstream task. 3) In
NPT, we use the exactly same architecture but with no pre-trained weights for
the encoder.

###  4.1 Setup

Encoder for simulation experiment consists of $4$ $1$D convolutional layers
with output features $(32,64,128,64)$, kernel sizes $(4,4,3,2)$ respectively,
followed by ReLU [41] after each layer followed by a linear layer with $256$
units. For real data experiments, we use $3$ $1$D convolutional layers with
output features $(64,128,200)$, kernel sizes $(4,4,3)$ respectively, followed
by ReLU after each layer followed by a linear layer with $256$ units. We use
stride $1$ for all of the convolution layers. For the autoencoder based pre-
training for simulation experiment, we use the same encoder as for ST-DIM in
the reduction phase. For the decoder, we use the reverse architecture of the
encoder. Precisely, we use a linear layer with 704 units followed by $4$ $1$D
transpose convolution layers with output features $(128,64,32,10)$, kernel
sizes $(2,3,4,4)$ respectively and all with stride $1$, that result in
$10\times 20$ windows at the output.

In ST-DIM based pre-training, for all possible pairs in the batch, we take
features $c^{3}$ and $z$ after the $3$rd convolutional layer and output layer
respectively. We embed both $c^{3}$ and $z$ using $\psi$ and $\phi$
respectively to a $128$ dimensional vector to compute the score of a critic
function $f_{global}$ or $f_{local}$. Using these scores, we compute the loss.
The neural networks are trained using Adam optimizer [42]. The weights for
encoder are initialized using orthogonal initialization. However, the LSTM and
the final fully connected layer are initialized using Xavier [43].

As we are more interested in subjects for classification task, we feed each
time series (ICA time courses) into the encoder in the form of a sequence of
windows. The encoder encodes the windows of input data into latent
representation. The latent representation of the entire time series is then
concatenated and passed to a biLSTM with hidden dimension of size $200$. The
output of biLSTM is then used as input to a feed forward network of two linear
layers with $200$ and $2$ units to perform binary classification.

For experiments on the downstream tasks, a hold out is selected for testing
and is never used through the training/validation phase. For each downstream
task, the number of subjects used for supervised training is gradually
increased within a range to observe the effectiveness of pre-training in
downstream task with varied number of training subjects . For each experiment,
10 trials are performed to ensure random selection of training subjects and,
in each case, the performance is evaluated on the hold out dataset (test
data).

For brain data, each of the models (FPT, UFPT, NPT) yields the best results
based on its respective gain value of Xavier [43] initialization used for
biLSTM. To find the best gain value for each model, $20$ values in the range
$0-1$ are tried with an increment of $0.05$. For each value, $10$ experiments
are performed and best value is chosen based on the results on validation
dataset. Refer to Table 1 for more parametric details of the models.

| Simulations | Real Data
---|---|---
Training Batch Size | $16$ | min(32, no. of subjects)
Validation Batch Size | 200 | no. of subjects
Test Batch Size | 200 | no. of subjects
Initial Learning Rate | 3e-4 | 3e-4
Learning Rate Scheduler | None | Reduce LR on Plateau
Max Epochs | 500 | 3000

Table 1: Parameter details of models used for experiments

###  4.2 Simulation

To simulate the data, we generate multiple $10$-node graphs with $10\times 10$
stable transition matrices. Using these we generate multivariate time series
with autoregressive (VAR) and structural vector autoregressive (SVAR) models
[44].

First, we generate $50$ VAR times series with size $10\times 20000$. Then we
split our dataset to $50\times 10\times 14000$ samples for training, $50\times
10\times 4000$ —for validation and $50\times 10\times 2000$ — for testing.
Using these samples, We pre-train an encoder to learn consecutive windows
(positive examples) from the same VAR time series. As mentioned in Section
3.1.2, we also use autoencoder for pre-training the same encoder and show the
effectiveness of ST-DIM to learn time-series dynamics in self-supervised
manner. After pre-training, we use our pre-trained encoder for complete-time
series classification.

In the final downstream task, we classify the whole time-series whether it is
generated by VAR or SVAR (undersampled VAR at rate 2). We create $400$ graphs
with corresponding stable transition matrices and generate $2000\times
10\times 4000$ samples ($5$ for each) and split as $1600\times 10\times 4000$
for training, $200\times 10\times 4000$ for validation and $200\times 10\times
4000$ for hold-out test. Here we also use $10\times 20$ windows as a single
time-point input.

First, we train our encoder to learn consecutive $10\times 20$ windows from
the VAR time series using InfoNCE based loss, and secondly, we train a
supervised classifier based on windows. This window-based classification
provides promising results. However, in solving similar real problems, we are
more interested in subjects, i.e., entire time series, rather than a single-
window for classification. Hence, we perform classification based on the whole
time-series. In this setting, the entire time-series is encoded as a sequence
of representations and fed through a biLSTM. Then, as described in Section
4.1, using the feed forward network on top of the last hidden state of the
biLSTM, we map the representation to classification scores.

We observe that the ST-DIM based pre-trained models can easily be fine-tuned
only with small amount of downstream data. In this classification task, our
model can classify a randomly chosen time-series as a sample of VAR or SVAR.
Note, with very few samples, models based on the pre-trained encoder (FPT and
UFPT) outperform supervised models. However, as the number of samples grows,
the accuracy achieved with or without pre-training levels out. We also notice
that autoencoder based self-supervised pre-training does not assist in VAR vs.
SVAR classification. Consequently, we use only ST-DIM based pre-training for
all the real data experiments. Refer to Figure 4 for the results of simulation
experiments.

Figure 4: Area Under Curve (AUC) scores for VAR vs. SVAR time-series
classification using ST-DIM and autoencoder based pre-training methods. ST-DIM
based pre-training greatly improves the performance of downstream task with
small datasets. On the other side, autoencoder based pre-training fails to
learn dynamics and thus exhibits poor performance.

###  4.3 Keyword Detection

To show the broad implications of unsupervised pre-training, we first apply it
to a simple problem of keyword detection in audio files. We choose this
problem as it has many practical applications (e.g., virtual assistants in
smart phones, robots). We use LibriSpeech ASR corpus [45] for pre-training and
Speech Commands Dataset [46] for supervised training. The audio files of both
datasets are combined with a background noise of coffee shop collected from
[47] to make pre-training and classification harder.

Figure 5: Left: For pre-training, we combine audio files from LibriSpeech
dataset with background noise of coffee shop. T is the length of audio files
which ranges from $1$ to $20$ seconds. Right: For training, we superimpose a
T/2 length audio of word "cat" padded with T/2 zeros onto a background noise
of coffee shop of length T ($T=2$). For both pre-training and training, we
calculate the mel-spectrogram of the combined audio files that results in a
matrix of size _components $\times$ time courses_ for each audio file.

For pre-training, we calculate mel-spectrogram for all audio files from each
subject and then concatenate together, resulting in a matrix of dimensions
$C\times N$, where N, ranging anywhere from $10000$ to $60000$ depending on
the length and number of audio files for each subject, is the time courses for
the subject and $C=128$ is the number of components in mel-spectrogram
calculated at each time point. Figure 5 left shows the pre-processing of one
audio file from a subject used for pre-training.

Out of all the available subjects, we select 416 which have large number of
time points $(N\geqslant 20k)$. We use 300 subjects for training and 116 for
validation. For pre-training, we use the algorithm as described in section
3.1.1 and achieve accuracy of $\thicksim 0.95$ on the validation dataset.

For the downstream classification task, we collect samples of audios for the
keyword "cat" from [46] which has 1515 audio files for the keyword. To create
"cat" class examples, we superimpose each of the keyword audios having length
of 1 second onto a 2 seconds long background noise at random location,
resulting into a two seconds long audio consisting of background noise and
keyword cat. Figure 5 Right shows the pre-processing used in downstream
classification for a "cat" example. To create "no cat" class examples, we use
another 1515, 2 seconds long audio files containing only background noise of
coffee shop. Thus, we use 3030 audio files for the downstream task ("cat"/"no
cat" classification). For all of the 3030 audio files, we compute mel-
spectrogram to convert each of them into a matrix of size _components $\times$
time courses_.

We perform 20 cross validations, each having a hold out of size 128 for test
subjects. For each hold out, we conduct 5 different experiments varying number
of samples (ranging from $15$ to $100$) from each class to be used for
training. For each experiment, we perform 10 trials resulting into 200 trials
to get accurate estimates of the performance.

As observed in the simulation experiments, Figure 6 shows the importance of
ST-DIM based unsupervised pre-training for keyword detection problem. Even
with very few training subjects, FPT and UFPT outperform NPT by a noticeable
margin. As the number of samples increases, the performance of NPT improves
but still remains much lower than the pre-trained counterparts.

Figure 6: AUC scores for all the models on keyword detection task. Notice that
pre-trained models provide noticeably better performance than NPT even with
just $15$ subjects for training. The pre-trained models continue to improve
significantly reaching more than $0.9$ with just 100 subjects, whereas NPT
fails to achieve even $0.8$ with the same number of training subjects.

###  4.4 Brain Imaging

####  4.4.1 Datasets

Next, we apply the same unsupervised pre-training method to brain imagining
data. For encoder pre-training, we use HCP [5] consortium dataset and apply
the pre-trained encoder for further downstream tasks. We apply the same pre-
trained encoder for three different types of diseases spanning four datasets
to classify schizophrenia, autism and Alzheimer’s diseases. We use resting
fMRI data for all brain data experiments. Refer to Figure 1 for the details of
the datasets for disease classification.

Figure 7: Datasets used for pre-training and classification tasks. Healthy
controls from the HCP [5] are used for encoder pre-training guided by data
dynamics alone111Human silhouettes are by Natasha Sinegina for Creazilla.com
without modifications, https://creativecommons.org/licenses/by/4.0/. The pre-
trained encoder is then used in downstream classification tasks of $3$
different diseases, $4$ independently collected datasets, many of which
contain data from a number of sites, and consist of populations with
significant age difference. The age distributions in the datasets have the
following means, medians and standard deviations: HCP: 29.31, 29.00, 3.667;
ABIDE: 17.04, 15.40, 7.29; COBRE: 37.96, 37, 12.90; FBIRN: 37.87, 38, 11.25;
OASIS: 67.67, 68, 8.92.

Four datasets used in this study are collected from the FBIRN (Function
Biomedical Informatics Research Network222These data were downloaded from the
Function BIRN Data Repository, Project Accession Number 2007-BDR-6UHZ1.) [48]
project, from the COBRE (Center of Biomedical Research Excellence) [49]
project, from the release 1.0 of ABIDE (Autism Brain Imaging Data
Exchange333http://fcon_1000.projects.nitrc.org/indi/abide/) [50] and from
release 3.0 of OASIS (Open Access Series of Imaging
Studies444https://www.oasis-brains.org/) [51]. Written informed consent was
obtained from all participants of each dataset under protocols approved by the
institutional review board (IRB).

####  4.4.2 Preprocessing

We preprocessed the fMRI data using statistical parametric mapping (SPM12,
http://www.fil.ion.ucl.ac.uk/spm/) under MATLAB 2016 environment. A rigid body
motion correction was performed using the toolbox in SPM to correct subject
head motion, followed by the slice-timing correction to account for timing
difference in slice acquisition. The fMRI data were subsequently warped into
the standard Montreal Neurological Institute (MNI) space using an echo planar
imaging (EPI) template and were slightly resampled to $3\times 3\times 3$ mm3
isotropic voxels. The resampled fMRI images were finally smoothed using a
Gaussian kernel with a full width at half maximum (FWHM) = $6$ mm. After the
preprocessing. We included subjects in the analysis if the subjects have head
motion $\leq 3^{\circ}$ and $\leq 3$ mm, and with functional data providing
near full brain successful normalization [52].

Figure 8: ICA time courses are computed from the resting state fMRI data.
Results contain statistically independent spatial maps (top) and their
corresponding time courses.

For each dataset, $100$ ICA components as shown in Figure 8 are acquired using
the same procedure described in [52]. However, only $53$ non-noise components
as determined per slice (time point) are used in pre-training of encoder and
on downstream task. For experiments, including both pre-training and
classification the fMRI sequence is divided into windows of $20$ time points.

####  4.4.3 Schizophrenia

For schizophrenia classification, we conduct experiments on two different
datasets, FBIRN [48] and COBRE [49]. The datasets contain labeled
Schizophrenia (SZ) and Healthy Control (HC) subjects.

FBIRN The dataset has total $311$ subjects consisting of $150$ HC and $161$
affected with SZ. Each subject has $53$ non-noise components with $140$ time
points. We use two hold-out sets with sizes $32$ and $64$ respectively for
validation and test. The remaining data are used for supervised training. With
$140$ time points, we create $53\times 20$-sized windows with $50\%$ overlap
along time dimension resulting in $13$ windows for each subject. The details
of the results are shown in Figure 9.

Figure 9: AUC scores for all the three models (Refer to Figure 3) on FBIRN
dataset. It is noticeable that even with only $15$ subjects for supervised
training, the median AUC scores of FPT and NPT differ by a large margin
$(10\%)$.

As we can see, the pre-trained models (FPT and UFPT) outperform NPT. Even with
very few subjects (15) for training, the difference between AUC scores
obtained respectively with FPT and NPT is reasonably large ($\simeq 0.1$).

COBRE The dataset has total $157$ subjects — a collection of $68$ HC and $89$
affected with SZ. Like FBIRN, each subject has $53$ non-noise components in
its ICA time courses with $140$ time points. We use two hold-out sets of size
$32$ each respectively for validation and test. The remaining data are used
for supervised training. With $140$ time points, we create $53\times 20$-sized
windows with no overlapping resulting in $7$ windows for each subject. Unlike
FBIRN, it has been impossible to increase the number of subjects for
downstream training due to insufficiency of data.

Figure 10: AUC scores for all the three models on COBRE dataset. It is obvious
that even with $15$ subjects for training, FPT outperforms NPT noticeably,
that is, the difference between two median AUC scores is remarkable ($\simeq
0.15$).

The results in Figure 10 adheres to the results of FBIRN. That is, with only
$15$ training subjects, FPT and UFPT perform significantly better than NPT
having $\simeq 0.15$ difference in their median AUC scores.

####  4.4.4 Autism

Figure 11: AUC scores for all the three models on ABIDE dataset. Like
experiments on FBIRN and COBRE, it is evident that the pre-trained models
consistently perform better than NPT.

The dataset ABIDE has total $569$ subjects, of which, $255$ are HC and $314$
are affected with autism. Like other datasets, each subject has $53$ non-noise
components with $140$ time points. We use two hold-out sets of size $100$ each
respectively for validation and test purpose. The remaining data are used for
downstream training i.e., autism vs. HC classification. Like COBRE dataset,
with $140$ time points, we create $53\times 20$-sized windows with no
overlapping resulting in $7$ windows for each subject. Refer to Figure 11 for
the details of the experimental results.

As seen in the figure, same pre-trained encoder performs reasonably better
than NPT for autism vs. HC classification and thus reinforces our hypothesis
that unsupervised pre-training learns signal dynamics useful for downstream
tasks. Note the difference between age ranges of ABIDE and HCP datasets. The
age range of ABIDE is much lower than that of HCP dataset used for pre-
training. HCP dataset contains subjects of different ages with means 30.01 and
28.48, medians 30 and 28, and standard deviations 3.522 and 3.665 years
respectively for female and male, whereas ABIDE dataset has overall mean
17.04, median 15.40 and standard deviation 7.29 years. Refer to Figure 1 for
the demographic information of all the datasets. The dissimilarity in the age
range is supposed to cause significant difference between these two datasets
as the structure of brain and thought process of children is obviously
different than adults.

####  4.4.5 Alzheimer’s disease

The dataset OASIS [51] has total $372$ subjects with equal number ($186$) of
HC and AZ patients. We use two hold-out sets each of size $64$ respectively
for validation and test purpose. The remaining are used for supervised
training. Unlike other datasets described earlier, it has only $120$ time
points though the number of non-noise componets is same ($53$) as other
datasets. With $120$ time points, we use six $53\times 20$-sized non-
overlapping windows for each subject. Refer to Figure 12 for the details of
the experiments.

As we see, the AUC scores of all the three models for 15 subjects is
$\thicksim 0.5$, which can be treated as merely random guess. However, as the
number of subjects increases, pre-trained models gradually start performing
better than NPT which, in fact, even with $120$ subjects fails to learn from
the data. We suspect that the reason why pre-trained models do not work for
$15$ subjects is that the data set is much different than HCP. The big age gap
between subjects of HCP and OASIS is a major difference and $15$ subjects are
even not enough for pre-trained models.

Figure 12: AUC scores for all the models on OASIS dataset. As we continue
increasing the number of subjects, the pre-trained models start learning and
thus improve their respective scores. However, notice that the NPT model even
with $120$ subjects didn’t significantly improve its predictability.

##  5 Conclusions and Future Work

As we have demonstrated, self-supervised pre-training of a spatiotemporal
encoder gives significant improvement on the downstream tasks in both keyword
detection and brain imaging datasets. Pre-training on fMRI of healthy subjects
provides benefits that transfer across datasets, collection sites, and
multiple disease classification with varying age gap. Learning dynamics of
fMRI helps to improve classification results for schizophrenia, autism,
Alzheimer’s dieseases and speed up the convergence of the algorithm on small
datasets, that otherwise do not provide reliable generalizations. Although the
utility of these results is highly promising by itself, we conjecture that
direct application to spatio-temporal data will warrant benefits beyond
improved classification accuracy in the future work. Working with ICA
components is a smaller and thus easier to handle space that exhibits all
dynamics of the signal, in future we will move beyond ICA pre-processing and
work with fMRI data directly. We expect model introspection to yield insight
into the spatio-temporal biomarkers of schizophrenia. In future work, we will
use the same analogously pre-trained encoder on datasets with various other
mental disorders such as MCI and bipolar. We are optimistic about the outcome
because the proposed pre-training is oblivious of the downstream use and is
done in a manner quite different from the classifier’s work. It may indeed be
learning crucial information about dynamics that might contain important clues
into the nature of mental disorders.

##  6 Acknowledgement

This study was supported by startup funds to SMP and in part by NIH grant
R01EB020407.

Data for healthy subjects was provided [in part] by the Human Connectome
Project, WU-Minn Consortium (Principal Investigators: David Van Essen and
Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that
support the NIH Blueprint for Neuroscience Research; and by the McDonnell
Center for Systems Neuroscience at Washington University.

Data for Schizophrenia classification was used in this study were downloaded
from the Function BIRN Data Repository
(http://bdr.birncommunity.org:8080/BDR/), supported by grants to the Function
BIRN (U24-RR021992) Testbed funded by the National Center for Research
Resources at the National Institutes of Health, U.S.A. and from the
COllaborative Informatics and Neuroimaging Suite Data Exchange tool (COINS;
http://coins.trendscenter.org) and data collection was performed at the Mind
Research Network, and funded by a Center of Biomedical Research Excellence
(COBRE) grant 5P20RR021938/P20GM103472 from the NIH to Dr.Vince Calhoun.

Data for Alzheimer’s was provided by OASIS-3: Principal Investigators: T.
Benzinger, D. Marcus, J. Morris; NIH P50AG00561, P30NS09857781, P01AG026276,
P01AG003991, R01AG043434, UL1TR000448, R01EB009352. AV-45 doses were provided
by Avid Radiopharmaceuticals, a wholly owned subsidiary of Eli Lilly.

Autism data was provided by ABIDE. We acknowledge primary support for the work
by Adriana Di Martino provided by the (NIMH K23MH087770) and the Leon Levy
Foundation and primary support for the work by Michael P. Milham and the INDI
team was provided by gifts from Joseph P. Healy and the Stavros Niarchos
Foundation to the Child Mind Institute, as well as by an NIMH award to MPM
(NIMH R03MH096321).

## References

  * [1] David P Goldberg and Peter Huxley.  Common mental disorders: a bio-social model. Tavistock/Routledge, 1992.
  * [2] Vince D Calhoun, Robyn Miller, Godfrey Pearlson, and Tulay Adalı.  The chronnectome: time-varying connectivity networks as the next frontier in fmri data discovery.  Neuron, 84(2):262–274, 2014.
  * [3] Ali Khazaee, Ata Ebrahimzadeh, and Abbas Babajani-Feremi.  Application of advanced machine learning methods on resting-state fMRI network for identification of mild cognitive impairment and Alzheimer’s disease.  Brain Imaging and Behavior, 10(3):799–817, Sep 2016.
  * [4] Weizheng Yan, Sergey Plis, Vince D Calhoun, Shengfeng Liu, Rongtao Jiang, Tian-Zi Jiang, and Jing Sui.  Discriminating schizophrenia from normal controls using resting state functional network connectivity: A deep neural network and layer-wise relevance propagation method.  In 2017 IEEE 27th International Workshop on Machine Learning for Signal Processing (MLSP), pages 1–6. IEEE, 2017.
  * [5] David C Van Essen, Stephen M Smith, Deanna M Barch, Timothy EJ Behrens, Essa Yacoub, Kamil Ugurbil, Wu-Minn HCP Consortium, et al.  The WU-Minn human connectome project: an overview.  Neuroimage, 80:62–79, 2013.
  * [6] Dumitru Erhan, Yoshua Bengio, Aaron Courville, Pierre-Antoine Manzagol, Pascal Vincent, and Samy Bengio.  Why does unsupervised pre-training help deep learning?  Journal of Machine Learning Research, 11(Feb):625–660, 2010.
  * [7] Geoffrey E Hinton and Russ R Salakhutdinov.  Using deep belief nets to learn covariance kernels for gaussian processes.  In J. C. Platt, D. Koller, Y. Singer, and S. T. Roweis, editors, Advances in Neural Information Processing Systems 20, pages 1249–1256. Curran Associates, Inc., 2008.
  * [8] Pascal Vincent, Hugo Larochelle, Yoshua Bengio, and Pierre-Antoine Manzagol.  Extracting and composing robust features with denoising autoencoders.  In Proceedings of the 25th international conference on Machine learning, pages 1096–1103. ACM, 2008.
  * [9] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  In Advances in neural information processing systems, pages 1097–1105, 2012.
  * [10] Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever.  Language models are unsupervised multitask learners.  OpenAI Blog, 1(8), 2019.
  * [11] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  arXiv preprint arXiv:1810.04805, 2018.
  * [12] Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean.  Distributed representations of words and phrases and their compositionality.  In Advances in neural information processing systems, pages 3111–3119, 2013.
  * [13] Jonas Gehring, Yajie Miao, Florian Metze, and Alex Waibel.  Extracting deep bottleneck features using stacked auto-encoders.  In 2013 IEEE international conference on acoustics, speech and signal processing, pages 3377–3381. IEEE, 2013.
  * [14] Dong Yu, Li Deng, and George Dahl.  Roles of pre-training and fine-tuning in context-dependent dbn-hmms for real-world speech recognition.  In Proc. NIPS Workshop on Deep Learning and Unsupervised Feature Learning, 2010.
  * [15] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep learning.  MIT press, 2016.
  * [16] Kaiming He, Ross Girshick, and Piotr Dollar.  Rethinking imagenet pre-training.  In The IEEE International Conference on Computer Vision (ICCV), October 2019.
  * [17] Aaron van den Oord, Yazhe Li, and Oriol Vinyals.  Representation learning with contrastive predictive coding.  arXiv preprint arXiv:1807.03748, 2018.
  * [18] R Devon Hjelm, Alex Fedorov, Samuel Lavoie-Marchildon, Karan Grewal, Phil Bachman, Adam Trischler, and Yoshua Bengio.  Learning deep representations by mutual information estimation and maximization.  arXiv preprint arXiv:1808.06670, 2018.
  * [19] Olivier J Hénaff, Ali Razavi, Carl Doersch, SM Eslami, and Aaron van den Oord.  Data-efficient image recognition with contrastive predictive coding.  arXiv preprint arXiv:1905.09272, 2019.
  * [20] Philip Bachman, R Devon Hjelm, and William Buchwalter.  Learning representations by maximizing mutual information across views.  arXiv preprint arXiv:1906.00910, 2019.
  * [21] Kaiming He, Haoqi Fan, Yuxin Wu, Saining Xie, and Ross Girshick.  Momentum contrast for unsupervised visual representation learning, 2019\.
  * [22] Alex Fedorov, R Devon Hjelm, Anees Abrol, Zening Fu, Yuhui Du, Sergey Plis, and Vince D Calhoun.  Prediction of progression to Alzheimers disease with deep InfoMax.  arXiv preprint arXiv:1904.10931, 2019.
  * [23] Ankesh Anand, Evan Racah, Sherjil Ozair, Yoshua Bengio, Marc-Alexandre Côté, and R Devon Hjelm.  Unsupervised state representation learning in Atari.  arXiv preprint arXiv:1906.08226, 2019.
  * [24] Mirco Ravanelli and Yoshua Bengio.  Learning speaker representations with mutual information.  arXiv preprint arXiv:1812.00271, 2018.
  * [25] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei.  ImageNet: A Large-Scale Hierarchical Image Database.  In CVPR09, 2009.
  * [26] Meenakshi Khosla, Keith Jamison, Gia H Ngo, Amy Kuceyeski, and Mert R Sabuncu.  Machine learning in resting-state fMRI analysis.  Magnetic resonance imaging, 2019.
  * [27] Sergey M Plis, Devon Hjelm, Ruslan Salakhutdinov, Elena A Allen, Henry Jeremy Bockholt, Jeffrey D Long, Hans J Johnson, Jane Paulsen, Jessica A Turner, and Vince D Calhoun.  Deep learning for neuroimaging: a validation study.  Frontiers in Neuroscience, 8(229), 2014.
  * [28] Vince D Calhoun, Tulay Adali, Godfrey D Pearlson, and JJ Pekar.  A method for making group inferences from functional MRI data using independent component analysis.  Human brain mapping, 14(3):140–151, 2001.
  * [29] Harini Eavani, Theodore D Satterthwaite, Raquel E Gur, Ruben C Gur, and Christos Davatzikos.  Unsupervised learning of functional network dynamics in resting state fmri.  In International Conference on Information Processing in Medical Imaging, pages 426–437. Springer, 2013.
  * [30] R Devon Hjelm, Vince D Calhoun, Ruslan Salakhutdinov, Elena A Allen, Tulay Adali, and Sergey M Plis.  Restricted Boltzmann machines for neuroimaging: an application in identifying intrinsic networks.  NeuroImage, 96:245–260, 2014.
  * [31] R Devon Hjelm, Eswar Damaraju, Kyunghyun Cho, Helmut Laufs, Sergey M Plis, and Vince D Calhoun.  Spatio-temporal dynamics of intrinsic networks in functional magnetic imaging data using recurrent neural networks.  Frontiers in neuroscience, 12:600, 2018.
  * [32] Meenakshi Khosla, Keith Jamison, Amy Kuceyeski, and Mert R Sabuncu.  Detecting abnormalities in resting-state dynamics: An unsupervised learning approach.  In International Workshop on Machine Learning in Medical Imaging, pages 301–309. Springer, 2019.
  * [33] Arthur Mensch, Julien Mairal, Danilo Bzdok, Bertrand Thirion, and Gaël Varoquaux.  Learning neural representations of human cognition across many fMRI studies.  In Advances in Neural Information Processing Systems, pages 5883–5893, 2017.
  * [34] Hailong Li, Nehal A. Parikh, and Lili He.  A novel transfer learning approach to enhance deep neural network classification of brain functional connectomes.  Frontiers in Neuroscience, 12:491, 2018.
  * [35] Armin W Thomas, Klaus-Robert Müller, and Wojciech Samek.  Deep transfer learning for whole-brain fMRI analyses.  arXiv preprint arXiv:1907.01953, 2019.
  * [36] Alvaro Ulloa, Sergey Plis, and Vince Calhoun.  Improving classification rate of schizophrenia using a multimodal multi-layer perceptron model with structural and functional MR.  arXiv preprint arXiv:1804.04591, 2018.
  * [37] Aaron van den Oord, Yazhe Li, and Oriol Vinyals.  Representation learning with contrastive predictive coding.  arXiv preprint arXiv:1807.03748, 2018.
  * [38] Michael Tschannen, Josip Djolonga, Paul K Rubenstein, Sylvain Gelly, and Mario Lucic.  On mutual information maximization for representation learning.  arXiv preprint arXiv:1907.13625, 2019.
  * [39] Mark A Kramer.  Nonlinear principal component analysis using autoassociative neural networks.  AIChE journal, 37(2):233–243, 1991.
  * [40] Mike Schuster and Kuldip K Paliwal.  Bidirectional recurrent neural networks.  IEEE Transactions on Signal Processing, 45(11):2673–2681, 1997\.
  * [41] Xavier Glorot, Antoine Bordes, and Yoshua Bengio.  Deep sparse rectifier neural networks.  In Proceedings of the fourteenth international conference on artificial intelligence and statistics, pages 315–323, 2011.
  * [42] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [43] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pages 249–256, 2010.
  * [44] Helmut Lütkepohl.  New introduction to multiple time series analysis.  Springer Science & Business Media, 2005.
  * [45] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur.  Librispeech: an ASR corpus based on public domain audio books.  ICASSP, pages 5206–5210, 2015.
  * [46] Pete Warden.  Speech commands: A dataset for limited-vocabulary speech recognition.  2018\.
  * [47] Duff The Psych.  One hour of hq coffee shop background noise.  2014\.
  * [48] David B Keator, Theo GM van Erp, Jessica A Turner, Gary H Glover, Bryon A Mueller, Thomas T Liu, James T Voyvodic, Jerod Rasmussen, Vince D Calhoun, Hyo Jong Lee, et al.  The function biomedical informatics research network data repository.  Neuroimage, 124:1074–1079, 2016.
  * [49] Mustafa S Çetin, Fletcher Christensen, Christopher C Abbott, Julia M Stephen, Andrew R Mayer, José M Cañive, Juan R Bustillo, Godfrey D Pearlson, and Vince D Calhoun.  Thalamus and posterior temporal lobe show greater inter-network connectivity at rest and across sensory paradigms in schizophrenia.  Neuroimage, 97:117–126, 2014.
  * [50] Adriana Di Martino, Chao-Gan Yan, Qingyang Li, Erin Denio, Francisco X Castellanos, Kaat Alaerts, Jeffrey S Anderson, Michal Assaf, Susan Y Bookheimer, Mirella Dapretto, et al.  The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism.  Molecular psychiatry, 19(6):659, 2014.
  * [51] Eugene H Rubin, Martha Storandt, J Philip Miller, Dorothy A Kinscherf, Elizabeth A Grant, John C Morris, and Leonard Berg.  A prospective study of cognitive function and onset of dementia in cognitively healthy elders.  Archives of neurology, 55(3):395–401, 1998.
  * [52] Zening Fu, Arvind Caprihan, Jiayu Chen, Yuhui Du, John C Adair, Jing Sui, Gary A Rosenberg, and Vince D Calhoun.  Altered static and dynamic functional network connectivity in alzheimer’s disease and subcortical ischemic vascular disease: shared and specific brain connectivity abnormalities.  Human Brain Mapping, 2019.