ChemNetworks: a complex network analysis tool for chemical systems

Networks of Hydrogen Bond Networks in Water Clusters

Water clusters play a prominent role in atmospheric and solution chemistry. Numerous arrangements of protons, H-bond configurations or networks, shape the cluster properties. Studies of small water clusters by cryogenic scanning tunneling microscopy and high-resolution rovibrational spectroscopy have established proton rearrangement mechanisms forming pathways between H-bond networks. The mechanisms, concerted tunneling in particular, describe the local processes connecting pairs of configurations. Here, proton rearrangement networks mapping these transformations are defined and explored to provide a global view of the H-bond configurations in clusters. The networks are constructed for clusters of different sizes and structures. Their analysis reveals an odd-even effect with respect to the number of water molecules, exponential growth of the small-world character, bimodality of the degree distributions, and gapped assortativity of the networks. The last two properties signify the unexpected division of H-bond configurations into two classes according to their network connectivity. The results demonstrate qualitative differences between proton rearrangement mechanisms, suggest a strong influence of the cluster structure. The generated networks are of interest as real-world models for network rewiring; they establish an alternative platform for studies of proton rearrangements in H-bonded systems.

Time-Resolved Graphs of Polymorphic Cycles for H-Bonded Network Identification in Flexible Biomolecules

A novel approach based on a coarse-grained representation of topological graphs is proposed for the automatic analysis of molecular dynamics (MD) trajectories of hydrogen-bonded (H-Bonded) flexible biomolecules. Herein, our approach models an H-Bonded biomolecule by its H-Bonded cycles and its graph of cycles in which the vertices and links represent the intersections between these cycles. We propose a methodology in which each identified conformer/isomer from the MD is represented by a well-chosen set of H-Bonded cycles called a minimum cycle basis. The key component is the "polycycles" that distinguish the cycles that play the same polymorphic role in the molecule from the ones that lead to an actual conformational change of the molecule. The relevance of our proposed method is evaluated on MD trajectories of gas-phase biomolecules, for which the covalent bonds are unchanged over time and only the hydrogen bonds change over time. The polygraphs and their time evolution are shown to reveal the dynamicity of the metastructure(s) of the H-Bonded biomolecules while providing polymorphic information on the cycles. Such information on the dynamics and changes in the H-bond network, as some cycles change identity while retaining the same role in the overall structure, is not easily captured at the atomic level of representation. Such information can instead be captured by polymorphic cycles.

Complex Networks and Interacting Particle Systems

Complex networks is a growing discipline aimed at understanding large interacting systems. One of its goals is to establish a relation between the interactions of a system and the networks structure that emerges. Taking a Lennard-Jones particle system as an example, we show that when interactions are governed by a potential, the notion of structure given by the physical arrangement of the interacting particles can be interpreted as a binary approximation to the interaction potential. This approximation simplifies the calculation of the partition function of the system and allows to study the stability of the interaction structure. We compare simulated results with those from the approximated partition function and show how the network and system perspective complement each other. With this, we draw a direct connection between the interactions of a molecular system and the network structure it forms and assess the degree to which it describes the system. We conclude by discussing the advantages and limitations of this method for weighted networks, as well as how this concept might be extended to more general systems.

Multicomponent Cholesky Decomposition: Application to Nuclear-Electronic Orbital Theory

The Cholesky decomposition technique is commonly used to reduce the memory requirement for storing two-particle repulsion integrals in quantum chemistry calculations that use atomic orbital bases. However, when quantum methods use multicomponent bases, such as nuclear-electronic orbitals, additional challenges are introduced due to asymmetric two-particle integrals. This work proposes several multicomponent Cholesky decomposition methods for calculations using nuclear-electronic orbital density functional theory. To analyze the errors in different Cholesky decomposition components, benchmark calculations using water clusters are carried out. The largest benchmark calculation is a water cluster (H2O)27 where all 54 protons are treated quantum mechanically. This study provides energetic and complexity analyses to demonstrate the accuracy and performance of the proposed multicomponent Cholesky decomposition method.

Algorithmic Graph Theory, Reinforcement Learning and Game Theory in MD Simulations: From 3D Structures to Topological 2D-Molecular Graphs (2D-MolGraphs) and Vice Versa

This paper reviews graph-theory-based methods that were recently developed in our group for post-processing molecular dynamics trajectories. We show that the use of algorithmic graph theory not only provides a direct and fast methodology to identify conformers sampled over time but also allows to follow the interconversions between the conformers through graphs of transitions in time. Examples of gas phase molecules and inhomogeneous aqueous solid interfaces are presented to demonstrate the power of topological 2D graphs and their versatility for post-processing molecular dynamics trajectories. An even more complex challenge is to predict 3D structures from topological 2D graphs. Our first attempts to tackle such a challenge are presented with the development of game theory and reinforcement learning methods for predicting the 3D structure of a gas-phase peptide.

ReNeGate: A Reaction Network Graph-Theoretical Tool for Automated Mechanistic Studies in Computational Homogeneous Catalysis

Exploration of the chemical reaction space of chemical transformations in multicomponent mixtures is one of the main challenges in contemporary computational chemistry. To remove expert bias from mechanistic studies and to discover new chemistries, an automated graph-theoretical methodology is proposed, which puts forward a network formalism of homogeneous catalysis reactions and utilizes a network analysis tool for mechanistic studies. The method can be used for analyzing trajectories with single and multiple catalytic species and can provide unique conformers of catalysts including multinuclear catalyst clusters along with other catalytic mixture components. The presented three-step approach has the integrated ability to handle multicomponent catalytic systems of arbitrary complexity (mixtures of reactants, catalyst precursors, ligands, additives, and solvents). It is not limited to predefined chemical rules, does not require prealignment of reaction mixture components consistent with a reaction coordinate, and is not agnostic to the chemical nature of transformations. Conformer exploration, reactive event identification, and reaction network analysis are the main steps taken for identifying the pathways in catalytic systems given the starting precatalytic reaction mixture as the input. Such a methodology allows us to efficiently explore catalytic systems in realistic conditions for either previously observed or completely unknown reactive events in the context of a network representing different intermediates. Our workflow for the catalytic reaction space exploration exclusively focuses on the identification of thermodynamically feasible conversion channels, representative of the (secondary) catalyst deactivation or inhibition paths, which are usually most difficult to anticipate based solely on expert chemical knowledge. Thus, the expert bias is sought to be removed at all steps, and the chemical intuition is limited to the choice of the thermodynamic constraint imposed by the applicable experimental conditions in terms of threshold energy values for allowed transformations. The capabilities of the proposed methodology have been tested by exploring the reactivity of Mn complexes relevant for catalytic hydrogenation chemistry to verify previously postulated activation mechanisms and unravel unexpected reaction channels relevant to rare deactivation events.

Design principles of PI(4,5)P2 clustering under protein-free conditions: Specific cation effects and calcium-potassium synergy

SignificanceClustering of phosphatidylinositol 4,5-bisphosphate (PIP₂) with proteins into what are known as "PIP₂ rafts" is a critical component of intracellular signaling, yet little is known about PIP₂ clusters at the atomic level. Using molecular dynamics simulations and network theory, this paper shows that Ca²⁺ generates large clusters by linking PIP₂ dimers already formed by doubly charged P4/P5 phosphates, while monovalent cations form smaller and less-stable clusters by adding PIP₂ monomers preferentially via weaker interactions with P4/P5 (for Na⁺) or with glycerol P1 (for K⁺). Synergy arises between K⁺ and Ca²⁺ because each ion forms linkages with different phosphates, thereby giving clusters more ways to grow. This explains why Ca²⁺ is pumped into cells by ion channels to form PIP₂ rafts.

Thermal conductivity temperature dependence of water confined in nanoporous silicon

Recently, it has been shown that high density nanoconfined water was the reason of the important enhancement of the effective thermal conductivity up to a factor of 50% of a nanoporous silicon filled with water. In this work, using molecular dynamics simulations, we further investigate the role of the temperatureT(from 285 to 360 K) on the thermal conductivity enhancement of nanohybrid porous silicon and water system. Furthermore, by studying and analysing several structural and dynamical parameters of the nanoconfined water, we give physical insights of the observed phenomena. Upon increasing the temperature of the system, the thermal conductivity of the hybrid system increases reaching a maximum forT= 300 K. With this article, we prove the existence of new heat flux channels between a solid matrix and a nanoconfined liquid, with clear signatures both in the radial distribution function, mean square displacements, water molecules orientation, hydrogen bond networks and phonon density of states.

Identifying nonadditive contributions to the hydrophobicity of chemically heterogeneous surfaces via dual-loop active learning

Hydrophobic interactions drive numerous biological and synthetic processes. The materials used in these processes often possess chemically heterogeneous surfaces that are characterized by diverse chemical groups positioned in close proximity at the nanoscale; examples include functionalized nanomaterials and biomolecules, such as proteins and peptides. Nonadditive contributions to the hydrophobicity of such surfaces depend on the chemical identities and spatial patterns of polar and nonpolar groups in ways that remain poorly understood. Here, we develop a dual-loop active learning framework that combines a fast reduced-accuracy method (a convolutional neural network) with a slow higher-accuracy method (molecular dynamics simulations with enhanced sampling) to efficiently predict the hydration free energy, a thermodynamic descriptor of hydrophobicity, for nearly 200 000 chemically heterogeneous self-assembled monolayers (SAMs). Analysis of this dataset reveals that SAMs with distinct polar groups exhibit substantial variations in hydrophobicity as a function of their composition and patterning, but the clustering of nonpolar groups is a common signature of highly hydrophobic patterns. Further molecular dynamics analysis relates such clustering to the perturbation of interfacial water structure. These results provide new insight into the influence of chemical heterogeneity on hydrophobicity via quantitative analysis of a large set of surfaces, enabled by the active learning approach.

Protonation of Serine in Gas and Condensed and Microsolvated States in Aqueous Solution

Identification of molecules and elucidation of their chemical structure are ubiquitous problems in chemistry. Mass spectrometry (MS) can be used due to its sensitivity and versatility. For detection to occur, analytes must be ionized and transferred to the gas phase. Soft ionization processes such as electrospray ionization are popular; however, resulting microsolvated phases can alter the chemistry of analytes and therefore detection and identification. To understand these processes, we use computational methods to probe the ionization propensity of serine in the gas phase, aqueous microsolvated clusters, and aqueous solution. We show that the tautomeric form of serine is altered by the presence of water, as five water molecules can stabilize the zwitterionic tautomer. Inclusion of cosolutes such as ions can stabilize the zwitterion with as few as one or two water molecules present. We demonstrate that ionization propensity, as measured by gas phase bacisity, can increase by over 100 kJ/mol when placed in a small water-serine cluster, showing the sensitivity of the chemistry of microsolvated analytes. Finally, detailed analysis reveals that small droplets (less than seven water molecules) are extremely sensitive to addition of further water molecules. Beyond this limit, structural and electronic properties change little with droplet size.

Liquid Structure of a Water-in-Salt Electrolyte with a Remarkably Asymmetric Anion

Water-in-salt systems, i.e., super-concentrated aqueous electrolytes, such as lithium bis(trifluoromethanesulfonyl)imide (21 mol/kg_water), have been recently discovered to exhibit unexpectedly large electrochemical windows and high lithium transference numbers, thus paving the way to safe and sustainable charge storage devices. The peculiar transport features in these electrolytes are influenced by their intrinsically nanoseparated morphology, stemming from the anion hydrophobic nature and manifesting as nanosegregation between anions and water domains. The underlying mechanism behind this structure–dynamics correlation is, however, still a matter of strong debate. Here, we enhance the apolar nature of the anions, exploring the properties of the aqueous electrolytes of lithium salts with a strongly asymmetric anion, namely, (trifluoromethylsulfonyl)(nonafluorobutylsulfonyl) imide. Using a synergy of experimental and computational tools, we detect a remarkable level of structural heterogeneity at a mesoscopic level between anion-rich and water-rich domains. Such a ubiquitous sponge-like, bicontinuous morphology develops across the whole concentration range, evolving from large fluorinated globules at high dilution to a percolating fluorous matrix intercalated by water nanowires at super-concentrated regimes. Even at extremely concentrated conditions, a large population of fully hydrated lithium ions, with no anion coordination, is detected. One can then derive that the concomitant coexistence of (i) a mesoscopically segregated structure and (ii) fully hydrated lithium clusters disentangled from anion coordination enables the peculiar lithium diffusion features that characterize water-in-salt systems.

Some opinions on MD-based vibrational spectroscopy of gas phase molecules and their assembly: An overview of what has been achieved and where to go

We hereby review molecular dynamics simulations for anharmonic gas phase spectroscopy and provide some of our opinions of where the field is heading. With these new directions, the theoretical IR/Raman spectroscopy of large (bio)-molecular systems will be more easily achievable over longer time-scale MD trajectories for an increase in accuracy of the MD-IR and MD-Raman calculated spectra. With the new directions presented here, the high throughput 'decoding' of experimental IR/Raman spectra into 3D-structures should thus be possible, hence advancing e.g. the field of MS-IR for structural characterization by spectroscopy. We also review the assignment of vibrational spectra in terms of anharmonic molecular modes from the MD trajectories, and especially introduce our recent developments based on Graph Theory algorithms. Graph Theory algorithmic is also introduced in this review for the identification of the molecular 3D-structures sampled over MD trajectories.

Unexpected inverse correlations and cooperativity in ion-pair phase transfer

Liquid/liquid extraction is one of the most widely used separation and purification methods, where a forefront of research is the study of transport mechanisms for solute partitioning and the relationships that these have to solution structure at the phase boundary. To date, organized surface features that include protrusions, water-fingers, and molecular hinges have been reported. Many of these equilibrium studies have focused upon small-molecule transport – yet the extent to which the complexity of the solute, and the competition between different solutes, influence transport mechanisms have not been explored. Here we report molecular dynamics simulations that demonstrate that a metal salt (LiNO₃) can be transported via a protrusion mechanism that is remarkably similar to that reported for H₂O by tri-butyl phosphate (TBP), a process that involves dimeric assemblies. Yet the LiNO₃ out-competes H₂O for a bridging position between the extracting TBP dimer, which in-turn changes the preferred transport pathway of H₂O. Examining the electrolyte concentration dependence on ion-pair transport unexpectedly reveals an inverse correlation with the extracting surfactant concentration. As [LiNO₃] increases, surface adsorbed TBP becomes a limiting reactant in correlation with an increased negative surface charge induced by excess interfacial NO₃⁻, however the rate of transport is enhanced. Within the highly dynamic interfacial environment, we hypothesize that this unique cooperative effect may be due to perturbed surface organization that either decreases the energy of formation of transporting protrusion motifs or makes it easier for these self-assembled species to disengage from the surface.

A forefront of research in separations science (specifically liquid–liquid extraction) is the study of transport mechanisms for solute partitioning, and the relationships that these have to solution structure at the phase boundary.

Pressure Effects on the Relaxation of an Excited Ethane Molecule in High-Pressure Bath Gases

We use molecular dynamics to calculate the rotational and vibrational energy relaxation of C₂H₆ in Ar, Kr, and Xe bath gases over a pressure range of 10-400 atm and at temperatures of 300 and 800 K. The C₂H₆ is instantaneously excited by 80 kcal/mol randomly distributed into both vibrational and rotational modes. The computed relaxation rates show little sensitivity to the identity of the noble gas in the bath. Vibrational relaxation rates show a nonlinear pressure dependence at 300 K. At 800 K the reduced range of bath gas densities covered by the range of pressures does not yet show any nonlinearity in the pressure dependence. Rotational relaxation is characterized with two relaxation rates. The slower rate is comparable to the vibrational relaxation rate. The faster rate has a linear pressure dependence at 300 K but an irregular, nonlinear pressure dependence at 800 K. To understand this, a model was developed based on approximating the periodic box used in the molecular dynamics simulations by an equal-volume collection of cubes where each cube is sized to allow only single occupancy by the noble gas or the molecule. Combinatorial statistics then leads to a pressure- and temperature-dependent analytic distribution of the bath gas species the molecule encounters in a collision. This distribution, the dissociation energy of molecule/bath gas complexes and bath gas clusters, and the computed energy release per collision combine to show that only at 300 K is the energy release sufficient to dissociate likely complexes and clusters. This suggests that persistent and pressure-dependent clusters and complexes at 800 K may be responsible for the nonlinear pressure dependence of rotational relaxation.

Determination of kinetic properties in unimolecular dissociation of complex systems from graph theory based analysis of an ensemble of reactive trajectories

In this paper, we report how graph theory can be used to analyze an ensemble of independent molecular trajectories, which can react during the simulation time-length, and obtain structural and kinetic information. This method is totally general and here is applied to the prototypical case of gas phase fragmentation of protonated cyclo-di-glycine. This methodology allows us to analyze the whole set of trajectories in an automatic computer-based way without the need of visual inspection but by getting all the needed information. In particular, we not only determine the appearance of different products and intermediates but also characterize the corresponding kinetics. The use of colored graph and canonical labeling allows for the correct characterization of the chemical species involved. In the present case, the simulations consist of an ensemble of unimolecular fragmentation trajectories at constant energy such that from the rate constants at different energies, the threshold energy can also be obtained for both global and specific pathways. This approach allows for the characterization of ion-molecule complexes, likely through a roaming mechanism, by properly taking into account the elusive nature of such species. Finally, it is possible to directly obtain the theoretical mass spectrum of the fragmenting species if the reacting system is an ion as in the specific example.

Shear stress dependence of force networks in 3D dense suspensions

The geometric organization and force networks of 3D dense suspensions that exhibit both shear thinning and thickening have been examined as a function of varying strength of interparticle attractive interactions using lubrication flow discrete element simulations. Significant rearrangement of the geometric topology does not occur at either the local or global scale as these systems transition across the shear thinning and shear thickening regimes. In contrast, massive rearrangements in the balance of attractive, lubrication, and contact forces are observed with interesting behavior of network growth and competition. In agreement with prior work, in shear thinning regions the attractive force is dominant, however as the shear thickening region is approached there is growth of lubrication forces. Lubrication forces oppose the attraction forces, but as viscosity continues to increase under increasing shear stress, the lubrication forces are dominated by contact forces that also resist attraction. Contact forces are the dominant interactions during shear thickening and are an order of magnitude higher than their values in the shear-thinning regime. At high attractive interaction strength, contact networks can form even under shear thinning conditions, however high shear stress is still required before contact networks become the driving mechanism of shear thickening. Analysis of the contact force network during shear thickening generally indicates a uniformly spreading network that rapidly forms across empty domains; however the growth patterns exhibit structure that is significantly dependent upon the strength of interparticle interactions, indicating subtle variations in the mechanism of shear thickening.

The Middle Science: Traversing Scale In Complex Many-Body Systems

A roadmap is developed that integrates simulation methodology and data science methods to target new theories that traverse the multiple length- and time-scale features of many-body phenomena.

Impact of Water Extraction on Malonamide Aggregation: A Molecular Dynamics and Graph Theoretic Approach

Solution structure in liquid-liquid extraction affects the efficacy of separation; however, even for simplified organic phases, structural characterization and attribution of aggregation to intermolecular interactions are fundamental challenges. We investigate water uptake into organic phases for two malonamides commonly applied to actinide and lanthanide separations. Extracted water induces reorganization of the amphiphilic extractant molecules, although we find this rearrangement is not strongly manifested in small-angle X-ray scattering making it challenging to probe without methods such as atomistic simulation. Using a graph theoretic approach to define hydrogen bonded water/malonamide aggregates from molecular dynamics simulations, we find evidence of a characteristic aggregate size by water number that results from geometric accommodation of the surrounding malonamide molecules. This implies a degree of size selectivity inherent to these water-in-oil aggregates. Conversely, we find no evidence of a characteristic size of the aggregates with respect to their malonamide number. By defining a separate graphical representation of self-association of the amphiphilic malonamides, we quantify how water affects the local and nonlocal topology of the malonamide network, providing a basis for characterization of the structure and impact of polar solutes in increasingly complex organic phases.

GGL-Tox: Geometric Graph Learning for Toxicity Prediction

Toxicity analysis is a major challenge in drug design and discovery. Recently significant progress has been made through machine learning due to its accuracy, efficiency, and lower cost. US Toxicology in the 21st Century (Tox21) screened a large library of compounds, including approximately 12 000 environmental chemicals and drugs, for different mechanisms responsible for eliciting toxic effects. The Tox21 Data Challenge offered a platform to evaluate different computational methods for toxicity predictions. Inspired by the success of multiscale weighted colored graph (MWCG) theory in protein-ligand binding affinity predictions, we consider MWCG theory for toxicity analysis. In the present work, we develop a geometric graph learning toxicity (GGL-Tox) model by integrating MWCG features and the gradient boosting decision tree (GBDT) algorithm. The benchmark tests of the Tox21 Data Challenge are employed to demonstrate the utility and usefulness of the proposed GGL-Tox model. An extensive comparison with other state-of-the-art models indicates that GGL-Tox is an accurate and efficient model for toxicity analysis and prediction.

Molecular-scale origins of solution nanostructure and excess thermodynamic properties in a water/amphiphile mixture

The molecular and nanoscale origins of nonideality in excess thermodynamic properties are essential to understanding cosolvent mixtures, yet they remain challenging to determine. Here, we consider a binary mixture of water and an amphiphile, N,N,N',N'-tetramethylmalonamide (TMMA), which is characterized by strong hydrogen bonding between the two components and no hydrogen bonding between amphiphiles. Using molecular dynamics simulation, validated with excess volume measurements and X-ray scattering, we identify three distinct solution regimes across the composition range of the binary mixture and find that the transition between two of these regimes, marked by the water percolation threshold, is closely correlated with minima in the excess volume and excess enthalpy. Structural analysis of the simulations reveals an interplay between local interactions and solution nanostructure, determined by the relative strength of the water-water and water-amphiphile hydrogen bonding interactions. By comparison with other amphiphiles, such as linear alcohols, the relative strength of like and unlike interactions between water and amphiphile affects the relationship between thermodynamics and structural regimes. This provides insight into how molecular forces of mutual solvation interact across length scales and how they manifest in excess thermodynamic properties.

Thermodynamics of Malonamide Aggregation Deduced from Molecular Dynamics Simulations

The aggregation of malonamide extractants diluted in an aliphatic solvent phase has been studied in the presence of water by molecular dynamics simulation. Using association criteria based on distances between molecules and graphs theory, the aggregate distribution has been computed and the corresponding Gibbs energy of aggregates and mass action law constants have been determined. Finally, a model allowing us to the compute critical micelle concentration and osmotic data for a variable concentration of extractants, with or without a correction of the organic phase activity, was developed. It appears however that the accurate depiction of the aggregation allows modeling the thermodynamics of the solution even without an explicit calculation of the activity: both models give results in good agreement with the experiments.

An octanol hinge opens the door to water transport

Despite their prevalent use as a surrogate for partitioning of pharmacologically active solutes across lipid membranes, the mechanism of transport across water/octanol phase boundaries has remained unexplored. Using molecular dynamics, graph theoretical, cluster analysis, and Langevin dynamics, we reveal an elegant mechanism for the simplest solute, water. Self-assembled octanol at the interface reversibly binds water and swings like the hinge of a door to bring water into a semi-organized second interfacial layer (a “bilayer island”). This mechanism is distinct from well-known lipid flipping and water transport processes in protein-free membranes, highlighting important limitations in the water/octanol proxy. Interestingly, the collective and reversible behavior is well-described by a double well potential energy function, with the two stable states being the water bound to the hinge on either side of the interface. The function of the hinge for transport, coupled with the underlying double well energy landscape, is akin to a molecular switch or shuttle that functions under equilibrium and is driven by the differential free energies of solvation of H₂O across the interface. This example successfully operates within the dynamic motion of instantaneous surface fluctuations, a feature that expands upon traditional approaches toward controlled solute transport that act to avoid or circumvent the dynamic nature of the interface.

Despite their pharmacological relevance, the mechanism of transport across water/octanol phase boundaries has remained unexplored. Octanol molecular assemblies are demonstrated to reversibly bind water and swing like the hinge of a door.

The "Hole" Story in Ionized Water from the Perspective of Ehrenfest Dynamics

The radiolysis of liquid water and the radiation-matter interactions that happen in aqueous environments are important to the fields of chemistry, materials, and environmental sciences, as well as the biological and physiological response to extreme conditions and medical treatments. The initial stage of radiolysis is the ultrafast response, or hole dynamics, that triggers chemical processes within complex energetic landscapes that may include reactivity. A fundamental understanding necessitates the use of theoretical methods that are capable of simulating both ultrafast coherence and non-adiabatic energy transfer pathways. In this work, we carry out an ab initio Ehrenfest dynamics study to provide a more complete description of the ultrafast dynamics and reactive events initiated by photoionization of water. After sudden ionization, a range of processes, including hole trapping and transfer, large OH oscillations, proton transfer and subsequent relay, formation of the metastable Zundel complex, and long-lived coherence, are identified and new insight into their driving forces is elucidated.

A review of mathematical representations of biomolecular data

Recently, machine learning (ML) has established itself in various worldwide benchmarking competitions in computational biology, including Critical Assessment of Structure Prediction (CASP) and Drug Design Data Resource (D3R) Grand Challenges. However, the intricate structural complexity and high ML dimensionality of biomolecular datasets obstruct the efficient application of ML algorithms in the field. In addition to data and algorithm, an efficient ML machinery for biomolecular predictions must include structural representation as an indispensable component. Mathematical representations that simplify the biomolecular structural complexity and reduce ML dimensionality have emerged as a prime winner in D3R Grand Challenges. This review is devoted to the recent advances in developing low-dimensional and scalable mathematical representations of biomolecules in our laboratory. We discuss three classes of mathematical approaches, including algebraic topology, differential geometry, and graph theory. We elucidate how the physical and biological challenges have guided the evolution and development of these mathematical apparatuses for massive and diverse biomolecular data. We focus the performance analysis on protein-ligand binding predictions in this review although these methods have had tremendous success in many other applications, such as protein classification, virtual screening, and the predictions of solubility, solvation free energies, toxicity, partition coefficients, protein folding stability changes upon mutation, etc.

Interfacial structure in the liquid-liquid extraction of rare earth elements by phosphoric acid ligands: a molecular dynamics study

Solvent extraction (SX), wherein two immiscible liquids, one containing the extractant molecules and the other containing the solute to be extracted are brought in contact to effect the phase transfer of the solute, underpins metal extraction and recovery processes. The interfacial region is of utmost importance in the SX process, since besides thermodynamics, the physical and chemical heterogeneity at the interface governs the kinetics of the process. Yet, a fundamental understanding of this heterogeneity and its implications for the extraction mechanism are currently lacking. We use molecular dynamics (MD) simulations to study the liquid-liquid interface under conditions relevant to the SX of Rare Earth Elements (REEs) by a phosphoric acid ligand. Simulations revealed that the extractant molecules and varying amounts of acid and metal ions partitioned to the interface. The presence of these species had a significant effect on the interfacial thickness, hydrogen bond life times and orientations of the water molecules at the interface. Deprotonation of the ligands was essential for the adsorption of the metal ions at the interface, with these ions forming a number of different complexes at the interface involving one to three extractant molecules and four to eight water molecules. Although the interface itself was rough, no obvious 'finger-like' water protrusions penetrating the organic phase were seen in our simulations. While the results of our work help us gain fundamental insights into the sequence of events leading to the formation of a variety of interfacial complexes, they also emphasize the need to carry out a more detailed atomic level study to understand the full mechanism of extraction of REEs from the aqueous to organic phases by phosphoric acid ligands.

Characterization of Specific Ion Effects on PI(4,5)P2 Clustering: Molecular Dynamics Simulations and Graph-Theoretic Analysis

Numerous cellular functions mediated by phosphatidylinositol (4,5)-bisphosphate (PI(4,5)P2; PIP2) involve clustering of the lipid as well as colocalization with other lipids. Although the cation-mediated electrostatic interaction is regarded as the primary clustering mechanism, the ion-specific nature of the intermolecular network formation makes it challenging to characterize the clusters. Here we use all-atom molecular dynamics (MD) simulations of PIP2 monolayers and graph-theoretic analysis to gain insight into the phenomenon. MD simulations reveal that the intermolecular interactions preferentially occur between specific cations and phosphate groups (P1, P4, and P5) of the inositol headgroup with better-matched kosmotropic/chaotropic characters consistent with the law of matching water affinities (LMWA). Ca2+ is strongly attracted to P4/P5, while K+ preferentially binds to P1; Na+ interacts with both P4/P5 and P1. These specific interactions lead to the characteristic clustering patterns. Specificially, the size distributions and structures of PIP2 clusters generated by kosmotropic cations Ca2+ and Na+ are bimodal, with a combination of small and large clusters, while there is little clustering in the presence of only chaotropic K+; the largest clusters are obtained in systems with all three cations. The small-world network (a model with both local and long-range connections) best characterizes the clusters, followed by the random and the scale-free networks. More generally, the present results interpreted within the LMWA are consistent with the relative eukaryotic intracellular concentrations Ca2+ ≪ Na+ < Mg2+ < K+; that is, concentrations of Ca2+ and Na+ must be low to prevent damaging aggregation of lipids, DNA, RNA and phosphate-containing proteins.

MathDL: mathematical deep learning for D3R Grand Challenge 4

We present the performances of our mathematical deep learning (MathDL) models for D3R Grand Challenge 4 (GC4). This challenge involves pose prediction, affinity ranking, and free energy estimation for beta secretase 1 (BACE) as well as affinity ranking and free energy estimation for Cathepsin S (CatS). We have developed advanced mathematics, namely differential geometry, algebraic graph, and/or algebraic topology, to accurately and efficiently encode high dimensional physical/chemical interactions into scalable low-dimensional rotational and translational invariant representations. These representations are integrated with deep learning models, such as generative adversarial networks (GAN) and convolutional neural networks (CNN) for pose prediction and energy evaluation, respectively. Overall, our MathDL models achieved the top place in pose prediction for BACE ligands in Stage 1a. Moreover, our submissions obtained the highest Spearman correlation coefficient on the affinity ranking of 460 CatS compounds, and the smallest centered root mean square error on the free energy set of 39 CatS molecules. It is worthy to mention that our method on docking pose predictions has significantly improved from our previous ones.

Properties of the tetravalent actinide series in aqueous phase from a microscopic simulation self-consistent engine

In the context of nuclear fuel recycling and environmental issues, the understanding of the properties of radio-elements with various approaches remains a challenge regarding their dangerousness. Moreover, experimentally, some issues are also of importance; first, it is imperative to work at sufficiently high concentrations to reach the sensitivities of the analytical tools, however this condition often leads to precipitation for some of them; second, stabilizing specific oxidation states of some actinides remains a challenge, thus making it difficult to extract general trends across the actinide series. Complementary to experiments, modeling can be used to unbiasedly probe the actinide's properties in an aquatic environment and offers a predictive tool. We report the first molecular dynamics simulations based on homogeneously built force fields for the whole series of the tetravalent actinides in aqueous phase from ThIV to BkIV and including PuIV. The force fields used to model the interactions among the constituents include polarization and charge donation microscopic effects. They are built from a self-consistent iterative ab initio based engine that can be included in future developments as an element of a potential machine learning procedure devoted to generating accurate force fields. The comparison of our simulated hydrated actinide properties to available experimental data shows the model robustness and the relevance of our parameter assignment engine. Moreover, our simulated structural, dynamical and evolution of the hydration free energy data show that, apart from AmIV and CmIV, the actinide properties change progressively along the series.

AGL-Score: Algebraic Graph Learning Score for Protein-Ligand Binding Scoring, Ranking, Docking, and Screening

Although algebraic graph theory-based models have been widely applied in physical modeling and molecular studies, they are typically incompetent in the analysis and prediction of biomolecular properties, confirming the common belief that "one cannot hear the shape of a drum". A new development in the century-old question about the spectrum-geometry relationship is provided. Novel algebraic graph learning score (AGL-Score) models are proposed to encode high-dimensional physical and biological information into intrinsically low-dimensional representations. The proposed AGL-Score models employ multiscale weighted colored subgraphs to describe crucial molecular and biomolecular interactions in terms of graph invariants derived from graph Laplacian, its pseudo-inverse, and adjacency matrices. Additionally, AGL-Score models are integrated with an advanced machine learning algorithm to predict biomolecular macroscopic properties from the low-dimensional graph representation of biomolecular structures. The proposed AGL-Score models are extensively validated for their scoring power, ranking power, docking power, and screening power via a number of benchmark datasets, namely CASF-2007, CASF-2013, and CASF-2016. Numerical results indicate that the proposed AGL-Score models are able to outperform other state-of-the-art scoring functions in protein-ligand binding scoring, ranking, docking, and screening. This study indicates that machine learning methods are powerful tools for molecular docking and virtual screening. It also indicates that spectral geometry or spectral graph theory has the ability to infer geometric properties.

The Role of Excited-State Proton Relays in the Photochemical Dynamics of Water Nanodroplets

In this work, we applied nonadiabatic excited-state molecular dynamics in tandem with ab initio electronic structure theory to illustrate a complete mechanistic landscape underpinning the ultraviolet absorption-initiated photochemical dynamics in water nanodroplets. The goal is to understand the nonequilibrium excited-state molecular dynamics initiated by the relaxation of a solvated photoelectron and consequential photochemical processes. The lowest-lying excited state shows the proton dissociation for a single water molecule forming intermediate hydronium complexes through a proton relay. At approximately 100 fs, the proton relay process gives rise to the relaxation of the excited state accompanied by a rapid increase in the nonadiabatic coupling strength with the ground state, and the nanodroplet nonradiatively decays. The nonadiabatic transition to the ground state produces excited vibrational states that facilitate the recombination of the dissociated proton and hydroxyl group, eventually leading to the desorption of water molecules from the nanodroplet. Additionally, lifetimes of transient photochemical events are also resolved for the relaxation of a solvated electron, excited-state proton relay, and nonradiative transition.

Chokepoints in Mechanical Coupling Associated with Allosteric Proteins: The Pyruvate Kinase Example

Although the critical role of allostery in controlling enzymatic processes is well appreciated, there is a current dearth in our understanding of its underlying mechanisms, including communication between binding sites. One potential key aspect of intersite communication is the mechanical coupling between residues in a protein. Here, we introduce a graph-based computational approach to investigate the mechanical coupling between distant parts of a protein, highlighting effective pathways via which protein motion can transfer energy between sites. In this method, each residue is treated as a node on a weighted, undirected graph, in which the edges are defined by locally correlated motions of those residues and weighted by the strength of the correlation. The method was validated against experimental data on allosteric regulation in the human liver pyruvate kinase as obtained from full-protein alanine-scanning mutagenesis (systematic mutation) studies, as well as computational data on two G-protein-coupled receptors. The method provides semiquantitative information on the regulatory importance of specific structural elements. It is shown that these elements are key for the mechanical coupling between distant parts of the protein by providing effective pathways for energy transfer. It is also shown that, although there are a multitude of energy transfer pathways between distant parts of a protein, these pathways share a few common nodes that represent effective "chokepoints" for the communication.

Global topology of contact force networks: Insight into shear thickening suspensions

Highly concentrated particle suspensions (also called slurries) can undergo a sharp increase in viscosity, or shear thickening, under applied stress. Understanding the fundamental features leading to such rheological change is crucial to optimize flow conditions or to design flow modifiers for slurry processing. While local changes to the particle environment under applied shear can be related to changes in viscosity, there is a broader need to connect the shear thickening transition to the fundamental organization of particle-interaction forces which lead to long-range organization. In particular, at a high volume fraction of particles, recent evidence indicates frictional forces between contacting particles is of importance. Herein, the network of frictional contact forces is analyzed within simulated two-dimensional shear thickening suspensions. Two topological metrics are studied to characterize the response of the contact force network (CFN) under varying applied shear stress. The metrics, geodesic index and the void parameter, reflect complementary aspects of the CFN: One is the connectedness of the contact network and the second is the distribution of spatial areas devoid of particle-particle contacts. Considered in relation to the variation of the viscosity, the topological metrics show that the network grows homogeneously at large scales but with many local regions devoid of contacts, indicating clearly the role of CFN growth in causing the large change in the rheological response at the shear thickening transition.

Surfactant-enhanced heterogeneity of the aqueous interface drives water extraction into organic solvents

Liquid/liquid extraction (LLE) is one of the most industrially relevant separations methods. Adsorbed surfactant is demonstrated to enhance interfacial heterogeneity and lead to water protrusions that form the basis for transport into the organic phase.

Metanetwork Transmission Model for Predicting a Malaria-Control Strategy

Background: Mosquitoes are the primary vectors responsible for malaria transmission to humans, with numerous experiments having been conducted to aid in the control of malaria transmission. One of the main approaches aims to develop malaria parasite resistance within the mosquito population by introducing a resistance (R) allele. However, when considering this approach, some critical factors, such as the life of the mosquito, female mosquito fertility capacity, and human and mosquito mobility, have not been considered. Thus, an understanding of how mosquitoes and humans affect disease dynamics is needed to better inform malaria control policymaking. Methods: In this study, a method was proposed to create a metanetwork on the basis of the geographic maps of Gambia, and a model was constructed to simulate evolution within a mixed population, with factors such as birth, death, reproduction, biting, infection, incubation, recovery, and transmission between populations considered in the network metrics. First, the same number of refractory mosquitoes (RR genotype) was introduced into each population, and the prevalence of the R allele (the ratio of resistant alleles to all alleles) and malaria were examined. In addition, a series of simulations were performed to evaluate two different deployment strategies for the reduction of the prevalence of malaria. The R allele and malaria prevalence were calculated for both the strategies, with 10,000 refractory mosquitoes deployed into randomly selected populations or selection based on nodes with top-betweenness values. The 10,000 mosquitoes were deployed among 1, 5, 10, 20, or 40 populations. Results: The simulations in this paper showed that a higher RR genotype (resistant-resistant genes) ratio leads to a higher R allele prevalence and lowers malaria prevalence. Considering the cost of deployment, the simulation was performed with 10,000 refractory mosquitoes deployed among 1 or 5 populations, but this approach did not reduce the original malaria prevalence. Thus, instead, the 10,000 refractory mosquitoes were distributed among 10, 20, or 40 populations and were shown to effectively reduce the original malaria prevalence. Thus, deployment among a relatively small fraction of central nodes can offer an effective strategy to reduce malaria. Conclusion: The standard network centrality measure is suitable for planning the deployment of refractory mosquitoes. Importance: Malaria is an infectious disease that is caused by a plasmodial parasite, and some control strategies have focused on genetically modifying the mosquitoes. This work aims to create a model that takes into account mosquito development and malaria transmission among the population and how these factors influence disease dynamics so as to better inform malaria-control policymaking.

Competitive Interactions Within Cm(III) Solvation in Binary Water/Methanol Solutions

Competitive forces exist in multicomponent solutions, and within electrolytes they consist of both ion-solvent and solvent-solvent interactions. These can influence a myriad of processes, including ligand complexation. In the case of water/alcohol solutions, recent work revealed an interesting dilemma regarding the overall solution dynamics and organization as compared to solute-solvent interactions. This is particularly true for highly charged ions in solution, whose ion-solvent interactions were demonstrated to be highly sensitive to the composition of the immediate solvation environment. Faster solvent exchange should be observed about the ion, considering that second-order Møller-Plesset perturbation theory predicts an average decrease in ion-solvent dissociation energy when methanol enters the first solvation shell of Cm³⁺_(aq). Yet the addition of methanol to water causes the dynamic features of the hydrogen-bond network of the entire solution to slow. The apparent competition between these contrary forces was examined using a combination of electronic structure calculations with both ab initio and classical molecular dynamics simulations, using binary water/methanol solutions and Cm³⁺ as a representative solute. This combination of theoretical methods predicts that, among the competitive effects of the solvent-solvent and ion-solvent interactions, the solution-phase dynamics imparted by the addition of methanol to water kinetically restricts the solvation exchange rates about Cm³⁺ in these binary solutions.

Molecular Dynamics on Wood-Derived Lignans Analyzed by Intermolecular Network Theory

The dynamics of interactions to a solvent is a key factor in the proper characterization of new molecular structures. In molecular dynamics simulations, the solvent molecules are explicitly present, thereby defining a more accurate description on how the solvent molecules affect the molecular conformation. Intermolecular interactions in chemical systems, e.g., hydrogen bonds, can be considered as networks or graphs. Graph theoretical analyses can be an outstanding tool in analyzing the changes in interactions between solvent and solute. In this study, the software ChemNetworks is applied to interaction studies between TIP4P solvent molecules and organic solutes, i.e., wood-derived lignan-based ligands called LIGNOLs, thereby supporting the research of interaction networks between organic molecules and solvents. This new approach is established by careful comparisons to studies using previously available tools. In the hydration studies, tetramethyl 1,4-diol is found to be the LIGNOL which was most likely to form hydrogen bonds to the TIP4P solvent.

Graph-Theoretic Analysis of Monomethyl Phosphate Clustering in Ionic Solutions

All-atom molecular dynamics simulations combined with graph-theoretic analysis reveal that clustering of monomethyl phosphate dianion (MMP^2-) is strongly influenced by the types and combinations of cations in the aqueous solution. Although Ca²⁺ promotes the formation of stable and large MMP^2- clusters, K⁺ alone does not. Nonetheless, clusters are larger and their link lifetimes are longer in mixtures of K⁺ and Ca²⁺. This "synergistic" effect depends sensitively on the Lennard-Jones interaction parameters between Ca²⁺ and the phosphorus oxygen and correlates with the hydration of the clusters. The pronounced MMP^2- clustering effect of Ca²⁺ in the presence of K⁺ is confirmed by Fourier transform infrared spectroscopy. The characterization of the cation-dependent clustering of MMP^2- provides a starting point for understanding cation-dependent clustering of phosphoinositides in cell membranes.

Water phase transitions from the perspective of hydrogen-bond network analysis

Analysis of the water phase transitions from the perspective of hydrogen bond networks.

Structural and Dynamic Heterogeneity of Capillary Wave Fronts at Aqueous Interfaces

Using a unique combination of slab-layering analyses and identification of truly interfacial molecules, this work examines water/vapor and water/n-hexane interfaces, specifically the structural and dynamic perturbations of the interfacial water molecules at different locations within the surface capillary waves. From both the structural and dynamic properties analyzed, it is found that these interfacial water molecules dominate the perturbations within the interfacial region, which can extend deep into the water phase relative to the Gibbs dividing surface. Of more importance is the demonstration of structural and dynamic heterogeneity of the interfacial water molecules at the capillary wave front, as indicated by the dipole orientation and the structural and dynamic behavior of hydrogen bonds and their networks.

Network analysis and percolation transition in hydrogen bonded clusters: nitric acid and water extracted by tributyl phosphate

Network analysis of hydrogen bonded clusters formed in simulation by extraction of nitric acid and water by TBP interprets cluster topologies and identifies the mechanism for third phase formation.