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 theory3,4.

In the past, chemical reaction networks have been presented for small numbers of reactants5, 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 project6,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 cartography8. One network approach has been based on the representation of materials phase diagrams9,10. A different approach was to analyze a set of materials as a network, according to the cross-correlation of the electronic density of states11. 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 compounds12, and the second, AFLOW, is based on inorganic material compounds13 (see “Methods”). Each of the datasets contains \(n_C\) compounds and \(n_E\) elements; the specific values are shown in Table 1.

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.
Full size table

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 layers14,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
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.

Full size image

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'}=1\) if \(\sum _e B_{ce}B_{c'e}>0\), \((A_E)_{ee'}=1\) if \(\sum _c B_{ce}B_{ce'}>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'} (A_C)_{cc'}\) and \(2L_E = \sum _{e,e'} (A_E)_{ee'}\). The degree of compounds and elements of the monopartite networks are respectively given by \(k_c = \sum _{c'} (A_C)_{cc'}\) and \(k_e = \sum _{e'} (A_E)_{ee'}\).

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

Full size image

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

$$\begin{aligned} x^b_{e} = \frac{ a_{e} }{ \sum _{e'} a_{e'} }, \quad y^b_{c} = \frac{ \ell _{c} }{ \sum _{c'} \ell _{c'} } \end{aligned}$$
(1)

as effective parameters of a maximum-entropy fitness model15,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

$$\begin{aligned} f(\delta ,x^b_e,y^b_c) = \frac{ \delta x^b_e y^b_c }{1 + \delta x^b_e y^b_c} \end{aligned}$$
(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:

$$\begin{aligned} L = \sum _{e,c} f(\delta ^{*},x^b_e,y^b_c) \end{aligned}$$
(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 multipliers23,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} = \sum _{e \in c} a_{e}\). Hence

$$\begin{aligned} x^m_e=\frac{ a_{e} }{ \sum _{e'} a_{e'} },\ \ y^m_{c} = \frac{ a_{c}}{ \sum _{c'} a_{c'} }. \end{aligned}$$
(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

$$\begin{aligned} f(\delta _E,x^m_e,x^m_{e'}) = \frac{ \delta _E x^m_e x^m_{e'} }{1 + \delta _E x^m_e x^m_{e'}}, \qquad f(\delta _C,y^m_c,y^m_{c'}) = \frac{ \delta _C y^m_c y^m_{c'} }{1 + \delta _C y^m_c y^m_{c'}}. \end{aligned}$$
(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:

$$\begin{aligned} 2L_E = \sum _{e,e'} f(\delta _E^{*},x^{m}_e,x^{m}_{e'}), \qquad 2L_C = \sum _{c,c'} f(\delta _C^{*},y^{m}_c,y^{m}_{c'}) \end{aligned}$$
(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'} f(\delta _E^{*},x^{m}_e,x^{m}_{e'})\) and \({\tilde{k}}_c = \sum _{c'} f(\delta _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 and S4.

Community analysis

We further analyze the community structure emerging from this way of exploring the chemical space. We use the Louvain greedy algorithm27, 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 communities28,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
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.

Full size image

As expected, there is a community of compounds of large degree, which has the highest abundance of oxygen (\(Z=8\)), (Fig. 3a,d). The oxygen community has a high overlap with the oxygen club, but they are not identical, (Fig. 3b,e). The rest of the communities are centered around other, not necessarily prominent elements, (Fig. 3c,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. 1e,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 element33 (Oxygen has also been found to play a central role in biochemical networks and the complexity of life34). 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 networks35, recommendation networks36, disease phenome-genome37, countries-products15,38,39, food ingredients-flavors40, social networks14,16, ecological networks17, biological and medical networks18,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 cheminformatics2. 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 graphs43,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 database13 we downloaded all the compounds that also belong to the ICSD catalogue (similarly to11), 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 CHEMnetBASE12.

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 \le \ell _c \le 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).