# Maintenance scheduling of manufacturing systems based on optimal price of
the network

Pegah Rokhforoz1 and  Olga Fink2 Chair of Intelligent Maintenance Systems, ETH
Zurich, Switzerland, and School of Electrical and Computer Engineering,
University of Tehran, Tehran, Iran.  Chair of Intelligent Maintenance Systems,
ETH Zurich, Switzerland. Corresponding author, email<EMAIL_ADDRESS>
###### Abstract.

Goods can exhibit positive externalities impacting decisions of customers in
socials networks. Suppliers can integrate these externalities in their pricing
strategies to increase their revenue. Besides optimizing the prize, suppliers
also have to consider their production and maintenance costs. Predictive
maintenance has the potential to reduce the maintenance costs and improve the
system availability. To address the joint optimization of pricing with network
externalities and predictive maintenance scheduling based on the condition of
the system, we propose a bi-level optimization solution based on game theory.
In the first level, the manufacturing company decides about the predictive
maintenance scheduling of the units and the price of the goods. In the second
level, the customers decide about their consumption using an optimization
approach in which the objective function depends on their consumption, the
consumption levels of other customers who are connected through the graph, and
the price of the network which is determined by the supplier. To solve the
problem, we propose the leader-multiple-followers game where the supplier as a
leader predicts the strategies of the followers. Then, customers as the
followers obtain their strategies based on the leader’s and other followers’
strategies. We demonstrate the effectiveness of our proposed method on a
simulated case study. The results demonstrate that knowledge of the social
network graph results in an increased revenue compared to the case when the
underlying social network graph is not known. Moreover, the results
demonstrate that obtaining the predictive maintenance scheduling based on the
proposed optimization approach leads to an increased profit compared to the
baseline decision-making (perform maintenance at the degradation limit).

##  1\. introduction

Social networks can have a significant influence on the economic system [1,
2]. In a social network, users connect and communicate with each other. Hence,
the decision of a user has the potential to affect the decisions of other
users, particularly for goods that exhibit externalities. In such cases, users
receive different extents of externalities from the consumption of other users
to whom they are connected in their social network. The information about the
consumption behavior of users spreads in the network. If suppliers know the
graph network and the consumption behaviour of the users, they can try to
integrate this information in their optimization. Aiming to maximize their
revenue, suppliers determine the price of the good for each user such that
their revenue is maximized.

Besides considering the prizing of the manufactured goods, manufacturers also
need to consider the production and the maintenance costs. For example,
deterioration of manufacturing units can significantly decrease the production
quality [3]. The quality of the production depends, among others, on the
maintenance strategy of the manufacturing equipment to prevent production
disruptions [4], [5]. In fact, the goal of maintenance scheduling is to
maximize the availability of the system, while fulfilling the demand of the
users or customers. Traditionally, preventive maintenance is scheduled
regularly to restore the units to their healthy states. However, preventive
maintenance does not consider the system’s condition and the dynamic changes
[6]. Predictive maintenance has the potential to significantly reduce the
maintenance costs. Furthermore, maintenance scheduling does typically not
consider the consumption behaviour of its customers which influences the
prizing but also the produced amount which in turn influences the requirements
on the availability of the system but also the degradation of the system and
the demand for maintenance actions. Hence, in order to maximize the
availability and reward function, it is advisable for a supplier to schedule
the maintenance of its manufacturing units using the underlying graph of the
social network of the users and the users’ consumption function.

In this paper, we address the maintenance scheduling problem of a
manufacturing systems of goods with positive externalities. In this context,
the goal of the manufacturer is to find the optimal price of the goods and the
optimal maintenance scheduling for the manufacturing units such that its
revenue and availability are maximized. In order to obtain the maintenance
scheduling, we consider the framework of predictive maintenance in which the
remaining useful lifetime (RUL) of the units and the associated uncertainties
can be predicted. We define the optimization model for the manufacturer which
comprises two main elements: 1) the revenue that can be obtained by selling
the produced goods and 2) maintenance cost of the manufacturing units. To
model the positive externality effects, we assume that users or customers are
connected through a graph and that they choose their consumption level of the
good based on their utility function and the price of the good. Hence, we face
a bi-level optimization problem with the two levels representing the supplier
and the users. In order to solve the optimization problem, we propose to
formulate the problem as a leader-follower game between the supplier and the
users. In this game, the supplier is the leader and can predict the strategies
(consumption function) of the followers. Then, based on this prediction, it
determines the maintenance scheduling and optimal price of the network and
sets the price for each of the users. We investigate the efficiency of the
proposed approach based on a case study.

To the best of our knowledge, this is the first research that obtains the
maintenance scheduling of different units in the manufacturing system by
considering the price of the network and the positive externality effects of
the goods using a leader-follower game among the manufacturer and the
customers.

The rest of the paper is organized as follows. The review of related work is
presented in Section 2. The preliminaries on predicting the remaining useful
lifetime (RUL) and the leader-multiple-followers game are introduced in
Section 3. The supplier’s and customers’ objective functions are formulated in
Section 4. The solution method based of the leader-follower game is proposed
in Section 5. The simulation results are presented in Section 6. The
conclusion remarks are made in Section 7.

##  2\. Related work

Network price of goods with positive externalities. Several approaches have
been proposed to determine the optimal price in networks with positive
externalities. The authors of [2] study the optimal pricing strategies of a
supplier of goods with positive externalities assuming that consumers are
connected by a social network graph. Their usage levels depend directly on the
usage of their neighbors within a social network [2]. Since the users’
consumption behaviour affects that of other users, it has been modeled based
on the game theory concept [7], [8]. The optimal contracting between the
seller and buyers is proposed in [7]. The authors develop the solution method
using the Stackelberg game where the seller acts as a leader and the buyers
are followers. The strategic interaction and networks have been addressed
using game theory in [8]. In this study, the decisions of some suppliers
affect those of other suppliers in the network.

In all of these works, the manufacturers solely consider the pricing of the
goods and do not take the impact of the production on the operation and
maintenance costs into account. In fact, sudden failures causing an
interruption of the production can severally impact the ability of the
supplier to respond to the demand. Furthermore, not considering the
interruptions caused by planned maintenance in the pricing will result in sub-
optimal prices and revenues. To maximize their revenue, the suppliers should
strive to perform maintenance on their manufacturing units when the price of
the network is low and they sell less.

Maintenance scheduling. Maintenance scheduling with the aim of increasing the
system reliability and availability has been addressed in many applications
such as electrical market [9], [10], network transportation [11], [12], and
manufacturing system [13], [14].

Several studies have proposed different approaches for maintenance scheduling
and optimal management of manufacturing systems. The authors of [13] propose a
method based on a Genetic algorithm to obtain the maintenance scheduling of
manufacturing equipment using the predicted level of degradation. The
maintenance scheduling for manufacturing systems proposed in [15] considers
the product quality and the maintenance cost including the repair and
preventive maintenance cost. Dynamic maintenance scheduling which obtains the
optimal predictive maintenance using dynamic programming method is developed
in [16], [17], and [18]. In these studies, the deterioration process of the
units is formulated as a Markov decision process. An optimization-based
approach for dynamic maintenance scheduling of different units which considers
the operating conditions is proposed in [14]. In all of these studies, the
maintenance scheduling is obtained only based on the degradation level of
different units without considering the price of the goods and the effect of
maintenance on the revenue of the manufacturer.

####  Notation:

Given $D\in{\mathbb{R}^{n\times{n}}}$, $D^{-1}$ denotes the inverse of $D$. We
define the column augmentation of $Z_{n}(t)$ for $t=1,\cdots,T$, as
${Z}_{n}=\textbf{col}(Z_{n}(1),\cdots,Z_{n}(T)):=[Z_{n}(1),\cdots,Z_{n}(T)]$.

##  3\. Preliminaries

###  3.1. Remaining useful life (RUL) and predictive maintenance

The remaining useful life (RUL) is defined as the amount of time that an asset
will continue to satisfy its desired operating conditions [19]. Predicting RUL
is one of the core tasks in predictive maintenance applications. Predicting
the RUL enables on the one hand to perform maintenance before the failure
occurs and by that also to improve the availability of the system while on the
other hand, the lifetime of the system can be fully exploited resulting in
less frequent replacements or repairs. Several methods have been proposed in
the literature to estimate RUL. These methods can be categorized in three
categories: 1) model-based approaches, where the physics principle models of
the degradation of the asset are applied [20]; 2) data-driven methods, in
which the RUL is estimated based on the condition monitoring data only [21];
3) knowledge-based approaches that depend on the domain knowledge of an expert
[22]. RUL estimations are always subject to uncertainties. While not all of
the RUL prediction methods also estimate the uncertainty of the predictions,
uncertainty quantification is an integral part of decision support systems and
is particularly desirable by domain experts using such systems. For many
systems with multiple-units that are jointly fulfilling the demand, performing
predictive maintenance at the end of life of the system may not the optimal.
Since, the maintenance decision of each unit also depends on the decision of
other units and the production requirements. Hence, predictive maintenance
scheduling is pivotal for fulfilling the demand while minimizing the
maintenance costs. In our problem formulation, we assume that we are able to
predict the RUL and can quantify the associated uncertainty. We assume that
the distribution function of RUL is given. The manufacturer should, therefore,
take into account the uncertainty of RUL and schedule the maintenance based on
the production requirements and the associated costs of the unavailability,
while considering that a part of the lifetime of the unit may be lost when
replacing the unit before the end of life.

###  3.2. Leader-multiple-followers game

Let us consider N-players, $i=1,\cdots,N$, as followers and one player as a
leader. Let $x_{i}$ denote the follower $i$’s strategy and $y$ the leader’s
strategy, respectively. Let $U^{f}_{i}(x_{i},x_{-i},y)$ denote the utility
function of the follower $i$ where $x_{-i}$ is the strategy of all followers
except follower $i$. Let $U^{l}(y,x)$, in which
$x=\textbf{col}(x_{1},\cdots,x_{n})$, denote the utility function of the
leader. The equilibrium point of the leader-multiple-followers game can be
defined as follows:

######  Definition 1.

Let $(x^{*},y^{*})$ be the equilibrium point among the leader and the
followers, then,

$\displaystyle
U^{f}_{i}(x^{*}_{i},x^{*}_{-i},y^{*})\geq{U^{f}_{i}(x_{i},x^{*}_{-i},y^{*})},\quad
i=1,\cdots,N,$ (1) $\displaystyle U^{l}(y^{*},x^{*})\geq{U^{l}(y,x^{*})}$

In the leader-follower game, the leader predicts the followers’ strategies and
implements his strategy first, while the followers react to the leader’s
decisions.

##  4\. Problem formulation

In this section, we propose the optimization model for obtaining the
maintenance scheduling of the units in the manufacturing system based on the
optimal price of the network.

Let us define $\mathcal{J}={\\{1,\cdots,J}\\}$ as the set of units in the
manufacturing system that are manufacturing the good that has positive
externalities to the customers. Let us define $\mathcal{N}={1,\cdots,N}$ as
the set of customers who are connected through the graph matrix $W$ in the
social network. The $il$ element of matrix $W$ denoted by $w_{il}\geq{0}$
represents the strength of the connection between the two customers in the
network and represents concurrently also the influence of customer $l$ on
customer $i$. The manufacturer aims to obtain the optimal price of the network
while maintaining a high availability of the manufacturing units. We assume
that the manufacturer has implemented a predictive maintenance strategy to
improve the availability of the system and is able to predict the RUL of the
systems and the associated uncertainty. Predicting the RUL enables
manufacturer to schedule the maintenance before the end of life ant thereby,
prevent unscheduled down times while fully exploiting the lifetimes of the
units. The manufacturers aims on the one hand to maintain a high availability
of the units and exploit the useful lifetime of the units as much as possible,
while performing the maintenance at points in time when it has the least
possible impact on its revenue. The manufacturer is, therefore, seeking to
obtain the maintenance scheduling of manufacturing units for the decision
horizon time $\mathcal{T}={\\{1,\cdots,T}\\}$ while maximizing the revenue
that can be obtained by selling the products to the customers. In the
following, we explain the proposed optimization model.

###  4.1. Manufacturer’s model: maintaining and pricing

In the following, we propose the manufacturer’s optimization problem which
concurrently optimizes the price of the network and the maintenance scheduling
of the units.

Maintaining: In this study, we assume that the deterioration state of the unit
increases with time. The maintenance restores the state of the unit to its
initial condition (as good as new). Hence, we model the deterioration and the
restoration as follows:

$s_{j}(t+1)=(1-x_{j}(t))s_{j}(t)+1,$ (2)

where $s_{j}(t)$ and $s_{j}(t-1)$ denote the deterioration state of unit $j$
at instants $t$ and $t-1$, respectively. $x_{j}(t)$ is a binary variable and
$x_{j}(t)=1$ denotes that unit $j$ performs maintenance at instant $t$.

When the deterioration state reaches the degradation threshold, it reaches the
end of life which is equivalent to the remaining useful life (RUL) reaching
zero. If the unit has not been maintained before the end of life, it will
fail. As discussed above, the prediction of the RUL is always associated with
uncertainties. It can, therefore, be considered as a random variable with a
known distribution function. We assume that the predicted RUL and the
associated uncertainty are given. To avoid failure, the units should perform
maintenance before the end of life. To model this aim, we apply the chance
constraint problem as follows:

${\mathbb{P}}{\Big{(}S_{2,j}\in\mathcal{A}|\quad
s_{j}(t)-{S_{2,j}}\leq{0}\Big{)}}\geq{1-\alpha},$ (3)

where $\mathbb{P}$ is a probability measure defined over $\mathcal{A}$.
$S_{2,j}$ is the degradation threshold of the deterioration state of unit $j$.
$s_{j}(t)$ which satisfies the chance constraint (3) can be considered as the
$\alpha$-level feasible solution.

Another constraint of this problem is that the total consumption of all
customers should be less than the total production of all units. We assume
that the demand cannot be postponed and is satisfied at all times.
Furthermore, we assume that the good is perishable and cannot be stored to
cover the demand during the maintenance down time. Moreover, we assume that
when the unit performs maintenance, its production is zero. Hence, we can
model this constraint as follows:

$\sum\limits_{i\in\mathcal{N}}q_{i}(t)\leq{\sum\limits_{j\in\mathcal{J}}(1-x_{j}(t))q_{j,max}},\\\
$ (4)

where $q_{i}(t)$ is the amount of the consumption of customer $i$ at instant
$t$. $q_{j,\max}$ is the maximum production output by manufacturing unit $j$.

The supplier seeks to obtain the maintenance scheduling which satisfies (3)
and (4), while concurrently minimizing the maintenance cost of the units which
is modeled as follows:

$C=\sum\limits_{j\in\mathcal{J}}\sum\limits_{t\in\mathcal{T}}c_{j}x_{j}(t),\\\
$ (5)

where $c_{j}$ is the cost of unit $j$.

Pricing: The supplier aims to obtain the price of the network which maximizes
the following reward:

$R=\sum\limits_{i\in\mathcal{N}}\sum\limits_{t\in\mathcal{T}}\phi_{i}(t){q_{i}(t)},\\\
$ (6)

where $\phi_{i}(t)\geq{0}$ is the price of the consumed good of customer $i$
at instant $t$. $q_{i}(t)\geq{0}$ is the consumption of customer $i$ at
instant $t$.

Supplier’s objective function: Using (5) and (6), the utility function of the
supplier can be expressed as follows:

$U=\sum\limits_{i\in\mathcal{N}}\sum\limits_{t\in\mathcal{T}}\phi_{i}(t){q_{i}(t)}-\sum\limits_{j\in\mathcal{J}}\sum\limits_{t\in\mathcal{T}}c_{j}x_{j}(t),$
(7)

Hence, by applying (7), and by considering constraints (2), (3), and (4), we
can formulate the supplier’s objective function as follows:

$\displaystyle\max\limits_{{\phi},{x},{s}}\sum\limits_{i\in\mathcal{N}}\sum\limits_{t\in\mathcal{T}}\phi_{i}(t){q_{i}(t)}-\sum\limits_{j\in\mathcal{J}}\sum\limits_{t\in\mathcal{T}}c_{j}x_{j}(t)$
(8) $\displaystyle\text{S.b.}\quad\text{C}_{1}:\quad
s_{j}(t+1)=(1-x_{j}(t))s_{j}(t)+1,$
$\displaystyle\quad\qquad\text{C}_{2}:\quad{\mathbb{P}}{\Big{(}S_{2,j}\in\mathcal{A}|\quad
s_{j}(t)-{S_{2,j}}\leq{0}\Big{)}}\geq{1-\alpha},$
$\displaystyle\quad\qquad\text{C}_{3}:\quad\sum\limits_{i\in\mathcal{N}}q_{i}(t)\leq{\sum\limits_{j\in\mathcal{J}}(1-x_{j}(t))q_{j,max}},$
$\displaystyle\quad\qquad\text{C}_{4}:\quad
q_{i}(t)\geq{0},\quad\phi_{i}(t)\geq{0},$

where ${\phi}=\textbf{Col}{\\{\phi_{1},\cdots,\phi_{N}}\\}$,
${\phi}_{i}=\textbf{Col}{\\{\phi_{i}(1),\cdots,\phi_{i}(T)}\\}$,
${x}=\textbf{Col}{\\{x_{1},\cdots,x_{N}}\\}$,
${x}_{i}=\textbf{Col}{\\{x_{i}(1),\cdots,x_{i}(T)}\\}$,
${s}=\textbf{Col}{\\{s_{1},\cdots,s_{N}}\\}$,
${s}_{i}=\textbf{Col}{\\{s_{i}(1),\cdots,s_{i}(T)}\\}$, $i\in\mathcal{N}$.

The problem (8) is a nonlinear mixed-integer programming problem due to the
constraints $\text{C}_{1}$ and $\text{C}_{2}$. We convert (8) to the mixed-
integer linear programming by using the big-M method [23] and scenario-based
approach [24]. In the following, we explain these two approaches.

First, let us introduce a new variable $y_{j}(t)=x_{j}(t)s_{j}(t)$. By using
the big-$M$ method, we can convert constraint $\text{C}_{1}$ to the following
mixed-integer linear constraints:

$\displaystyle s_{j}(t+1)=s_{j}(t)-y_{j}(t)+1,$ (9) $\displaystyle
y_{j}(t)\geq{s_{j}(t)-(1-x_{j}(t))M},$ $\displaystyle
y_{j}(t)\leq{s_{j}(t)+(1-x_{j}(t))M},$
$\displaystyle{0}\leq{y_{j}(t)}\leq{x_{j}(t)}M.$

Now, we can substitute the chance constraint (3) with the finite number of
constraints using the scenario based approach. Each constraint corresponds to
a different realization of $S_{2,j}^{k}$, $k=1,\cdots,K$, for the uncertain
parameter $S_{2,j}$. Thus, we have:

$s_{j}(t)-S_{2,j}^{k}\leq{0},\quad k=1,\cdots,K.$ (10)

######  Remark 1.

The number of the scenarios $K$ should be chosen sufficiently large such that
the feasible solution of (10) is an $\alpha$-level feasible solution of (3).

Hence, by using (9) and (10), the optimization problem (8) leads to the
following mixed-integer linear programming problem:

$\displaystyle\max\limits_{{\phi},{x},{s},{y}}\sum\limits_{i\in\mathcal{N}}\sum\limits_{h\in\mathcal{H}}\phi_{i}(t){q_{i}(t)}-\sum\limits_{j\in\mathcal{J}}\sum\limits_{t\in\mathcal{T}}c_{j}x_{j}(t)$
(11) $\displaystyle\text{S.b.}\quad\text{C}^{\prime}_{1}:\quad
s_{j}(t+1)=s_{j}(t)-y_{j}(t)+1,$
$\displaystyle\quad\qquad\text{C}^{\prime}_{2}:\quad
y_{j}(t)\geq{s_{j}(t)-(1-x_{j}(t))M},$
$\displaystyle\quad\qquad\text{C}^{\prime}_{3}:\quad
y_{j}(t)\leq{s_{j}(t)+(1-x_{j}(t))M},$
$\displaystyle\quad\qquad\text{C}^{\prime}_{4}:{0}\leq{y_{j}(t)}\leq{x_{j}(t)}M.$
$\displaystyle\quad\qquad\text{C}^{\prime}_{5}:\quad
s_{j}(t)-S_{2,j}^{k}\leq{0},\quad k=1,\cdots,K.$
$\displaystyle\quad\qquad\text{C}^{\prime}_{6}:\quad\sum\limits_{i\in\mathcal{N}}q_{i}(t)\leq{\sum\limits_{j\in\mathcal{J}}(1-x_{j}(t))q_{j,max}},$
$\displaystyle\quad\qquad\text{C}^{\prime}_{7}:\quad
q_{i}(t)\geq{0},\quad\phi_{i}(t)\geq{0},$

where ${y}=\textbf{Col}{\\{y_{1},\cdots,y_{N}}\\}$,
${y}_{i}=\textbf{Col}{\\{y_{i}(1),\cdots,y_{i}(T)}\\}$.

The consumption $q_{i}(t)$ is obtained by the customers’ objective function,
which is explained in the following.

###  4.2. Customers’ model: consumption

The customers seek to obtain the optimal consumption which maximizes their
reward function. Let us model the customer $i$ objective function as follows:

$\max\limits_{{q}_{i}}\sum\limits_{t\in\mathcal{T}}\big{(}v_{i}(q_{i}(t))+\sum\limits_{l\in\mathcal{N}}w_{il}q_{l}(t)q_{i}(t)-\phi_{i}(t){q}_{i}(t)\big{)},$
(12)

where ${q}_{i}=\textbf{Col}{\\{q_{i}(1),\cdots,q_{i}(T)}\\}$. The objective
function consists of three parts. The first term denotes the reward that
customer $i$ can obtain by consuming the good. Inspired by the studies in [2]
and [1], we model the first reward term as follows:

$v_{i}(q_{i}(t))=-\frac{1}{2}a_{i}q_{i}^{2}(t)+b_{i}q_{i}(t),$ (13)

where $a_{i}\geq{0}$ and $b_{i}$ are the customer’s $i$ model. The second term
expresses the effect of other customers which are connected through the social
network graph to customer $i$. This term implies that the consumption
$q_{k}(t)$ of customer $k$, who is connected to customer $i$ through the
social network graph with the weight $w_{ik}$, influences the objective
function of customer $i$. The third term formulates the price that customer
$i$ is willing to pay to the supplier for obtaining the consumed good
$q_{i}(t)$.

##  5\. Proposed solution: leader and multiple-followers game

In this section, we propose the solution method using the concept of leader-
multiple-followers game [25], [26] for solving the bi-level optimization
problem ((11) and (12)). In this problem, we consider the customers as the
followers who are seeking to obtain their consumption and the supplier as the
leader who is responsible to obtain the price of the network and maintenance
scheduling of its manufacturing units. The schematic of the proposed framework
is shown in Figure 1.

Figure 1. Schematic of the leader-multiple-followers game for obtaining the
maintenance scheduling based on the price of the network.

###  5.1. Consumption equilibrium of multiple-follower

In the proposed solution, we assume that the customers are the followers of
the supplier. According to (12), since the customers’ objective function
depends on the strategies of other customers, we face a game theory problem
among the customers which can be described as
$G=(\mathcal{N},\mathcal{Q},\mathcal{U})$:

1) $\mathcal{N}$ denotes the set of customers as the players.

2) $\mathcal{Q}=\prod\limits_{i\in\mathcal{N}}q_{i}$ denotes the strategy
space of the players.

3) $\mathcal{U}={\\{U_{1},\cdots,U_{N}}\\}$ is the set of utilities, where the
utility of player $i$ is defined as follows:

$U_{i}(q_{i},q_{-i},\phi_{i})=\sum\limits_{t\in\mathcal{T}}\big{(}v_{i}(q_{i}(t))+\sum\limits_{l\in\mathcal{N}}w_{il}q_{l}(t)q_{i}(t)-\phi_{i}(t){q}_{i}(t)\big{)}$
(14)

where $q_{-i}={\\{q_{1},\cdots,q_{i-1},q_{i+1},\cdots,q_{N}}\\}$ are the
strategies of all players except player $i$.

The Nash equilibrium (NE) is one appropriate output solution of the game.
Using Definition (1), the NE among customers can be defined as follows:

$U_{i}(q^{*}_{i},q^{*}_{-i},\phi^{*}_{i})\geq{U_{i}(q_{i},q^{*}_{-i},\phi^{*}_{i})}$
(15)

In order to find the NE of the game among customers, let us define the
following assumption:

######  Assumption 1.

$a_{i}\geq{\sum\limits_{l\in\mathcal{N}}w_{il}}$, $i\in\mathcal{N}$.

######  Theorem 1.

Under Assumption (1), the game $G=(\mathcal{N},\mathcal{Q},\mathcal{U})$ has a
unique NE which is defined as follows:

$q^{*}(t)=(A-W)^{-1}(B-\Phi^{*}(t)),$ (16)

where $A=\mathop{\mathrm{diag}}(a_{1},\cdots,a_{N})$,
$W=\begin{pmatrix}w_{11}&\cdots&w_{1N}\\\ \vdots&\cdots&\vdots\\\
w_{N1}&\cdots&w_{NN}\end{pmatrix}$, $B=\begin{pmatrix}b_{1}\\\ \vdots\\\
b_{N}\end{pmatrix}$, $\Phi^{*}=\begin{pmatrix}\phi^{*}_{1}\\\ \vdots\\\
\phi^{*}_{N}\end{pmatrix}$

###### Proof.

The best response of player $i$ is as follows:

$q_{i}(t)=\beta_{i}(q_{-i}(t))=\frac{b_{i}-\phi^{*}_{i}}{a_{i}}+\frac{\sum\limits_{l\in\mathcal{N}}w_{il}q_{l}(t)}{a_{i}},\quad{i\in\mathcal{N}},$
(17)

Hence, (17) can be written in the matrix form as follows:

$Aq({t})=B-\Phi^{*}(t)+Wq(t)$ (18)

Under Assumption (1), matrix $A-W$ is invertible. Hence, the NE of the game
can be obtained as follows:

$q^{*}(t)=(A-W)^{-1}(B-\Phi^{*}(t)).$ (19)

∎

###  5.2. Price and maintenance scheduling of the leader

As we mentioned above, in the proposed model, we assume that the supplier can
predict the strategies of followers. Based on this prediction, the supplier
optimizes its strategy. Moreover, the supplier acts first, in the sense that
it sends its strategy to the followers, which then decide about their actions.
We assume that the pricing is performed at the level of an individual
customer. Customers who are more connected within the social network (and are,
therefore, impacted by the network) can than be better targeted by the
supplier. The strategy of the supplier comprises the price of the network and
the maintenance scheduling of the manufacturing units. Hence, in the
following, we obtain the optimal price and maintenance scheduling of the
supplier relying on the fact that it knows the NE of the game (16) among the
customers.

Let us define the matrix $(A-W)^{-1}$ as follows:

$(A-W)^{-1}=(R_{il})_{i,l\in\mathcal{N}}.$ (20)

In view of equality (16) and using (20), the optimization problem of supplier
(11) can be defined as the following mixed-integer quadratic programming
problem:

$\displaystyle\max\limits_{{\phi},{x},{s},{y}}\sum\limits_{t\in\mathcal{T}}\phi^{\top}(t)(A-W)^{-1}(B-\Phi(t))-\sum\limits_{j\in\mathcal{J}}\sum\limits_{t\in\mathcal{T}}c_{j}x_{j}(t)$
(21)
$\displaystyle\text{S.b.}\quad\text{C}^{\prime}_{1},\quad\text{C}^{\prime}_{2}\quad,\quad\text{C}^{\prime}_{3},\quad\text{C}^{\prime}_{4},\quad\text{C}^{\prime}_{5},$
$\displaystyle\qquad\quad\text{C}^{\prime\prime}_{6}:\quad\sum\limits_{l\in\mathcal{N}}R_{il}(b_{l}-\phi_{l})\leq{\sum\limits_{j\in\mathcal{J}}(1-x_{j}(t))q_{j,max}},$
$\displaystyle\quad\qquad\text{C}^{\prime\prime}_{7}:\quad(A-W)^{-1}(B-\Phi(t))\geq{0},\quad\phi(t)\geq{0}.$

By solving (21), the optimal price $\phi^{*}_{i}$, $i\in\mathcal{N}$, can be
obtained and sent to the customers. Based on the price, they can obtain their
NE using (16).

##  6\. Simulation results

In this section, we implement the proposed method on a manufacturing system
where the supplier has five manufacturing units and aims to obtain their
maintenance scheduling for $30$ days. We assume a normal distribution for the
deterioration threshold value. The mean and the standard deviation of the
deterioration state thresholds and the maintenance costs of each of the
manufacturing units are displayed in Table 1.

Table 1. Deterioration and maintenance parameters of the manufacturing units Unit | Mean of deterioration | Standard deviation of deterioration | Maintenance cost
---|---|---|---
| state thresholds $(S_{2,j})$ | state thresholds $(S_{2,j})$ |  $(c_{j})$ $[MU]$
1 | 12 | 1.4 | 20.48
2 | 10 | 3.2 | 21.39
3 | 11.2 | 2.5 | 22.73
4 | 9.4 | 1.1 | 24.78
5 | 11.8 | 2.1 | 24.82

Moreover, we assume $10$ customers who are connected through a social network
graph. The profile of the customers’ demand is depicted in Figure 2.

Figure 2. The total demand of customers during $30$ days.

We consider the following graph weights among the customers in the network:

$W=\begin{pmatrix}0&0.54&0.82&0.28&0.11&0.41&0.16&0.05&0.76&0.86\\\
0.54&0&0.79&0&0&0&0&0&0&0\\\ 0.82&0.79&0&0.36&0&0&0&0&0&0\\\
0.28&0&0.36&0&0&0&0&0&0\\\ 00.11&0&0&0.08&0&0.62&0&0&0&0\\\
0.41&0&0&0&0.62&0&0&0&0\\\ 0.16&0&0&0&0&0.63&0&0.75&0&0\\\
0.05&0&0&0&0&0&0.75&0&0.54&0\\\ 0.76&0&0&0&0&0&0&0.54&0&0.25\\\
0.86&0&0&0&0&0&0&0&0.25&0\end{pmatrix}$ (22)

The weight matrix shows that customer $1$ is connected to all other customers.
Hence, his/her decisions have a high impact on the decisions of other
customers. Figure 3 shows the consumption level of the customers. As we can
see, the consumption level of customer $1$ is higher than that of other
customers since the supplier offers a lower price to this customer to
encourage him/her to consume more goods which has a high impact on other
customers’ strategies.

Figure 3. The total consumption of customers during $30$ days.

The network price and maintenance scheduling of the units are shown in Figures
4 and 5.

Figure 4. The profile of the network price during $30$ days.

As we can see from Figure 4 and 5, the manufacturing units perform maintenance
more frequently when the price of the network is low (in the first $15$ days)
compared to the time periods when the network price is high (between day $15$
and $25$).

Figure 5. Maintenance scheduling of all supplier’s manufacturing units during
$30$ days.

The deterioration state of the units is shown in Figure 6. The figure shows
that the units perform maintenance before their failure time

Figure 6. Deterioration state of all supplier’s units during $30$ days.

###  6.1. The effect of network structure

In this section, we investigate the effect of the social network structure on
the supplier’s objective function. In case that the weights of the network are
not known to the supplier, the supplier predicts the NE of the network as
follows:

$q^{*}(t)=A^{-1}(B-\phi^{*}(t))$ (23)

Then, the price and maintenance scheduling are obtained as follows:

$\displaystyle\max\limits_{{\phi},{x},{s},{y}}\sum\limits_{t\in\mathcal{T}}\phi^{\top}(t)(A)^{-1}(B-\Phi(t))-\sum\limits_{j\in\mathcal{J}}\sum\limits_{t\in\mathcal{T}}c_{j}x_{j}(t)$
(24)
$\displaystyle\text{S.b.}\quad\text{C}^{\prime}_{1},\quad\text{C}^{\prime}_{2}\quad,\quad\text{C}^{\prime}_{3},\quad\text{C}^{\prime}_{4},\quad\text{C}^{\prime}_{5},$
$\displaystyle\qquad\quad\text{C}^{\prime\prime}_{6}:\quad\sum\limits_{l\in\mathcal{N}}R_{il}(b_{l}-\phi_{l})\leq{\sum\limits_{j\in\mathcal{J}}(1-x_{j}(t))q_{j,max}},$
$\displaystyle\quad\qquad\text{C}^{\prime\prime}_{7}:\quad(A)^{-1}(B-\Phi(t))\geq{0},\quad\phi(t)\geq{0}.$

In our case study, the profit that the supplier can obtain by solving (21)
(knowing the network structure) for the entire duration of the period over
$30$ days is $1.2893e+05$ $MU$ (monetary unit) and the supplier’s objective
function without knowing the network structure (24) is $5.3622e+04$ $MU$.
Hence, by knowing the network structure, a supplier can obtain the strategy
that results in a $7.5308e+04$ $MU$ higher profit.

###  6.2. Comparison of the proposed method to the baseline maintenance
decision

The result of the reward function is compared to the baseline solution where
the units perform maintenance at the deterioration threshold
$\min\limits_{k=1,\cdots,K}S^{k}_{2,j}$, $j\in\mathcal{J}$. This baseline
decision makes sure that the units can perform maintenance before they fail.
The maintenance decision and the degradation cost are shown in Figures 7 and
8.

Figure 7. Baseline maintenance scheduling of all supplier’s manufacturing
units during $30$ days. Figure 8. Deterioration state of all supplier’s units
for baseline maintenance decision during $30$ days .

When the units decide to perform maintenance when the threshold is reached,
the profit of the supplier is $7.4919e+04$ $MU$, which is $5.4011$ $MU$ less
than the profit when they obtain their maintenance decision using (11). In
this case, the maintenance scheduling of units does not depend on the price of
the network but solely on the deterioration state.

##  7\. Conclusion

In this paper, we address the problem of maintenance scheduling for
manufacturing units of goods with positive externalities by concurrently
optimizing the pricing and the maintenance schedule. The customers are
connected through a social network graph. In order to solve the problem, we
propose a bi-level optimization approach where at the first level, the
supplier obtains its maintenance scheduling and the price at the level of
individual customers. At the second level, the customers obtain their
strategies based on the price offered by the supplier. In order to solve the
bi-level optimization problem, we propose a solution based on a leader-
multiple-followers game where the customers are the followers and the supplier
acts as a leader that optimizes its decisions based on the predicted
strategies of the customers. The proposed approach can be extended and
implemented to the case of more complex network with more agents. The
numerical results of the case study show that when the network structure is
known to the supplier, the obtained profit of the supplier is higher compared
to the case that the connections between the customers in the network are not
known. Moreover, we demonstrated that obtaining the maintenance decision using
the proposed bi-level optimization method leads to more profit for the
supplier than base-line solution.

As future research, multiple suppliers can be considered in the proposed
framework where the strategies of the suppliers are related to each other and
they also play a game. Moreover, one can consider uncertainties in the RUL
prediction and the customers’ models. This will enable to develop a robust
optimization method which can take the uncertainties into account.

## Acknowledgement

The contribution of Olga Fink was funded by the Swiss National Science
Foundation (SNSF) Grant no. PP00P2_176878.

## References

  * [1] A. Jadbabaie, A. Kakhbod, Optimal contracting in networks, Journal of Economic Theory 183 (2019) 1094–1153.
  * [2] O. Candogan, K. Bimpikis, A. Ozdaglar, Optimal pricing in networks with externalities, Operations Research 60 (4) (2012) 883–905.
  * [3] B. Lu, X. Zhou, Quality and reliability oriented maintenance for multistage manufacturing systems subject to condition monitoring, Journal of Manufacturing Systems 52 (2019) 76–85.
  * [4] A. Sarkar, S. Panja, B. Sarkar, Survey of maintenance policies for the last 50 years, International Journal of Software Engineering and Applications 2 (3) (2011) 130–148.
  * [5] İ. Erozan, A fuzzy decision support system for managing maintenance activities of critical components in manufacturing systems, Journal of Manufacturing Systems 52 (2019) 110–120.
  * [6] E. Omshi, A. Grall, S. Shemehsavar, A dynamic auto-adaptive predictive maintenance policy for degradation with unknown parameters, European Journal of Operational Research 282 (1) (2020) 81–92.
  * [7] D. J. Wu, P. R. Kleindorfer, J. E. Zhang, Optimal bidding and contracting strategies for capital-intensive goods, European Journal of Operational Research 137 (3) (2002) 657–676.
  * [8] Y. Bramoullé, R. Kranton, M. D’amours, Strategic interaction and networks, American Economic Review 104 (3) (2014) 898–930.
  * [9] P. Rokhforoz, B. Gjorgiev, G. Sansavini, O. Fink, Multi-agent maintenance scheduling based on the coordination between central operator and decentralized producers in an electricity market, Reliability Engineering & System Safety 210 (2021) 107495.
  * [10] A. Volkanovski, B. Mavko, T. Boševski, A. Čauševski, M. Čepin, Genetic algorithm optimisation of the maintenance scheduling of generating units in a power system, Reliability Engineering & System Safety 93 (6) (2008) 779–789.
  * [11] A. Amiri, Designing a distribution network in a supply chain system: Formulation and efficient solution procedure, European journal of operational research 171 (2) (2006) 567–576.
  * [12] P. Feng, Y. Liu, F. Wu, C. Chu, Two heuristics for coordinating production planning and transportation planning, International Journal of Production Research 56 (21) (2018) 6872–6889.
  * [13] Z. M. Yang, D. Djurdjanovic, J. Ni, Maintenance scheduling in manufacturing systems based on predicted machine degradation, Journal of intelligent manufacturing 19 (1) (2008) 87–98.
  * [14] M. Celen, D. Djurdjanovic, Integrated maintenance and operations decision making with imperfect degradation state observations, Journal of Manufacturing Systems 55 (2020) 302–316.
  * [15] B. Lu, X. Zhou, Opportunistic preventive maintenance scheduling for serial-parallel multistage manufacturing systems with multiple streams of deterioration, Reliability Engineering & System Safety 168 (2017) 116–127.
  * [16] A. Ghasemi, S. Yacout, M. Ouali, Optimal condition based maintenance with imperfect information and the proportional hazards model, International journal of production research 45 (4) (2007) 989–1012.
  * [17] A. Ghasemi, S. Yacout, M. Ouali, Optimal strategies for non-costly and costly observations in condition based maintenance, International Journal of Applied Mathematics 38 (2) (2008).
  * [18] S. Kröning, B. Denkena, Dynamic scheduling of maintenance measures in complex production systems, CIRP Journal of Manufacturing Science and Technology 6 (4) (2013) 292–300.
  * [19] J. Lee, F. Wu, W. Zhao, M. Ghaffari, L. L. D. Siegel, Prognostics and health management design for rotary machinery systems—reviews, methodology and applications, Mechanical systems and signal processing 42 (1-2) (2014) 314–334.
  * [20] G. Liao, H. Yin, M. Chen, Z. Lin, Remaining useful life prediction for multi-phase deteriorating process based on wiener process, Reliability Engineering & System Safety 207 (2021) 107361.
  * [21] X. Si, W. Wang, C. Hu, D. Zhou, M. Pecht, Remaining useful life estimation based on a nonlinear diffusion degradation process, IEEE Transactions on reliability 61 (1) (2012) 50–67.
  * [22] M. Djeziri, S. Benmoussa, R. Sanchez, Hybrid method for remaining useful life prediction in wind turbine systems, Renewable Energy 116 (2018) 173–187.
  * [23] J. Fortuny-Amat, B. McCarl, A representation and economic interpretation of a two-level programming problem, Journal of the operational Research Society 32 (9) (1981) 783–792.
  * [24] K. Margellos, P. Goulart, J. Lygeros, On the road between robust optimization and the scenario approach for chance constrained optimization problems, IEEE Transactions on Automatic Control 59 (8) (2014) 2258–2263.
  * [25] M. Hu, M. Fukushima, Existence, uniqueness, and computation of robust nash equilibria in a class of multi-leader-follower games, SIAM Journal on optimization 23 (2) (2013) 894–916.
  * [26] J. Cruz, Leader-follower strategies for multilevel systems, IEEE Transactions on Automatic Control 23 (2) (1978) 244–255.


# Network topology mapping of Chemical Compounds Space

Georgios Tsekenis  Institute for Complex Systems, National Research Council,
Rome, Italy  Department of Molecular Sciences and Nanosystems (DMSN), “Ca’
Foscari” University of Venice, Venice, Italy<EMAIL_ADDRESS>Giulio Cimini
Physics Department and INFN, University of Rome Tor Vergata, Rome, Italy
Marinos Kalafatis  Department of Microbiology, University of Illinois at
Urbana-Champaign, Urbana, Illinois, USA  Achille Giacometti  Department of
Molecular Sciences and Nanosystems (DMSN), “Ca’ Foscari” University of Venice,
Venice, Italy  European Centre of Living Technologies (ECLT), “Ca’ Foscari”
University of Venice, Venice, Italy  Tommaso Gili  Networks Unit, IMT School
for Advanced Studies Lucca, Lucca, Italy  Guido Caldarelli  Department of
Molecular Sciences and Nanosystems (DMSN), “Ca’ Foscari” University of Venice,
Venice, Italy  Institute for Complex Systems, National Research Council, Rome,
Italy  European Centre of Living Technologies (ECLT), “Ca’ Foscari” University
of Venice, Venice, Italy  Rara Foundation - Sustainable Materials and
Technologies ETS, Venice, Italy

###### Abstract

We define bipartite and monopartite relational networks of chemical elements
and compounds using two different datasets of inorganic chemical and material
compounds, as well as study their topology. We discover that the connectivity
between elements and compounds is distributed exponentially for materials, and
with a fat tail for chemicals. Compounds networks show similar distribution of
degrees, and feature a highly-connected club due to oxygen. Chemical compounds
networks appear more modular than material ones, while the communities
detected reveal different dominant elements specific to the topology. We
successfully reproduce the connectivity of the empirical chemicals and
materials networks by using a family of fitness models, where the fitness
values are derived from the abundances of the elements in the aggregate
compound data. Our results pave the way towards a relational network-based
understanding of the inherent complexity of the vast chemical knowledge atlas,
and our methodology can be applied to other systems with the ingredient-
composite structure.

## Introduction

The space of chemical compounds comprises hundreds of thousands of different
combinations of the over one hundred chemical elements. Such an ample volume
was produced by employing several experimental and computational techniques
developed for the study of Chemistry over the past centuries. Navigating the
vast chemical space is a formidable task and has been the topic of previous
research (e.g. [1, 2]). Motivated by the need to harness the burgeoning
complexity of the ever-growing chemicals and materials fields, in this
manuscript, we present a constitutive relational network study of inorganic
chemistry and materials science, relying on the toolbox of complex networks
theory [3, 4].

In the past, chemical reaction networks have been presented for small numbers
of reactants [5], without addressing the overall complexity of the problem.
Furthermore, in materials science, recent efforts have concentrated on faster
and cheaper targeted engineering of materials, the so-called Materials Genome
project [6, 7]. Such an approach customarily relies on aggregate statistics.
However, incorporating meaningful relational networks can significantly
improve the inferential power of statistical approaches, such as materials
cartography [8]. One network approach has been based on the representation of
materials phase diagrams [9, 10]. A different approach was to analyze a set of
materials as a network, according to the cross-correlation of the electronic
density of states [11]. Unfortunately, these methods produce fully, or almost
fully, connected graphs where all substances are related, which is not very
different from an aggregate approach.

Here, we construct and study element-compound networks of extensive
catalogues/libraries of chemicals and materials. Furthermore, we successfully
model these networks with versatile fitness models derived from maximum
entropy methodology. That way, we set large bodies of knowledge onto a new
frame of reference, providing novel points of view and enabling further future
utility.

## Networks from Data

We construct relational networks from two different datasets that we sample
from two separate databases. The first, CRC, is based on inorganic chemical
compounds [12], and the second, AFLOW, is based on inorganic material
compounds [13] (see Methods). Each of the datasets contains $n_{C}$ compounds
and $n_{E}$ elements; the specific values are shown in Table 1.

dataset | compounds | elements | biLnks | cLinks | eLinks
---|---|---|---|---|---
| $n_{C}$ | $n_{E}$ | $L$ | $L_{C}$ | $L_{E}$
CRC | $3149$ | $89$ | $9118$ | $2117169$ | $1118$
AFLOW | $32243$ | $85$ | $97489$ | $106344446$ | $2956$
Table 1: Information about the datasets used to build all the networks.
$n_{C}$ and $n_{E}$ are the number of compounds and elements. $L$ is the
number of bipartite links, while $L_{C}$ and $L_{E}$ are the number of
monopartite links for compounds and elements, respectively.

We build a bipartite network for each dataset by linking every compound $c$ to
the elements $e$ it contains (Fig. 1a). For each dataset, the resulting
bipartite network is composed of two layers: one consisting of the compounds
$c$ and the other of the elements $e$, and is characterized by a $n_{C}\times
n_{E}$ binary bi-adjacency matrix $B$ linking the two layers [14, 15, 16, 17,
18], where the matrix element $B_{ce}=1$ if $c$ contains $e$ and zero
otherwise . For each bipartite network, the total number of links, $L$, is
given by the sum of all $B$ matrix elements: $L=\sum_{c,e}{B_{ce}}$. The
degree of a node is the sum of its incident links: $d_{c}=\sum_{e}B_{ce}$ and
$d_{e}=\sum_{c}B_{ce}$ for compounds and elements, respectively.

The degree distributions for both layers and both datasets are shown in the
four left panels of Fig. 2(a,b,e,f). The degrees ($d_{c}$) of the compounds
layer are discrete since each compound is linked to as many distinct elements
it contains. Their overall distribution appears modulated by a Gaussian-like
curve. The degree ($d_{e}$) distributions of the elements appear dominated for
larger values by a fat-tail for the CRC network, and by an exponential decay
for the AFLOW network, indicating a different inherent complexity of inorganic
chemicals vs materials. Oxygen is the most connected element corresponding to
the maximum degree, a feature confirmed also by other analyses as we shall see
below.

Figure 1: Explanatory visualizations of the chemical relational networks, with
an illustration of the linking processes on a tiny dataset of six compounds
(shown). a) The bipartite network of compounds and elements comprises two
layers, layer $C$ for the compounds and layer $E$ for the elements. b) The
monopartite projection $C$ of the chemical compounds, which are linked through
their common elements. c) Analogously, the monopartite projection $E$ of the
elements, which are linked through their co-participation in compounds. d) The
relative abundance, $a_{e}$, vs atomic number, $Z$, of the element species in
the CRC (green squares) and AFLOW (purple circles) datasets shows that the
most prominent elements are O, H, C for chemicals (green dotted vertical
lines), and O, S for materials (purple dotted vertical lines). e,f) The
cumulative distribution of element abundances $P_{>}(a_{e})$ appears with a
fat tail in CRC and an exponential tail in AFLOW. Oxygen is the most abundant
element in both datasets.

We further consider the relationships between compounds or between elements,
by projecting the bipartite network on either layer to get the corresponding
monopartite network. In the compounds network (Fig. 1b), the nodes are the
compounds, and a pair of compounds are linked if they share a common element.
In the elements network (Fig. 1c), the nodes are the elements, with links
between elements that co-participate in a compound. The adjacency matrices of
the binary monopartite networks, $A_{C}$ and $A_{E}$, are obtained by the
binary bi-adjacency matrix $B$: $(A_{C})_{cc^{\prime}}=1$ if
$\sum_{e}B_{ce}B_{c^{\prime}e}>0$, $(A_{E})_{ee^{\prime}}=1$ if
$\sum_{c}B_{ce}B_{ce^{\prime}}>0$, and zero otherwise. Summing all non-zero
entries in the adjacency matrices gives the number of links in the monopartite
networks: $2L_{C}=\sum_{c,c^{\prime}}(A_{C})_{cc^{\prime}}$ and
$2L_{E}=\sum_{e,e^{\prime}}(A_{E})_{ee^{\prime}}$. The degree of compounds and
elements of the monopartite networks are respectively given by
$k_{c}=\sum_{c^{\prime}}(A_{C})_{cc^{\prime}}$ and
$k_{e}=\sum_{e^{\prime}}(A_{E})_{ee^{\prime}}$.

The degree distributions of the monopartite compounds ($k_{c}$) and elements
($k_{e}$) networks are shown in the four right panels of Fig. 2(c,d,g,h) for
both CRC and AFLOW. A striking feature of the degree distributions for
compounds is that they appear to be composed of two main modes. Further
investigation reveals that all the compounds in the high-degree bump contain
the oxygen element. We denote with a vertical black dotted line (upper panels)
the smallest degree of a compound that contains oxygen. Correspondingly, in
the elements networks oxygen has the maximum, or nearly maximum, degree, which
we denote with a vertical black dotted line (lower panels). In both datasets
we discover that compounds containing oxygen form an oxygen club. This is a
maximally interconnected community composed of compounds with a large degree.
The oxygen club is a result of oxygen’s prominence in the inorganic chemistry
and materials science datasets, as well as the rules of the network. This
particular feature is almost impossible to be captured by a specific model and
needs to be addressed by further analysis of the structure of the datasets.
The shape of the degree distributions of the elements networks appears to have
an exponential body for the CRC network and a linear decay for the AFLOW
network.

Figure 2: Degree distributions. Cumulative degree distribution, $P_{>}(d)$,
for the compounds (a,b) and elements (e,f) layers of the CRC and AFLOW
bipartite networks. Cumulative degree distribution, $P_{>}(k)$, for the
compounds (c,d) and elements (g,h) monopartite networks of CRC and AFLOW. In
(c,d) the vertical dotted lines indicate the smallest degree of a compound
with oxygen. In (e,f,g,h) the vertical dotted lines indicate the degree of the
oxygen element. The continuous green lines in (a,b) are cumulative normal
distributions, in (e) a power law $\sim x^{-1.0}$, and in (f) an exponential
curve $\sim\exp(-const*x)$. All green lines are visual guides. Black points
represent empirical data while red points are obtained from fitness model
estimates.

## Fitness Models

We model these networks by assuming that there is a hidden underlying process
where all the elements compete for prominence based on an unknown intrinsic
fitness. We discover that the abundance, $a_{e}$, of each element $e$, which
is simply the element occurrence in all compounds of a dataset shown in Fig.
1d, is an excellent quantity to consider as element fitness. Similarly, we
find that the fitness of a compound $c$ can be represented well by the number
of element species, $\ell_{c}$, it contains. We model the bipartite networks
by using normalized fitness values

$\displaystyle x^{b}_{e}=\frac{a_{e}}{\sum_{e^{\prime}}a_{e^{\prime}}},\quad
y^{b}_{c}=\frac{\ell_{c}}{\sum_{c^{\prime}}\ell_{c^{\prime}}}$ (1)

as effective parameters of a maximum-entropy fitness model [19, 20, 21, 15,
22, 23, 24].

Specifically, in our model, each pair of nodes from the two different layers
(i.e., an element $e$ and a compound $c$) is connected according to a linking
probability with a Fermionic form

$\displaystyle f(\delta,x^{b}_{e},y^{b}_{c})=\frac{\delta
x^{b}_{e}y^{b}_{c}}{1+\delta x^{b}_{e}y^{b}_{c}}$ (2)

where $\delta$ is a single tuning parameter for each network. The best fitting
value $\delta^{*}$ is extracted by matching the number of links, $L$, of the
real network with that of the model:

$\displaystyle L=\sum_{e,c}f(\delta^{*},x^{b}_{e},y^{b}_{c})$ (3)

Using $\delta^{*}$ from Eq. (3) and the normalized fitness $\\{x^{b}_{e}\\}$,
$\\{y^{b}_{c}\\}$, we calculate the expected model degrees,
$\tilde{d}_{c}=\sum_{e}f(\delta^{*},x^{b}_{e},y^{b}_{c})$ and
$\tilde{d}_{e}=\sum_{c}f(\delta^{*},x^{b}_{e},y^{b}_{c})$.

We remark that Eq. (2) derives from an entropy maximisation procedure with
degree constraints, where the fitness values replace the unspecified Lagrange
multipliers [25, 26, 23]. Hence our modelling is an effective maximum-entropy
procedure informed by a heuristic fitness ansatz, where the fitness of the
nodes generate the model degrees. The alternative route, which we do not
follow here, would be to find the values of the multipliers such that the
expected degrees match the empirical values, through e.g. likelihood
maximization.

For the monopartite networks, we follow a similar approach. We use abundance
for the fitnesses of the elements, while for the compounds we sum up the
abundances of the elements they contain, $a_{c}=\sum_{e\in c}a_{e}$. Hence

$\displaystyle x^{m}_{e}=\frac{a_{e}}{\sum_{e^{\prime}}a_{e^{\prime}}},\quad
y^{m}_{c}=\frac{a_{c}}{\sum_{c^{\prime}}a_{c^{\prime}}}.$ (4)

Links between nodes in the two monopartite networks are computed with a
linking function similar to the previous one. We have, for elements and
compounds networks respectively

$\displaystyle
f(\delta_{E},x^{m}_{e},x^{m}_{e^{\prime}})=\frac{\delta_{E}x^{m}_{e}x^{m}_{e^{\prime}}}{1+\delta_{E}x^{m}_{e}x^{m}_{e^{\prime}}},\qquad
f(\delta_{C},y^{m}_{c},y^{m}_{c^{\prime}})=\frac{\delta_{C}y^{m}_{c}y^{m}_{c^{\prime}}}{1+\delta_{C}y^{m}_{c}y^{m}_{c^{\prime}}}.$
(5)

$\delta_{E}$ and $\delta_{C}$ are still free parameters for each network,
whose values $\delta_{E}^{*}$ and $\delta_{C}^{*}$ are determined using the
number of links in the empirical networks:

$\displaystyle
2L_{E}=\sum_{e,e^{\prime}}f(\delta_{E}^{*},x^{m}_{e},x^{m}_{e^{\prime}}),\qquad
2L_{C}=\sum_{c,c^{\prime}}f(\delta_{C}^{*},y^{m}_{c},y^{m}_{c^{\prime}})$ (6)

Again, we calculate the expected degrees, $\tilde{k}_{e}$, $\tilde{k}_{c}$,
using $\delta^{*}_{E}$, $\delta^{*}_{C}$ from Eq. (6), and the normalized
fitness $\\{x^{m}_{e}\\}$, $\\{y^{m}_{c}\\}$, respectively, using
$\tilde{k}_{e}=\sum_{e^{\prime}}f(\delta_{E}^{*},x^{m}_{e},x^{m}_{e^{\prime}})$
and
$\tilde{k}_{c}=\sum_{c^{\prime}}f(\delta_{C}^{*},y^{m}_{c},y^{m}_{c^{\prime}})$.

We find very good or exceptional agreement between the real networks and the
fitness models regarding the degrees in all cases, as shown in Fig. 2. The
higher-order network measure of the degree assortativity exhibits stronger
fluctuations but is still captured on average as we report in the SI and Figs.
S3 S4.

## Community analysis

We further analyze the community structure emerging from this way of exploring
the chemical space. We use the Louvain greedy algorithm [27], a method based
on the maximization of the modularity $Q$ (a quantity related to how many
links tend to connect nodes within communities rather than nodes belonging to
different communities [28, 29, 30, 31, 32]). We identify between 3 and 5
communities in AFLOW, with $Q\approx 0.25$, and between 5 and 7 communities in
CRC, with a smaller $Q\approx 0.11$, as shown in Fig. 3. The small variability
of the results depends on the initialization of the algorithm; below we
discuss only the findings that are robust across multiple runs of the
algorithm.

Figure 3: (a,d): Community structure of the compounds monopartite networks,
for CRC and AFLOW; each community is labeled by its dominant element.
Communities are colored according to their size (blue to green representing
largest to smallest groups respectively), and the links are red. The most
connected compounds ordered by descending community size are, for CRC:
$Na_{3}PO_{3}S\cdot 12H_{2}O$, $(C_{2}H_{5}O)_{2}P(S)SNH_{4}$,
$NH_{4}SO_{3}F$, $CuCl_{2}\cdot 2NH_{4}Cl\cdot 2H_{2}O$, $NH_{4}BrO_{3}$,
$NH_{4}IO_{3}$, and for AFLOW: $C_{4}F_{12}N_{8}O_{12}S_{8}Se_{8}$
(5936f8e3995c49e8), $Co_{2}H_{48}Ni_{2}O_{40}S_{4}$ (b730426d3f44aa15),
$C_{12}Cl_{12}N_{12}O_{4}P_{4}S_{12}Sb_{4}$ (155cbacf430e2898),
$Ca_{12}F_{1}K_{1}O_{26}S_{2}Si_{4}$ (3d134e8d4260e63e). In parenthesis we
give the AFLOW unique identifiers for each compound. (b,e) Cumulative counts
of degree, $F_{>}(k)$, for each community; the vertical dotted line indicates
the smallest degree of compounds with oxygen. (c,f) Relative abundance of
element species in communities, where the most abundant elements are indicated
on the top of the plots.

As expected, there is a community of compounds of large degree, which has the
highest abundance of oxygen ($Z=8$), (Fig.3 a,d). The oxygen community has a
high overlap with the oxygen club, but they are not identical, (Fig.3 b,e).
The rest of the communities are centered around other, not necessarily
prominent elements, (Fig.3 c,f). More specifically, in the CRC network, the
second largest community is dominated by hydrogen ($Z=1$), and the third by
fluorine ($Z=9$). We notice that the three most prominent elements in the CRC
dataset overall, O, H, C ($Z=6$) (Fig.1d), and the most prominent elements of
the three largest communities, are all light elements (first row of their
groups in the periodic table) and are highly reactive. In the AFLOW network
there are communities that contain most of the oxygen ($Z=8$), sulfur ($Z=16$)
and silicon ($Z=14$). We notice that two most prominent elements in the AFLOW
dataset, O and S (Fig.1d), have their own communities, and are the first two
elements of the original group $VII$ or the newer group $16$ of the periodic
table, which are collectively called chalcogens.

## Discussion

In summary, we developed a simple but fundamentally effective way to delve
into the hidden complexity of large, aggregate, chemical datasets, and reveal
their higher-order correlations. We discovered that the connectivity of
elements to compounds follows a heterogeneous distribution with different
kinds of tails: a fat one for the CRC network and an exponential one for the
AFLOW network. We traced this significant difference to the corresponding
distributions of elemental abundance in the CRC and AFLOW datasets, as shown
in Fig.1(e,f). The connectivity analysis also revealed the special role of
oxygen in the networks as we found that it dominates all orders of correlation
amongst inorganic AFLOW and CRC, Fig.1. Therefore, we revealed that oxygen
holds a prominent position in the complexity of inorganic chemistry, beyond
simply being the most common element [33] (Oxygen has also been found to play
a central role in biochemical networks and the complexity of life [34]). A
further community analysis we performed revealed chemical knowledge of purely
topological origin. The largest communities in CRC compounds network are
dominated by light, highly reactive elements. The picture is starkly different
for AFLOW, where the most prevalent elements are somewhat heavier and less
reactive. The AFLOW compounds network is less modular, comprising more
communities, as compared to the CRC. All of the results presented in this
Report are obtained thanks to the network methodology we developed, and cannot
be derived from aggregate analyses.

In addition, we were able to formalize our findings through a maximum entropy
network approach. Our fitness models were tailored for the bipartite network
and its monopartite projections, employing a single-fitted parameter and novel
fitness values that are external to the network. Our analysis is able to
quantify self-consistently both networks of CRC and AFLOW, and reproduces
successfully their statistically different connectivity. The parsimonious
modelling methodology we developed can be applied to any bipartite network, or
to a pair of complementary networks, such as article-author networks [35],
recommendation networks [36], disease phenome-genome [37], countries-products
[38, 39, 15], food ingredients-flavors [40], social networks [16, 14],
ecological networks [17], biological and medical networks [41, 42, 18], and so
on.

Network science can benefit chemistry and materials science by reorganizing
its extensive body of knowledge through complex networks. Analyzing and
modeling chemistry networks allows us to systematize intrinsic behaviors and
emergent or occluded patterns into quantitative relations. Such informed
chemical/material graph atlases can accelerate decisions on
"synthesizability", and minimize costs for intelligent design of novel
composites with desired properties. This can be done by utilizing graphical
algorithms and network methods to complete tasks that are computationally
overwhelming or demanding to investigate as is the case when starting from raw
data, first principles, experimentally, or traditional cheminformatics [2].
Specifically, the network connectivity properties that we study here describe
the relation among existing substances, and can inform searches for
alternative or novel ones.

Our approach involves large networks of substances, different from approaches
that perform learning with neural networks on individual, modestly-sized,
molecular graphs [43, 44]. In its current form, it takes advantage solely of
the chemical composition of substances, but it can be systematically expanded
to include more material properties as node variables, such as crystal
structure, thermodynamic quantities, or mechanical properties. It can utilize
more sophisticated measures for weighted linking, e.g. the number of atoms in
common, or quantify the similarity of nodes with cosine/dot-product, for
further gains.

## Methods

Datasets creation from databases

From the AFLOW database [13] we downloaded all the compounds that also belong
to the ICSD catalogue (similarly to [11]), and have a value entry in the
following eleven properties: composition, species, density, volume atom,
pressure, valence cell iupac, spin atom, scintillation attenuation length,
energy atom, enthalpy atom, eentropy atom (electronic entropy).

We utilized the entire database of Physical Constants of Inorganic Compounds
of the CRC Handbook of Chemistry and Physics Online 102nd Edition (2021),
which is part of CHEMnetBASE [12].

Both databases may reflect the biases of their creators, historical trends in
chemistry, and/or the research interests, needs, and abilities of the
scientific and engineering community. The AFLOW database comprises solids,
while the CRC database contains compounds in all phases at standard conditions
complicating property annotation. The only implicit constraint we imposed on
the AFLOW database was the materials to be sufficiently well analysed/studied.
We presume that the CRC database was compiled with a similar criterion of most
commonly used substances. Our results are proven for our datasets, since the
global shapes of the distributions are preserved when we randomly sub-sampled
our datasets. As the databases expand and research becomes gradually more
systematic, we foresee that the potential benefit from our network analysis
would also expand past the space limited by current research.

Graph link density

The link density of the networks of elements is an increasing function of the
size of the compounds datasets considered, while the density of the networks
of compounds is independent of such size. This is due to the fact that for the
elements network, the total number of nodes (i.e. elements) is constant,
$n_{E}\sim O(100)$, while the number of links between elements increases as
more compounds are analyzed. For the compounds network, the length of
compounds. i.e. the number of elements species, is nearly constant,
$1\leq\ell_{c}\leq 8$, and as more compounds are added, both the number of
nodes (i.e. compounds) and links increase. This results in a constant link
density, which is roughly $\sim 0.20$ for the AFLOW and $\sim 0.42$ for the
CRC compounds networks (see SI, Fig. S1 ).

## Acknowledgements

We would like to thank Stefano Bonetti and Roberta Sinatra for useful
discussions. This work has been partially supported by the grant EU ’HumanE-
AI-Net’, no. 952026, and by the Italian Ministry of Foreign Affairs and
International Cooperation (“Mac2Mic”), and by MIUR PRIN-COFIN2022 grant
2022JWAF7Y.

## Author contributions statement

GT and GCa conceived the paper, designed research, performed research, wrote
the manuscript. GT conducted the analysis of the empirical data and the
numerical modeling. GCi, MK, AG, TG, contributed to research, data, and
writing of the paper.

## Additional information

Competing interests The authors declare no competing interests.

## Data Availability

The CRC dataset can be obtained from the table of Physical Constants of
Inorganic Compounds of the CRC Handbook of Chemistry and Physics Online 102nd
Edition (2021), which is part of CHEMnetBASE [12], at
https://hbcp.chemnetbase.com .
The AFLOW dataset can be obtained from the AFLOW library of crystallographic
prototypes [13], at http://www.aflowlib.org .

## References

  * [1] Kirkpatrick, P. & Ellis, C.  Chemical space.  _Nature_ 432, 823 (2004).
  * [2] Leach, A. R. & Gillet, V. J.  _An introduction to chemoinformatics_ (Springer, 2007).
  * [3] Caldarelli, G.  _Scale-Free Networks_ (Oxford University Press, 2007).
  * [4] Barabási, A.-L.  The network takeover.  _Nature Physics_ 8, 14–16 (2012).  URL https://doi.org/10.1038/nphys2188.
  * [5] Unsleber, J. P. & Reiher, M.  The exploration of chemical reaction networks.  _https://doi.org/10.1146/annurev-physchem-071119-040123_ 71, 121–142 (2020).  URL https://www.annualreviews.org/doi/abs/10.1146/annurev-physchem-071119-040123.
  * [6] Jain, A. _et al._ Commentary: The materials project: A materials genome approach to accelerating materials innovation.  _APL Materials_ 1, 011002 (2013).
  * [7] Batra, R., Song, L. & Ramprasad, R.  Emerging materials intelligence ecosystems propelled by machine learning.  _Nature Reviews Materials_ 6, 655–678 (2021).  URL https://doi.org/10.1038/s41578-020-00255-y.
  * [8] Isayev, O. _et al._ Materials cartography: Representing and mining materials space using structural and electronic fingerprints.  _Chemistry of Materials_ 27, 735–743 (2015).  URL https://doi.org/10.1021/cm503507h.
  * [9] Aykol, M. _et al._ Network analysis of synthesizable materials discovery.  _Nature Communications_ 10, 1–7 (2019).  URL https://doi.org/10.1038/s41467-019-10030-5.
  * [10] Hegde, V. I., Aykol, M., Kirklin, S. & Wolverton, C.  The phase stability network of all inorganic materials.  _Sci. Adv_ 6 (2020).  URL http://advances.sciencemag.org/.
  * [11] Veremyev, A. _et al._ Networks of materials: Construction and structural analysis.  _AIChE Journal_ 67, e17051 (2021).  URL https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.17051.  https://aiche.onlinelibrary.wiley.com/doi/pdf/10.1002/aic.17051.
  * [12] _CRC Handbook of Chemistry and Physics, Physical Constants of Inorganic Compounds_ (CRC Press, Taylor & Francis Group, an Informa Group company, 2021).  URL https://hbcp.chemnetbase.com/faces/contents/ContentsSearch.xhtml.
  * [13] AFLOW library of crystallographic prototypes.  URL http://www.aflowlib.org/.
  * [14] Holme, P., Liljeros, F., Edling, C. R. & Kim, B. J.  Network bipartivity.  _Physical Review E_ 68, 056107 (2003).  URL https://journals.aps.org/pre/abstract/10.1103/PhysRevE.68.056107.
  * [15] Saracco, F., Clemente, R. D., Gabrielli, A. & Squartini, T.  Randomizing bipartite networks: the case of the world trade web.  _Scientific Reports 2015 5:1_ 5, 1–18 (2015).  URL https://www.nature.com/articles/srep10595.
  * [16] Faust, K.  Centrality in affiliation networks.  _Social Networks_ 19, 157–191 (1997).
  * [17] Ings, T. C. _et al._ Review: Ecological networks – beyond food webs.  _Journal of Animal Ecology_ 78, 253–269 (2009).  URL https://onlinelibrary.wiley.com/doi/full/10.1111/j.1365-2656.2008.01460.xhttps://onlinelibrary.wiley.com/doi/abs/10.1111/j.1365-2656.2008.01460.xhttps://besjournals.onlinelibrary.wiley.com/doi/10.1111/j.1365-2656.2008.01460.x.
  * [18] Pavlopoulos, G. A. _et al._ Bipartite graphs in systems biology and medicine: a survey of methods and applications.  _GigaScience_ 7, giy014 (2018).  URL https://doi.org/10.1093/gigascience/giy014.
  * [19] Caldarelli, G., Capocci, A., De Los Rios, P. & Muñoz, M. a.  Scale-free networks from varying vertex intrinsic fitness.  _Physical review letters_ 89, 258702 (2002).  URL http://link.aps.org/doi/10.1103/PhysRevLett.89.258702.
  * [20] Garlaschelli, D. & Loffredo, M. I.  Fitness-dependent topological properties of the world trade web.  _Physical Review Letters_ 93, 188701 (2004).  URL https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.93.188701.
  * [21] Cimini, G., Squartini, T., Garlaschelli, D. & Gabrielli, A.  Systemic risk analysis on reconstructed economic and financial networks.  _Scientific Reports_ 5, 15758 (2015).
  * [22] Squartini, T. _et al._ Enhanced capital-asset pricing model for the reconstruction of bipartite financial networks.  _Physical Review E_ 96, 032315 (2017).  URL https://journals.aps.org/pre/abstract/10.1103/PhysRevE.96.032315.
  * [23] Cimini, G. _et al._ The statistical physics of real-world networks.  _Nature Reviews Physics 2019 1:1_ 1, 58–71 (2019).  URL https://www.nature.com/articles/s42254-018-0002-6.
  * [24] Cimini, G., Carra, A., Didomenicantonio, L. & Zaccaria, A.  Meta-validation of bipartite network projections.  _Communications Physics 2022 5:1_ 5, 1–12 (2022).  URL https://www.nature.com/articles/s42005-022-00856-9.
  * [25] Park, J. & Newman, M. E.  Origin of degree correlations in the internet and other networks.  _Physical Review E_ 68, 026112 (2003).  URL https://journals.aps.org/pre/abstract/10.1103/PhysRevE.68.026112.
  * [26] Park, J. & Newman, M. E.  Statistical mechanics of networks.  _Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics_ 70, 13 (2004).  URL https://journals.aps.org/pre/abstract/10.1103/PhysRevE.70.066117.
  * [27] Nguyen, L. V. _et al._ Fast unfolding of communities in large networks you may also like the baxter q operator of critical dense polymers alessandro nigro-a bayesian fusion model for space-time reconstruction of finely resolved velocities in turbulent flows from low resolution measurements fast unfolding of communities in large networks.  _J. Stat. Mech_ 10008 (2008).
  * [28] Newman, M. E. & Girvan, M.  Finding and evaluating community structure in networks.  _Physical Review E_ 69, 026113 (2004).  URL https://journals.aps.org/pre/abstract/10.1103/PhysRevE.69.026113.
  * [29] Newman, M. E.  Modularity and community structure in networks.  _Proceedings of the National Academy of Sciences_ 103, 8577–8582 (2006).  URL https://www.pnas.org/doi/abs/10.1073/pnas.0601602103.
  * [30] Good, B. H., Montjoye, Y. A. D. & Clauset, A.  Performance of modularity maximization in practical contexts.  _Physical Review E - Statistical, Nonlinear, and Soft Matter Physics_ 81, 046106 (2010).  URL https://journals.aps.org/pre/abstract/10.1103/PhysRevE.81.046106.
  * [31] Zhang, P. & Moore, C.  Scalable detection of statistically significant communities and hierarchies, using message passing for modularity.  _Proceedings of the National Academy of Sciences of the United States of America_ 111, 18144–18149 (2014).  URL https://www.pnas.org/doi/abs/10.1073/pnas.1409770111.
  * [32] Bongiorno, C., London, A., Miccichè, S. & Mantegna, R. N.  Core of communities in bipartite networks.  _Physical Review E_ 96, 022321 (2017).  URL https://journals.aps.org/pre/abstract/10.1103/PhysRevE.96.022321.
  * [33] Wang, H. C., Botti, S. & Marques, M. A.  Predicting stable crystalline compounds using chemical similarity.  _npj Computational Materials 2021 7:1_ 7, 1–9 (2021).  URL https://www.nature.com/articles/s41524-020-00481-6.
  * [34] Raymond, J. & Segrè, D.  The effect of oxygen on biochemical networks and the evolution of complex life.  _Science_ 311, 1764–1767 (2006).
  * [35] Newman, M. E.  The structure of scientific collaboration networks.  _Proceedings of the National Academy of Sciences_ 98, 404–409 (2001).  URL https://www.pnas.org/doi/abs/10.1073/pnas.98.2.404.
  * [36] Zhou, T., Ren, J., Medo, M. & Zhang, Y. C.  Bipartite network projection and personal recommendation.  _Physical Review E - Statistical, Nonlinear, and Soft Matter Physics_ 76, 046115 (2007).  URL https://journals.aps.org/pre/abstract/10.1103/PhysRevE.76.046115.
  * [37] Goh, K. I. _et al._ The human disease network.  _Proceedings of the National Academy of Sciences of the United States of America_ 104, 8685–8690 (2007).  URL https://www.pnas.org/doi/abs/10.1073/pnas.0701361104.
  * [38] Hidalgo, C. A. & Hausmann, R.  The building blocks of economic complexity.  _Proceedings of the National Academy of Sciences of the United States of America_ 106, 10570–10575 (2009).  URL https://www.pnas.org/doi/abs/10.1073/pnas.0900943106.
  * [39] Tacchella, A., Cristelli, M., Caldarelli, G., Gabrielli, A. & Pietronero, L.  A new metrics for countries’ fitness and products’ complexity.  _Scientific Reports 2012 2:1_ 2, 1–7 (2012).  URL https://www.nature.com/articles/srep00723.
  * [40] Ahn, Y. Y., Ahnert, S. E., Bagrow, J. P. & Barabási, A. L.  Flavor network and the principles of food pairing.  _Scientific Reports 2011 1:1_ 1, 1–7 (2011).  URL https://www.nature.com/articles/srep00196.
  * [41] Lee, D.-S. _et al._ The implications of human metabolic network topology for disease comorbidity.  _Proceedings of the National Academy of Sciences_ 105, 9880–9885 (2008).
  * [42] Zhou, X., Menche, J., Barabási, A.-L. & Sharma, A.  Human symptoms–disease network.  _Nature communications_ 5, 4212 (2014).
  * [43] Duvenaud, D. K. _et al._ Convolutional networks on graphs for learning molecular fingerprints.  In Cortes, C., Lawrence, N., Lee, D., Sugiyama, M. & Garnett, R. (eds.) _Advances in Neural Information Processing Systems_ , vol. 28 (Curran Associates, Inc., 2015).  URL https://proceedings.neurips.cc/paper_files/paper/2015/file/f9be311e65d81a9ad8150a60844bb94c-Paper.pdf.
  * [44] Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O. & Dahl, G. E.  Neural message passing for quantum chemistry.  In Precup, D. & Teh, Y. W. (eds.) _Proceedings of the 34th International Conference on Machine Learning_ , vol. 70 of _Proceedings of Machine Learning Research_ , 1263–1272 (PMLR, 2017).  URL https://proceedings.mlr.press/v70/gilmer17a.html.

##  1 Supplementary Information

###  1.1 Link density of networks

The link density of the compounds networks saturates at a finite value, while
the link density for the elements networks increases, as the number of
compounds in the analyzed dataset increases, as shown in Figure S1.

Figure S1: Link density of the compounds, $\rho_{C}=L_{C}/(n_{C}(n_{C}-1)/2)$
(blue dots), and elements, $\rho_{E}=L_{E}/(n_{E}(n_{E}-1)/2)$ (green x),
networks, for the CRC (left) and AFLOW (right) datasets, plotted against
dataset size of the number of compounds, $n_{C}$.

###  1.2 Distributions of Normalized Degree

The distributions of normalized degrees by the total number of nodes decay
faster for AFLOW than CRC. CRC has larger maximum degree relative to the
maximum degree possible, as compared to AFLOW, as we show in Figure S2.

Figure S2: Cumulative distributions of normalized degree, empirical (black
dots) and calculated (red x), for the compounds of the CRC (left) and AFLOW
(right) monopartite networks. As in Figure 2c,d of main, but the degrees are
divided by the total number of nodes/compounds, $n_{C}$ in the
dataset/network.

###  1.3 Degree assortativity of bipartite networks

In the bipartite network the nearest neighbor degrees of the opposite layer
are defined as

$\displaystyle d_{c}^{nn}=\frac{1}{d_{c}}\sum_{e=1}^{n_{E}}B_{ce}d_{e},\qquad
d_{e}^{nn}=\frac{1}{d_{e}}\sum_{c=1}^{n_{C}}B_{ce}d_{c}$ (7)

We can calculate the nearest neighbor degrees from the bipartite fitness model
with a double sum

$\displaystyle\tilde{d}_{c}^{nn}=\frac{1}{\tilde{d}_{c}}\sum_{e,c^{\prime}}f(\delta^{*},x^{b}_{e},y^{b}_{c})f(\delta^{*},x^{b}_{e},y^{b}_{c^{\prime}})$
(8)
$\displaystyle\tilde{d}_{e}^{nn}=\frac{1}{\tilde{d}_{e}}\sum_{c,e^{\prime}}f(\delta^{*},x^{b}_{e},y^{b}_{c})f(\delta^{*},x^{b}_{e^{\prime}},y^{b}_{c})$
(9)

which are plotted alongside the empirical data in Fig. S3.

The degrees of the elements that are nearest neighbors to compounds,
$d_{c}^{nn}$, are significantly spread so much so that they embrace both
assortative and dis-assortative behaviors for the whole range of compounds
degree, $d_{c}$, for both datasets, upper panels in Fig. S3. The degrees of
the compounds that are nearest neighbors to elements, $d_{e}^{nn}$, vs the
degrees of the elements, $d_{e}$, appear as much more disordered clouds with a
strong initial assortative behavior that dissipates or tends towards weak dis-
assortativty for larger $d_{e}$ values, lower panels in Fig. S3. Overall
assortativity is captured in on average for both the elements and compounds
layers in a partial success of our model of non-interacting fermions.

Figure S3: Nearest neighbor degree vs degree, empirical (black dots) and
calculated (purple x), for the compounds (top) and elements (bottom) layers of
the CRC (left) and AFLOW (right) bipartite networks. [An AFLOW dataset with a
number of compounds of $n_{C}=2016$ was used for computational purposes.]

###  1.4 Degree assortativity of monopartite networks

The nearest neighbor degree is defined as,

$\displaystyle
k_{c}^{nn}=\frac{1}{k_{c}}\sum_{c^{\prime}}(A_{C})_{cc^{\prime}}k_{c^{\prime}},\qquad
k_{e}^{nn}=\frac{1}{k_{e}}\sum_{e^{\prime}}(A_{E})_{ee^{\prime}}k_{e^{\prime}}$
(10)

In the mono-partite network the nearest neighbor degrees are calculated as,

$\displaystyle\tilde{k}_{e}^{nn}=\frac{1}{\tilde{k}_{e}}\sum_{e^{\prime},e^{\prime\prime}\neq
e}f(\delta_{E}^{*},x^{m}_{e},x^{m}_{e^{\prime}})f(\delta_{E}^{*},x^{m}_{e^{\prime}},x^{m}_{e^{\prime\prime}}),$
(11)
$\displaystyle\tilde{k}_{c}^{nn}=\frac{1}{\tilde{k}_{c}}\sum_{c^{\prime},c^{\prime\prime}\neq
c)}f(\delta_{C}^{*},y^{m}_{c},y^{m}_{c^{\prime}})f(\delta_{C}^{*},y^{m}_{c^{\prime}},y^{m}_{c^{\prime\prime}})$
(12)

The nearest neighbors degrees are plotted alongside the empirical data in
Figs. S4. We find an agreement between the real networks and the fermionic
fitness model as regards the dis-assortative behavior of the network of
elements, Fig. S4. The degree assortativity of the compounds networks is more
complicated, with a cloud for smaller degrees, and a more linear behavior at
larger degrees. The cloud exhibits stronger assortative, rather than dis-
assortative behavior, whereas the linear part is clearly dis-assortative. Our
model has qualitatively more in common with the latter, and only captures the
upper boundary of the former, in a partial success of our model of non-
interacting fermions in the compounds networks.

Figure S4: Nearest neighbor degree vs degree, empirical (black dots) and
calculated (purple x), for the compounds (top) and elements (bottom) of the
mono-partite networks of CRC (left) and AFLOW (right). [An AFLOW dataset with
a number of compounds of $n_{C}=2016$ was used for computational purposes.]

###  1.5 Modeled vs empirical degree of bipartite networks

Below we show more figures with complementary results from the bipartite
fitness modeling of the empirical networks. In Figure S5 we show the
calculated degrees versus the empirical degrees for both layers of the
bipartite networks. In Figure S6 we show the empirical and modeled degrees
versus the fitness for both layers of the bipartite networks.

Figure S5: The empirical degree vs the calculated degree (black dots), for the
compounds (top) and elements (bottom) layers of the CRC (left), and AFLOW
(right) bipartite networks. The green lines are diagonals $y=x$ as visual
guides.

###  1.6 Degrees vs fitness of bipartite networks

Figure S6: The empirical degree (black dots) and calculated degree (purple x)
vs fitness used in fermionic fitness model for the compounds (top) and
elements (bottom) layers of the CRC (left), and AFLOW (right) bipartite
networks.

###  1.7 Modeled vs empirical degree of monopartite networks

Below we show more figures with complementary results from the monopartite
fitness modeling of the empirical networks. In Figure S7 we show the
calculated degrees versus the empirical degrees for both monopartite networks.
In Figure S8 we show the empirical and modeled degrees versus the fitness for
both monopartite networks.

Figure S7: The empirical degree vs the simulated degree (black dots) for the
compounds (top) and elements (bottom) mono-partite networks of CRC (left) and
AFLOW (right). The green lines are diagonals $y=x$ as visual guides.

###  1.8 Degrees vs fitness of monopartite networks

Figure S8: The empirical degree (black dots) and simulated degree (purple x)
vs fitness for the compounds (top) and elements (bottom) mono-partite networks
of CRC (left) and AFLOW (right).