Network topology mapping of chemical compounds space

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  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 × 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 = ∑ c,e B ce .The degree of a node is the sum of its incident links: d c = ∑ e B ce and d e = ∑ 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.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: and zero otherwise.Summing all non-zero entries in the adjacency matrices gives the number of links in the monopartite networks: 2L C = ∑ c,c ′ (A C ) cc ′ and 2L E = ∑ e,e ′ (A E ) ee ′ .The degree of compounds and elements of the monopartite networks are respectively given by k 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.

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, ℓ c , it contains.We model the bipartite networks by using normalized fitness values as effective parameters of a maximum-entropy fitness model 15,[19][20][21][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 where δ is a single tuning parameter for each network.The best fitting value δ * is extracted by matching the number of links, L, of the real network with that of the model: ).We remark that Eq. ( 2) derives from an entropy maximisation procedure with degree constraints, where the fitness values replace the unspecified Lagrange multipliers 23,25,26 .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 = ∑ e∈c a e .Hence 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 δ E and δ C are still free parameters for each network, whose values δ * E and δ * C are determined using the number of links in the empirical networks: x m e ′ ), Again, we calculate the expected degrees, ke , kc , using δ * E , δ * C from Eq. ( 6), and the normalized fitness {x m e }, {y m c }, respectively, using ke = ∑ e ′ f (δ * E , x m e , x m e ′ ) and kc = ∑ c ′ f (δ * C , y m c , y m c ′ ).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 ≈ 0.25, and between 5 and 7 communities in CRC, with a smaller Q ≈ 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.
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 V II 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 15,38,39 , food ingredients-flavors 40 , social networks 14,16 , ecological networks 17 , biological and medical networks 18,41,42 , 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.

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 ∼ 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 ≤ ℓ c ≤ 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 ∼ 0.20 for the AFLOW and ∼ 0.42 for the CRC compounds networks (see SI, Fig. S1 ).
1 Supplementary Information

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.

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.

Degree assortativity of bipartite networks
In the bipartite network the nearest neighbor degrees of the opposite layer are defined as We can calculate the nearest neighbor degrees from the bipartite fitness model with a double sum which are plotted alongside the empirical data in Fig. S3.The degrees of the elements that are nearest neighbors to compounds, d nn c , 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 nn e , 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.

Degree assortativity of monopartite networks
The nearest neighbor degree is defined as, In the mono-partite network the nearest neighbor degrees are calculated as, 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.

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.

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 1 .
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.

Figure 2 .
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 ∼ x −1.0 , and in (f) an exponential curve ∼ exp(−const * x).All green lines are visual guides.Black points represent empirical data while red points are obtained from fitness model estimates.
Using δ * from Eq. (3) and the normalized fitness {x b e }, {y b c }, we calculate the expected model degrees, dc = ∑ e f (δ * , x b e , y b c ) and de = ∑ c f (δ * , x b e , y b c

Figure 3 .
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 • 12H 2 O, (C 2 H 5 O) 2 P(S)SNH 4 , NH 4 SO 3 F, CuCl 2 • 2NH 4 Cl • 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.

f
(δ * , x b e , y b c ) f (δ * , x b e ′ , y b c )

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.

11 / 17 Figure S3 . 17 Figure S4 .
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.]

Figure S5 .Figure S6 .
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.

Figure S7 .Figure S8 .
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.

Table 1 .
(see Methods).Each of the datasets contains n C compounds and n E elements; the specific values are shown in Table1.Information about the datasets used to build all the networks.n C and n E are the number of compounds and elements.