A general hydrogen-bond connectivity descriptor based on graph theory

Abstract

Hydrogen bonds play a central role in determining the structure and behaviour of liquids, particularly water, whose anomalous properties arise from its extended and highly heterogeneous hydrogen-bond network. Here, we introduce a graph-theoretical framework related to the Node Total Communicability (NTC) that enables a systematic, molecule-resolved description of hydrogen-bond networks from Molecular Dynamics (MD) simulations. In contrast to earlier related approaches with the NTC, the hydrogen-bond network is explicitly mapped onto a directed graph, allowing the asymmetric nature of hydrogen bonding to be retained. The method captures both local and longer-range connectivity, providing a unified metric to probe the structural organization of molecules. As a representative test case, we apply this approach to water and aqueous salt solutions, demonstrating how the presence of ions modifies hydrogen-bond connectivity. From this perspective, the analysis shows that ionic effects are confined to the first hydration shell, while the hydrogen-bond network beyond remains essentially unperturbed.

---

Introduction

Water is the universal solvent and the basis of countless physical, chemical, and biological processes. In the liquid state at room temperature, water forms a dynamic three-dimensional hydrogen-bond (HB) network in which some molecules adopt a nearly tetrahedral configuration, while others exhibit more distorted arrangements due to thermal fluctuations.¹ The structural organization of this network controls liquid water's thermodynamic and dynamic behavior and the competition between the tetrahedral and distorted local configurations has been identified as the source of the anomalies of liquid water, such as its high heat capacity, density maximum, and surface tension.^2,3 The local structure of liquid water has been the focus of decades of experimental and computational studies aiming at connecting the microscopic structure to macroscopic behavior.^2,4,5 From a computational standpoint, several structural order parameters, such as LSI (local structure index),⁶ q,⁷ d5,⁸ ζ,⁹ V4,¹⁰ V_4S^11,12 θ_avg,¹³ NTC (node total communicability),¹⁴ H_c,¹⁵ and Ψ¹⁶ have been introduced to characterize water's local topology from different features. The propensity of water molecules to organize into small structural units was also noted, and the formation of cyclic clusters by tetramers, pentamers, and hexamers has been investigated.^17–19

The addition of solutes introduces strong perturbations into water's HB network. Their presence changes viscosity, diffusion, and solvation structure and affects the stability of biomolecules.^20–22 Extensive work has shown that HBs in the first solvation shell are significantly altered, while the extent of the disturbance beyond that shell remains debated, especially for divalent ions and specific ion pairs.^23–25 Ion-specific effects are often summarized in the Hofmeister series, which ranks ions by their influence on water and biological systems. From a structural point of view, ions are usually divided into kosmotropes (“structure makers”), which strengthen hydrogen bonding, and chaotropes (“structure breakers”), which weaken it. However, at least for anions, this description has proven to be only partial. A very recent study revealed that two distinct mechanisms characterize anion–water interactions: soft anions (such as Cl⁻, Br⁻, and I⁻) modulate HB strength through partial transfer of electronic charges, resulting in weaker H-bonds; however, this effect is partially compensated by an increase in orientationally correlated H-bonds. In contrast, hard anions (F⁻, ClO₄⁻, and SO₄²⁻) perturb the network primarily through electrostriction.²⁶ An analogous dichotomy for cation–water interactions has not yet been reported. From a broader perspective, understanding the effects of the addition of ions on water's HB topology requires an analysis linking the local solvation structure to the collective network behavior. To this aim, graph theory based approaches appear as promising tools, because, by representing water molecules as nodes and bonds or interactions as edges, they allow to take into account the entire water–water HB network.²⁷ These kinds of approaches have been used to quantify the connectivity of water molecules, to identify the recurrence of specific structural motifs, and to characterize how solutes and nanoconfinement modify clustering behavior.^28–31

In our previous works,^14,32,33 we used a graph theory based structural descriptor, the node total communicability (NTC), to investigate the low- and high-density liquid states of water. Within that approach, we defined a graph, tailored to capture connectivity changes related to the density of the system, taking into account both the local density and the medium- to long-range structural organization of water molecules. Expanding on this concept, we now introduce a metric, derived from the same network-based information underlying the NTC, that allows a full explicit characterization of the whole water–water HB network, still retaining the information about the medium-to-long range environment. This metric, the HB score (S_HB), is applied here at first to bulk water, both in its crystalline state and under ambient conditions, and is then used to characterize how the whole HB network of water reorganizes in the presence of solutes such as NaCl. The binary solution is investigated across a range of concentrations to demonstrate the descriptive power of the approach. Since aqueous NaCl is a widely studied system, it serves as a controlled test case to validate the method and illustrate the richness of the information it provides. In particular, the HB score allows us to identify ion-specific effects on connectivity, resolve how water molecules behave and orient around different ions, and quantify the extent of the perturbation at each geodesic distance, providing a level of detail that is not readily accessible with standard post-processing tools.

Methods

Molecular dynamics simulations

All MD simulations were performed with the 5.1.2 version of the GROMACS software³⁴ using a cubic simulation box containing 4440 water molecules. The simulation protocol followed the same setup and parameters reported in Blazquez et al.³⁵ Temperature and pressure were kept constant by using the Nosé–Hoover temperature coupling³⁶ and the Parrinello–Rahman barostat with 2 ps relaxation times.³⁷ Periodic boundary conditions were used, long range electrostatic interactions were treated with the particle mesh Ewald method³⁸ with a real space cutoff of 1 nm; a cut-off radius of 1 nm was employed for short range interactions. All bonds were constrained using the LINCS algorithm³⁹ along with a 2 fs time step. Six simulations were performed with the same number of water molecules. One with only water, three including 32, 160, and 320 NaCl ion pairs (corresponding to 0 molal, 0.4 molal, 2 molal, and 4 molal salt concentrations, respectively) and two single-ion simulations (0.0125 molal). For each system a 15 ns NPT simulation at 1 bar and 300 K was run. The average density from the last 10 ns was used to initialize a 60 ns NVT simulation. The analyses reported in this work are performed over the last 50 ns of each NVT run. Further details of the density and water diffusion coefficients corrected for the size effect (Yeh-Hummer correction),⁴⁰ are provided in Table S1 of the SI.

Graph theory elements and graph construction

In this section, we present some concepts of graph theory,⁴¹ explaining how we construct the graphs given an MD trajectory. We refer to our previous works^14,32 for a more detailed description. A graph G = (V, E) consists of two sets: a set of nodes V = {v₁, v₂, …, v_N} and a set of edges E ⊆ V × V. If each edge in the graph has an orientation, the graph is called directed; otherwise, it is an undirected graph. In this work, we consider unweighted graphs (the weight associated with each edge is set to 1) without self-loops, i.e., edges of the form (v_i, v_i). Given a graph G, its associated adjacency matrix is a square matrix of dimension N × N such that the entry in position (i, j) is defined as

It is easy to prove that in an unweighted graph G, the number of walks of length k = 1, 2, … between two nodes v_i and v_j is equal to the entry (i, j) in the matrix A^k (the k-th power of the associated adjacency matrix).

In our previous works,^14,32,33 we constructed undirected graphs for bulk water, considering the oxygen atoms as the nodes of our graphs and defining an edge between two nodes as existing if the distance between two oxygen atoms is <0.35 nm (see Faccio et al.¹⁴ for explanations on the choice of this threshold). For the scope of the present work, we consider an edge between two nodes as existing if the two oxygen atoms are connected by a hydrogen bond (Fig. 1a). Therefore, the existence of a connection is further restricted by the hydrogen–donor–acceptor (H–D–A) angle, which must be <30°. With this definition of an edge, we can construct both undirected and directed graphs. For the latter, the edge is directed from donor to acceptor. For directed graphs, where edges have a defined orientation, the HB network of each water molecule can be represented as the sum of two contributions. From a graph-theory perspective, these correspond to the broadcaster and receiver roles (Fig. 1b), while from a chemical view-point, they reflect the donor and acceptor character of the molecule.

Fig. 1 (a) Hydrogen bonds for a reference molecule. (b) The same hydrogen bonds highlighted in the broadcaster and receiver. (c) Distribution of HB-NTC values for pure water at 300 K (red). The corresponding value for ice Ih is included as a reference (black dashed).

The geometric definition of hydrogen bonding we employ here is widely adopted in the literature.⁴² An alternative definition based on an energetic criterion could also be employed. A discussion of the sensitivity of our results to this choice, including a comparison with an energetic criterion based on the two-dimensional free-energy landscape along the H-bond distance and angle coordinates, is provided in the SI (see Fig. S1 and S2 and the corresponding text).

To calculate the distances and angles between the atoms, we use the minimum image convention. In this way, we take into account periodic boundary conditions. In order to focus on water–water connectivity, Na⁺ and Cl⁻ ions are excluded from the network, and only oxygen atoms of water are treated as nodes in all graphs considered in this work. Including Na⁺ and Cl⁻ in the graph would not pose a technical limitation. However, it would lead to a network that is no longer purely based on hydrogen bonding. In particular, Cl⁻ can act as an H-bond acceptor, whereas Na⁺ cannot act as a donor, thereby introducing an intrinsic asymmetry in the H-bond definition.

As a consequence of excluding water–ion interactions, water molecules coordinated to the ions appear less connected in the graph, since part of their water–water H-bonds are replaced by ion coordination. This implies that, in addition to the genuine structural modification of the HB network induced by the presence of the ion, an additional purely geometric excluded-volume effect also contributes to the observed differences in connectivity. To estimate the effect of the excluded volume of each ion, we compare a system containing one ion at the center of the box (50 ns-long single-ion simulations in the NVT ensemble, i.e. 0.0125 molal) and the corresponding pure water simulation, corrected for the excluded volume through the following procedure. For each pure water trajectory frame, we remove from the graph all water molecules whose oxygen–box-center distance is smaller than a threshold radius, defined as the first non-zero value of the oxygen–ion radial distribution function (1.86 Å for Na⁺ and 2.55 Å for Cl⁻). This effectively places a “fictitious” ion, a hard sphere of the same radius as the real ion's first coordination shell onset, at the box center, without introducing explicit ion–water interactions. This allows us to isolate the geometric excluded-volume effect from the direct perturbation of the HB network induced by ion coordination.

Node total communicability and HB score calculation

In our previous works,^14,32 we introduced an order parameter, the node total communicability (NTC), borrowed from graph theory.⁴³ Given a node v_i, its NTC is defined aswhere β is a positive parameter and 1 is the vector of all ones. Consequently, the NTC of a node v_i takes into account all the walks between v_i and the other nodes in the graph, but, through the factor , it gives less weight to the contribution of longer walks. In our experiments, we choose β = 1 (refer to our previous works^14,32 for the details of this choice). In our previous works, we considered two oxygen atoms to be connected if their distance was <0.35 nm, thereby treating interstitial water molecules, one of the typical structural features of ambient water, as connected nodes. Since our focus here is the characterization of the water HB network, we instead compute the NTC by considering two oxygen atoms as connected only if a hydrogen bond exists between them (see above). The distribution of the NTC defined in this way (hereafter denoted as HB-NTC) is reported for both ice Ih and ambient water in Fig. 1c.

For the regular lattice of Ice Ih, the NTC assumes a single value (i.e., e⁴ ≈ 54.6), consistent with a 4-regular graph.⁴⁴ For ambient water, the NTC distribution is broader, reflecting the enhanced topological and structural heterogeneity of the HB network. Compared with Ice Ih, the main peak is shifted toward lower values, consistent with the reduced average HB connectivity and the increased presence of network defects in the liquid phase. Secondary peaks and shoulders suggest the coexistence of distinct local HB environments within the network.

As mentioned before, the HB-NTC accounts for the global connectivity of each node through the summation over k. However, by truncating this summation at a determined value of k and using a directed graph in which edges correspond to hydrogen bonds, it is possible to gather information about the specific HB pattern of a node. This approach can be applied both by using the standard NTC formula (eqn (2)) to evaluate the magnitude of the connectivity of each molecule in terms of hydrogen bonds up to the geodesic distance k and by using an unweighted version of the NTC (i.e., without the factor ). With this second approach, the resulting data are more directly interpretable, as they provide a count of the number of HBs for each value of k (see Results and discussion section for the correspondence between k and the physical distance). With this approach we define the HB score, S_HB (v_i), that can be viewed as a k-resolved, unweighted counterpart of the HB-NTC, from which the HB-NTC can be directly constructed. S_HB (v_i) of node v_i can thus be defined as a graph-theory–based metric quantifying the hydrogen-bond connectivity of a molecule, with contributions resolved by k. Practically speaking, S_HB (v_i) is the string concatenation of the number of HBs n_HB (v_i,k) obtained for each k:

S_HB (v_i) = ||n_HB (v_i,k = 1)||n_HB (v_i,k = 2)||n_HB (v_i,k = 3)|| (3)

It has to be noted that, if within this count the number of HBs as a donor and as an acceptor has to be assessed (corresponding to the broadcaster and receiver nature of each node), the adjacency matrix derived from the geometric hydrogen-bond criteria cannot be used directly. A few adjustments are necessary to ensure coherent results; further details are provided in Section S3 of the SI.

Results and discussion

The aim of the graph-theory approach we propose is to provide a full characterization of the water–water HB pattern of molecules in a system (either as pure or as a solution), while also separating the contributions associated with distinct hydration shells. As discussed in the Methods section, this can be achieved by computing the HB score of node (v_i), S_HB (v_i), using a specific value of k. It is therefore necessary to evaluate the correlation between the geodesic distance k and the corresponding Euclidean distance in three-dimensional space. To this end, for each water molecule in a bulk water simulation at ambient pressure and temperature, we consider its neighbours for each k up to k = 5 and compute their physical distance from the reference molecule. The results of this analysis are presented in Fig. 2, which shows the oxygen–oxygen radial distribution function (RDF) alongside the distribution of the Euclidean distance for each geodesic distance k. As expected, k = 1 corresponds to the first hydration shell, while k = 2 is primarily localized in the second hydration shell. For k = 3, the chemical neighbours are spread into the first three hydration shells. At k = 4 and k = 5, the neighbours mostly occupy regions where long-range order is increasingly lost. k = 3 appears thus as a good compromise between capturing information beyond the local hydrogen bond network and avoiding noise arising from regions lacking long-range order. Therefore, all data presented in this work are obtained up to k = 3, including S_HB contributions.

Fig. 2 Oxygen–oxygen radial distribution function for the pure water simulation, showing the distribution of distances of molecules at geodesic distances up to k = 5.

The same analysis is also carried out for the ionic solutions, and the results do not show significant differences, as illustrated in Fig. S5 of the SI. It is also worth noting that a very small peak is observed in the first hydration shell for k = 4 and k = 5. This feature arises from the propensity of the hydrogen-bond network of water to form cyclic arrangements, so that molecules separated by four or five geodesic steps can still be located at short Euclidean distances as part of five- and six-membered hydrogen-bonded rings. In this context, the adjacency matrix representation of the network is particularly convenient, since information on rings of size n can be directly extracted from the trace of Aⁿ.

Once k = 3 is set as a convenient threshold, the HB score of each molecule can be determined according to eqn (3). While a single digit suffices for the k = 1 string, two digits are required for both the k = 2 and k = 3 strings. Consequently, a hydrogen bond configuration up to k = 3 can be represented using a five digit HB score. For example, each water molecule in a perfect ice Ih lattice can be described by an HB score of 41236 (i.e. 4 hydrogen bonds at k = 1, 12 at k = 2 and 36 at k = 3, Fig. 3a). Panel b of Fig. 3 reports the distribution of S_HB (v_i) for pure water at 300 K. The peaks in panel b report the different numbers of HB of water molecules in the system (i.e., they are the distribution of the values obtained at k = 1). However, each of these peaks is composed of multiple sub-peaks that arise from additional structural distinctions captured at k = 2, as shown in panel c. Analogously, these k = 2 sub-peaks further subdivide into even finer features reporting on the HB pattern at k = 3 (panel d). The value corresponding to ice Ih is also shown in each panel as a reference. The comparison clearly shows how the HB network of ambient water differs from the perfectly tetrahedral one of ice. In particular, the data for k = 1 show the competition between different local environments in ambient water that has been identified as the source of anomalies in liquid water. As a matter of fact, while the most pronounced peak is the same as Ice Ih (4 HBs, tetrahedral structure), a very relevant population of molecules with 3 HBs (distorted structure) is also present. The data for k = 2 and k = 3 show that this difference extends to longer distances: the peak of Ice Ih for k = 2,3 is only marginally populated in ambient water. This reflects the expected increase in hydrogen bond network distortion compared to the ordered structure of ice. These results also show that, with this approach, we can exactly quantify the higher degree of disorder in ambient water and its origin. In fact, panels c and d show that even water molecules involved in four hydrogen bonds feature a less structured second and third hydration shell. It is also worth highlighting that this information can be obtained with a single, quick calculation exploiting the high efficiency of graph-theory optimized algorithms. In contrast, information about the long range HB pattern of a molecule is not readily accessible with standard post-processing tools typically used to analyze MD simulations.

Fig. 3 (a) Example snapshot of the ice Ih structure, with neighbours from a reference molecule coloured according to their geodesic distance. (b)–(d): Distribution of S_HB (v_i) for bulk water at 300 K for k = 1, k = 2 and k = 3, respectively. S_HB of ice Ih is also shown as a dashed black line for comparison.

We then move to the investigation of aqueous NaCl solutions at different concentrations. In our previous works, we showed that the NTC is a useful order parameter to distinguish the two local structures of liquid water (ordered vs distorted). In that case, the lower NTC values observed in the low density (ordered) state of liquid water can be attributed to the lower connectivity of that state. Also for the NaCl solutions investigated here, the NTC shows a progressive decrease in the number of hydrogen bonds as the salt concentration increases, resulting in a lower connectivity in the graph (Fig. S6a in the SI). This trend becomes considerably more pronounced for the HB-NTC distribution (Fig. S6b in the SI), showing a greater sensitivity of the HB-based NTC in discriminating between solutions at different salt concentrations.

The same figure also shows that the HB-NTC distributions feature multiple peaks arising from the different HB configurations of the water molecules in the system. However, the interpretation of these peaks in terms of HB pattern is not straightforward, mainly because distinct hydrogen-bond configurations can produce the same NTC value. This ambiguity can be avoided by using the HB score to characterize the salt solutions. To characterize the concentration dependence of the hydrogen-bond pattern of water molecules surrounding the ions, we first analyze the behavior within the first solvation shell, whose extent is defined by the first minimum of the ion–oxygen radial distribution functions (Fig. S7 of the SI). For water molecules in the first shell, the connectivity at k = 3 decreases with increasing concentration. Panels (a) and (b) of Fig. 4 report the n_HB/n_HB0 ratio for sodium and chloride, respectively, where n_HB0 is defined as the average number of HBs at a given k of pure water. A systematic difference between the two ions is observed: water molecules around chloride are always more connected than those around sodium, reflecting different local structuring induced by the two ions (vide infra). As k decreases, the n_HB/n_HB0 ratio becomes more independent of concentration for both ions.

Fig. 4 (a) and (b): n_HB/n_HB0 ratio for molecules in the first shell of ions (Na⁺ and Cl⁻, respectively) for k = 1 (red curves), k = 2 (blue curves) and k = 3 (green curves) as a function of molal concentration. The number on the curves represents their slope. (c): n_HB/n_HB0 ratio for molecules in the first shell of ions (yellow, Na⁺; purple, Cl⁻) and all of the molecules outside the first shell of the ions (dark red) for k = 1, the n_HB0 value refers to that of pure water in all panels.

In Fig. 4, panel c, we compare, at k = 1, the n_HB of water molecules within the ionic solvation shells to that of water molecules outside the shells, using pure water as a reference (n_HB0). The analysis shows that the effect is almost constant across concentrations and that the n_HB values outside the solvation shells are indistinguishable from those of pure water, suggesting that the perturbation does not extend beyond the first solvation shell, in line with previous findings.²⁴

To investigate more thoroughly the possibility that ions perturb the water network beyond the first solvation shell, we perform single-ion simulations and compute the average number of hydrogen bonds (n_HB) at geodesic distances k = 1, k = 2, and k = 3 for water molecules in the first and second solvation shells of each ion, as well as for bulk water molecules. To disentangle the effects arising from the excluded volume of the ion from those associated with direct ion coordination, the same analysis was also performed on the “fictitious” ion configurations (see the Methods section). The results are reported in Fig. 5. For water molecules in the first solvation shell of ions, the comparison with the fictitious ion system reveals a clear difference between the two ions: for Na⁺ (panel c), the real ion produces a significantly lower n_HB than the fictitious one, indicating a strong direct coordination effect that goes beyond a simple excluded-volume contribution; for Cl⁻ (panel d), the real and fictitious curves are identical, showing that the small reduction in n_HB observed in the first shell is entirely accounted for by the excluded volume of the ion. This result shows that Cl⁻ integrates into the hydrogen-bond network in a manner that closely resembles water itself, consistent with its ability to adopt hydrogen-bond-like configurations with neighbouring water molecules. In the second solvation shell (panels e and f), real and fictitious ions yield indistinguishable n_HB values for both Na⁺ and Cl⁻ at all geodesic distances. These results demonstrate that the differences between bulk water and the solution for water molecules in the second solvation shell originate from the geometrical truncation of available network pathways due to the excluded volume of the ion, rather than from any propagation of the ion-induced perturbation beyond the first solvation shell.

Fig. 5 Representative configuration of the single-ion simulation (a) and of the pure water simulation used to construct the fictitious-ion reference (b). The grey sphere highlights the threshold radius defined as the first non-zero value of the oxygen-ion radial distribution function (1.86 Å for Na⁺ and 2.55 Å for Cl⁻). In the fictitious ion construction, all water molecules in the pure water simulation whose oxygen falls within the sphere centered at the box center (black dot) are removed from the graph; the connections eliminated by this procedure are indicated by red crosses. Panels (c)–(f) Average number of HBs for k = 1, k = 2 and k = 3 for the water molecules found in the first and second shells of the real and fictitious ions, bulk water values are also shown for reference.

We then investigate the structural origin of the different behavior of water molecules in the first solvation shell of the two ions. To determine the relative orientation of these water molecules and the ions, we compute the angle θ, defined as the angle between the ion–oxygen vector and the HOH angle bisector, as illustrated in Fig. 6, panels (a) and (b). The distribution of these angles was analyzed separately for molecules forming 1, 2, 3, or 4 hydrogen bonds.

Fig. 6 Distribution of the angle theta, for molecules found in the first hydration shell of sodium (a) and chloride (b), divided in the number of HBs in which molecules are involved. (c) Example HB configurations for water molecules in the first solvation shell of the ions engaging in 2 and 3 HBs. (d) Distribution of the number of hydrogen bonds at geodesic distance k = 1 of water molecules in the first solvation shell of sodium (black) and chloride (red). The inset represents the difference between the two distributions (Na⁺–Cl⁻). In particular, the yellow region marks where sodium contributes more than chloride, while the purple region indicates the opposite.

In an ideal configuration, sodium prefers to coordinate with the oxygen atom facing the ion and the hydrogens pointing outward, corresponding to a theta value of ≈150°–180°. Panel (a) shows that, to form a larger number of hydrogen bonds, water molecules around sodium adopt increasingly distorted conformations to expose the oxygen to potential hydrogen-bond donors. In contrast, water molecules around chloride orient themselves with a hydrogen pointing toward the ion, in a hydrogen-bond-like fashion, yielding a theta value of ≈50°–60°, resembling the configuration of an HB, as widely reported in the literature.^28,45 Panel (b) demonstrates that this orientation is largely preserved regardless of the number of hydrogen bonds in which the molecule participates. Panel (c) of Fig. 6 illustrates how water molecules around sodium adopt distorted orientations to accommodate more than two hydrogen bonds, in contrast to the more regular arrangement observed around chloride. From the HB distribution shown in panel d, it can be seen that, for both ions, the most frequent number of hydrogen bonds is three. However, relative to each other, water molecules around sodium more often feature one or two hydrogen bonds, while those around chloride more frequently engage in three or four. This reflects the fact that maintaining an additional hydrogen bond requires a partial reorientation for sodium. Fig. S8 of the SI, also shows that in order to engage in a higher number of HBs, water molecules in the first hydration shell of Na⁺ are found at a larger distance from the ions, partially losing the coordination with the ions.

Resolving the number of HBs for water molecules around the two ions in terms of broadcaster and receiver contributions (Fig. 7), the expected organization of water around NaCl can be observed. Fig. 7 shows that molecules around chloride and sodium display opposite behaviors in terms of broadcaster and receiver roles, consistent with the preferred orientation of water molecules, with the oxygen atom pointing toward Na⁺ and the hydrogen atoms toward Cl⁻.

Fig. 7 Broadcaster and receiver histograms, at a given hydrogen bond number, of water molecules found in the first solvation shell of ions in the lowest molality simulation. Each column represents the most sampled configurations for a given number of total HBs. For each column of the histogram, a reference snapshot is provided, with the reference molecule highlighted and the relative neighbours rendered transparently.

When examining the hydrogen bond distribution for k = 2 (reported in Fig. S9 and S10 of the SI), we define n₁ as the number of hydrogen bonds at k = 1. For n₁ values less than or equal to two, the sodium and chloride distributions follow the same trend. Noticeable differences appear only when n₁ is greater than two. This is because both sodium and chloride can engage in two HBs without relevant distortions. In Fig. 8, we show the distribution of the number of HBs for k = 2 for n₁ = 3 and n₁ = 4, respectively, it can be seen that for the same n₁ value, chloride is able to reach a slightly higher connectivity (insets in Fig. 8).

Fig. 8 Top panel: hydrogen-bond distribution at k = 2 for molecules with n₁ = 3, relative to the water molecules in the first solvation shell of the ions. Bottom panel: hydrogen-bond distribution at k = 2 for molecules with n₁ = 4, also relative to the first-shell molecules. The insets show the difference between the sodium and chloride distributions, with the yellow region marking higher contribution of sodium and the purple region marking the opposite. The x-axis uses a three digit HB score.

Moving to k = 3, if we consider a fixed n₁/n₂ couple (i.e., a given number of hydrogen bonds at k = 2 for a given n₁), no relevant difference between sodium and chloride emerges (see Fig. S11 in the SI for example). To understand this observation we explicitly inspected the ion–oxygen distances for molecules at k = 1, k = 2, and k = 3 geodesic distances from the molecules in the first solvation shell of ions, as shown in Fig. 9. For molecules at k = 1, the ion–oxygen distances fall squarely within the second solvation shell with a small shoulder in the first hydration shell, indicating that there's only a small percentage (1.9% for both Na⁺ and Cl⁻) of direct hydrogen bonds in the first hydration shell of ion distribution functions. The difference in this shoulder, however, increases in percentage when examining differences at k = 2 with 3.0% for sodium and 4.1% for chloride. This suggests a scenario in which some HB paths originating from the first hydration shell return to the same shell at k = 2 through an intermediate water molecule in the second solvation shell, which bridges two water molecules belonging to the first shell (Fig. 9). For sodium, this triangular-like configuration requires at least one of the two first-shell molecules to be distorted with respect to the ion, whereas the same configuration can be easily attained by chloride.

Fig. 9 Ion–oxygen radial distribution functions (RDFs), together with ion–oxygen distance distributions for water molecules at geodesic distances k = 1, k = 2, and k = 3 from molecules in the first solvation shell of the ions. For each RDF, molecules in the first solvation shell are identified, and the subsets at geodesic distances k = 1, 2, and 3 from these molecules are selected. The distributions report the distances between the ions and the selected water molecules. Both insets show representative snapshots of molecular orientations in the first hydration shell, with reference molecules in yellow, k = 1 molecules in red, k = 2 molecules in blue and k = 3 molecules in green.

On the other hand, at k = 3, the percentages in the first shell become similar again (4.2% for Na⁺ and 4.3% for Cl⁻), as the hydrogen-bond pattern forms a square-like arrangement (Fig. 9), resulting in a less strained structure where water molecules around sodium can participate without adopting a distorted conformation.

This indicates that, for the NaCl solution, the examination of the HB pattern at k = 3 is not crucial. Therefore, the higher connectivity of water molecules around chloride for k = 2 and n₁ = 3 and 4 originates from the same mechanism identified for k = 1: the formation of a hydrogen bond at k = 2 requires a distortion of the water molecule around sodium, but not around chloride.

Conclusions

In this work, we present a graph theory approach to encode and quantify the connectivity in any hydrogen bond network. We demonstrate that this method accurately captures variations in the hydrogen bond network of water in the presence of solutes, specifically NaCl at different concentrations. Beyond assessing overall connectivity, the method provides detailed insights into the distinct behavior of water molecules around each type of ion. By separately analyzing the broadcaster and receiver roles of each molecule, it is also possible to infer the relative orientation of water molecules with respect to the ions, revealing structural organization within the solvation shells. This approach proves to be highly sensitive, capable of detecting even subtle changes in the hydrogen bond network, and offers a complete characterization of the hydrogen bond environment surrounding each water molecule. Its generality and efficiency make it well suited for studying any process involving hydrogen bond networks and their dynamic rearrangements, including solvation phenomena, ion-specific effects, and phase transitions. For future developments, extending the framework to explicitly include water–ion interactions as additional edges, or integrating complementary descriptors such as energy-based metrics,^11,12 will enable a more complete representation of the connectivity within the solution network.

Author contributions

Nico Di Fonte: data curation, formal analysis, methodology, software, writing – original draft, and writing – review & editing. Isabella Daidone: conceptualization, methodology, writing – original draft, and writing – review & editing. Laura Zanetti-Polzi: conceptualization, methodology, writing – original draft, and writing – review & editing. All authors contributed equally to all aspects of the work.

Conflicts of interest

The authors have no conflicts of interest to declare.

Data availability

The code used to compute the HBScore, as well as a code to compute 3, 4 and 5 membered cycles in the HB network, is publicly available at: https://github.com/NDiFonte/HBScore.

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6cp00965d.

Acknowledgements

The authors acknowledge the CINECA award under the ISCRA initiative, for the availability of high-performance computing resources and support (Project Code: IsCd6). This work was supported by the European Union – NextGenerationEU under the Italian Ministry of University and Research (MUR), PRIN2022-P2022MC742PNRR and CUP E53D23018440001.

References

1. A. Nilsson and L. G. Pettersson, Nat. Commun., 2015, 6, 8998.
2. P. Gallo, K. Amann-Winkel, C. A. Angell, M. A. Anisimov, F. Caupin, C. Chakravarty, E. Lascaris, T. Loerting, A. Z. Panagiotopoulos and J. Russo, et al., Chem. Rev., 2016, 116, 7463–7500.
3. F. Martelli, J. Chem. Phys., 2019, 150, 094506.
4. P. H. Poole, F. Sciortino, T. Grande, H. E. Stanley and C. A. Angell, Phys. Rev. Lett., 1994, 73, 1632.
5. D. Bandyopadhyay, S. Mohan, S. Ghosh and N. Choudhury, J. Phys. Chem. B, 2013, 117, 8831–8843.
6. E. Shiratani and M. Sasai, J. Chem. Phys., 1998, 108, 3264–3276.
7. J. R. Errington and P. G. Debenedetti, Nature, 2001, 409, 318–321.
8. M. J. Cuthbertson and P. H. Poole, Phys. Rev. Lett., 2011, 106, 115706.
9. J. Russo and H. Tanaka, Nat. Commun., 2014, 5, 3556.
10. J. M. Montes de Oca, F. Sciortino and G. A. Appignanesi, J. Chem. Phys., 2020, 152, 244503.
11. N. A. Loubet, A. R. Verde and G. A. Appignanesi, J. Chem. Phys., 2024, 160, 144502.
12. N. A. Loubet, A. R. Verde and G. A. Appignanesi, J. Chem. Phys., 2025, 162, 244703.
13. A. V. Muthachikavil, G. M. Kontogeorgis, X. Liang, Q. Lei and B. Peng, Phys. Rev. E, 2022, 105, 034604.
14. C. Faccio, M. Benzi, L. Zanetti-Polzi and I. Daidone, J. Mol. Liq., 2022, 355, 118922.
15. A. Neophytou, D. Chakrabarti and F. Sciortino, Nat. Phys., 2022, 18, 1248–1253.
16. R. Foffi, J. Russo and F. Sciortino, J. Chem. Phys., 2021, 154, 184506.
17. A. Rahman and F. H. Stillinger, J. Am. Chem. Soc., 1973, 95, 7943–7948.
18. M. Matsumoto, A. Baba and I. Ohmine, J. Chem. Phys., 2007, 127, 134504.
19. M. Formanek and F. Martelli, AIP Adv., 2020, 10, 055205.
20. M. Laliberte and W. E. Cooper, J. Chem. Eng. Data, 2004, 49, 1141–1151.
21. M. Laliberte, J. Chem. Eng. Data, 2007, 52, 321–335.
22. Y. Zhang and P. S. Cremer, Curr. Opin. Chem. Biol., 2006, 10, 658–663.
23. K. Tielrooij, N. Garcia-Araez, M. Bonn and H. Bakker, Science, 2010, 328, 1006–1009.
24. C. Zhang, S. Yue, A. Z. Panagiotopoulos, M. L. Klein and X. Wu, Nat. Commun., 2022, 13, 822.
25. N. Di Fonte, G. Dell'Orletta, L. Zanetti-Polzi and I. Daidone, Phys. Chem. Chem. Phys., 2025, 27, 18901–18906.
26. M. Flór, V. Vorobev, V. Mandalaparthy, N. F. van der Vegt, P. S. Cremer and S. Roke, J. Am. Chem. Soc., 2025, 147, 37328–37336.
27. J.-H. Choi, H. Lee, H. R. Choi and M. Cho, Annu. Rev. Phys. Chem., 2018, 69, 125–149.
28. Y. Ding, A. A. Hassanali and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 3310–3315.
29. I. Bakó, T. Megyes, S. Bálint, T. Grósz and V. Chihaia, Phys. Chem. Chem. Phys., 2008, 10, 5004–5011.
30. Y. Feng, H. Fang, Y. Gao and K. Ni, Phys. Chem. Chem. Phys., 2022, 24, 9707–9717.
31. A. Semmeq, K. Anand, A. Carof, A. Bastida and F. Ingrosso, Molecules, 2025, 30, 1678.
32. C. Faccio, N. Di Fonte, I. Daidone and L. Zanetti-Polzi, J. Mol. Liq., 2023, 392, 123425.
33. N. Di Fonte, C. Faccio, L. Zanetti-Polzi and I. Daidone, J. Chem. Phys., 2024, 161, 034505.
34. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1, 19–25.
35. S. Blazquez, M. Conde and C. Vega, J. Chem. Phys., 2023, 158, 054505.
36. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
37. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
38. T. Darden, D. York and L. Pedersen, et al., J. Chem. Phys., 1993, 98, 10089.
39. B. Hess, H. Bekker, H. J. Berendsen and J. G. Fraaije, J. Comp. Chem., 1997, 18, 1463–1472.
40. I.-C. Yeh and G. Hummer, J. Phys. Chem. B, 2004, 108, 15873–15879.
41. E. Estrada, The structure of complex networks: theory and applications, Oxford University Press, 2012.
42. A. Luzar and D. Chandler, Nature, 1996, 379, 55–57.
43. M. Benzi and C. Klymko, J. Complex Netw., 2013, 1, 124–149.
44. M. Benzi, I. Daidone, C. Faccio and L. Zanetti-Polzi, J. Complex Networks, 2023, 11, cnad001.
45. A. P. Gaiduk and G. Galli, J. Phys. Chem. Lett., 2017, 8, 1496–1502.

---