Learning reduced-order models for cardiovascular simulations with graph neural networks.

Reduced-order models based on physics are a popular choice in cardiovascular modeling due to their efficiency, but they may experience loss in accuracy when working with anatomies that contain numerous junctions or pathological conditions. We develop one-dimensional reduced-order models that simulate blood flow dynamics using a graph neural network trained on three-dimensional hemodynamic simulation data. Given the initial condition of the system, the network iteratively predicts the pressure and flow rate at the vessel centerline nodes. Our numerical results demonstrate the accuracy and generalizability of our method in physiological geometries comprising a variety of anatomies and boundary conditions. Our findings demonstrate that our approach can achieve errors below 3% for pressure and flow rate, provided there is adequate training data. As a result, our method exhibits superior performance compared to physics-based one-dimensional models while maintaining high efficiency at inference time.

Memory-Augmented Generative Adversarial Networks for Anomaly Detection.

We propose a memory-augmented deep learning model for semisupervised anomaly detection (AD). While many traditional AD methods focus on modeling the distribution of normal data, additional constraints in the modeling process are needed to distinguish between normal and abnormal data. The proposed model, named memory augmented generative adversarial networks (MEMGAN), is coupled with external memory units through attentional operations. One property of MEMGAN in the latent space is such that encoded normal data are expected to reside in the convex hull of the memory units, while the abnormal ones are separated outside. This property makes the AD process of MEMGAN more robust and reliable. Experiments on AD datasets adapted from MVTec, MNIST, CIFAR10, and Arrhythmia demonstrate that MEMGAN notably improves over previous AD models. We also find that the decoded memory units in MEMGAN are more diverse and interpretable than those in previous memory-augmented models.

The multi-dimensional generalized Langevin equation for conformational motion of proteins.

Using the generalized Langevin equation (GLE) is a promising approach to build coarse-grained (CG) models of molecular systems since the GLE model often leads to more accurate thermodynamic and kinetic predictions than Brownian dynamics or Langevin models by including a more sophisticated friction with memory. The GLE approach has been used for CG coordinates such as the center of mass of a group of atoms with pairwise decomposition and for a single CG coordinate. We present a GLE approach when CG coordinates are multiple generalized coordinates, defined, in general, as nonlinear functions of microscopic atomic coordinates. The CG model for multiple generalized coordinates is described by the multidimensional GLE from the Mori-Zwanzig formalism, which includes an exact memory matrix. We first present a method to compute the memory matrix in a multidimensional GLE using trajectories of a full system. Then, in order to reduce the computational cost of computing the multidimensional friction with memory, we introduce a method that maps the GLE to an extended Markovian system. In addition, we study the effect of using a nonconstant mass matrix in the CG model. In particular, we include mass-dependent terms in the mean force. We used the proposed CG model to describe the conformational motion of a solvated alanine dipeptide system, with two dihedral angles as the CG coordinates. We showed that the CG model can accurately reproduce two important kinetic quantities: the velocity autocorrelation function and the distribution of first passage times.

Investigating the role of non-covalent interactions in conformation and assembly of triazine-based sequence-defined polymers.

Grate and co-workers at Pacific Northwest National Laboratory recently developed high information content triazine-based sequence-defined polymers that are robust by not having hydrolyzable bonds and can encode structure and functionality by having various side chains. Through molecular dynamics (MD) simulations, the triazine polymers have been shown to form particular sequential stacks, have stable backbone-backbone interactions through hydrogen bonding and π-π interactions, and conserve their cis/trans conformations throughout the simulation. However, we do not know the effects of having different side chains and backbone structures on the entire conformation and whether the cis or trans conformation is more stable for the triazine polymers. For this reason, we investigate the role of non-covalent interactions for different side chains and backbone structures on the conformation and assembly of triazine polymers in MD simulations. Since there is a high energy barrier associated with the cis-trans isomerization, we use replica exchange molecular dynamics (REMD) to sample various conformations of triazine hexamers. To obtain rates and intermediate conformations, we use the recently developed concurrent adaptive sampling (CAS) algorithm for dimers of triazine trimers. We found that the hydrogen bonding ability of the backbone structure is critical for the triazine polymers to self-assemble into nanorod-like structures, rather than that of the side chains, which can help researchers design more robust materials.

Efficiently sampling conformations and pathways using the concurrent adaptive sampling (CAS) algorithm.

Molecular dynamics simulations are useful in obtaining thermodynamic and kinetic properties of bio-molecules, but they are limited by the time scale barrier. That is, we may not obtain properties' efficiently because we need to run microseconds or longer simulations using femtosecond time steps. To overcome this time scale barrier, we can use the weighted ensemble (WE) method, a powerful enhanced sampling method that efficiently samples thermodynamic and kinetic properties. However, the WE method requires an appropriate partitioning of phase space into discrete macrostates, which can be problematic when we have a high-dimensional collective space or when little is known a priori about the molecular system. Hence, we developed a new WE-based method, called the "Concurrent Adaptive Sampling (CAS) algorithm," to tackle these issues. The CAS algorithm is not constrained to use only one or two collective variables, unlike most reaction coordinate-dependent methods. Instead, it can use a large number of collective variables and adaptive macrostates to enhance the sampling in the high-dimensional space. This is especially useful for systems in which we do not know what the right reaction coordinates are, in which case we can use many collective variables to sample conformations and pathways. In addition, a clustering technique based on the committor function is used to accelerate sampling the slowest process in the molecular system. In this paper, we introduce the new method and show results from two-dimensional models and bio-molecules, specifically penta-alanine and a triazine trimer.

Computing the non-Markovian coarse-grained interactions derived from the Mori-Zwanzig formalism in molecular systems: Application to polymer melts.

Memory effects are often introduced during coarse-graining of a complex dynamical system. In particular, a generalized Langevin equation (GLE) for the coarse-grained (CG) system arises in the context of Mori-Zwanzig formalism. Upon a pairwise decomposition, GLE can be reformulated into its pairwise version, i.e., non-Markovian dissipative particle dynamics (DPD). GLE models the dynamics of a single coarse particle, while DPD considers the dynamics of many interacting CG particles, with both CG systems governed by non-Markovian interactions. We compare two different methods for the practical implementation of the non-Markovian interactions in GLE and DPD systems. More specifically, a direct evaluation of the non-Markovian (NM) terms is performed in LE-NM and DPD-NM models, which requires the storage of historical information that significantly increases computational complexity. Alternatively, we use a few auxiliary variables in LE-AUX and DPD-AUX models to replace the non-Markovian dynamics with a Markovian dynamics in a higher dimensional space, leading to a much reduced memory footprint and computational cost. In our numerical benchmarks, the GLE and non-Markovian DPD models are constructed from molecular dynamics (MD) simulations of star-polymer melts. Results show that a Markovian dynamics with auxiliary variables successfully generates equivalent non-Markovian dynamics consistent with the reference MD system, while maintaining a tractable computational cost. Also, transient subdiffusion of the star-polymers observed in the MD system can be reproduced by the coarse-grained models. The non-interacting particle models, LE-NM/AUX, are computationally much cheaper than the interacting particle models, DPD-NM/AUX. However, the pairwise models with momentum conservation are more appropriate for correctly reproducing the long-time hydrodynamics characterised by an algebraic decay in the velocity autocorrelation function.

A comparison of weighted ensemble and Markov state model methodologies.

Computation of reaction rates and elucidation of reaction mechanisms are two of the main goals of molecular dynamics (MD) and related simulation methods. Since it is time consuming to study reaction mechanisms over long time scales using brute force MD simulations, two ensemble methods, Markov State Models (MSMs) and Weighted Ensemble (WE), have been proposed to accelerate the procedure. Both approaches require clustering of microscopic configurations into networks of "macro-states" for different purposes. MSMs model a discretization of the original dynamics on the macro-states. Accuracy of the model significantly relies on the boundaries of macro-states. On the other hand, WE uses macro-states to formulate a resampling procedure that kills and splits MD simulations for achieving better efficiency of sampling. Comparing to MSMs, accuracy of WE rate predictions is less sensitive to the definition of macro-states. Rigorous numerical experiments using alanine dipeptide and penta-alanine support our analyses. It is shown that MSMs introduce significant biases in the computation of reaction rates, which depend on the boundaries of macro-states, and Accelerated Weighted Ensemble (AWE), a formulation of weighted ensemble that uses the notion of colors to compute fluxes, has reliable flux estimation on varying definitions of macro-states. Our results suggest that whereas MSMs provide a good idea of the metastable sets and visualization of overall dynamics, AWE provides reliable rate estimations requiring less efforts on defining macro-states on the high dimensional conformational space.

AWE-WQ: fast-forwarding molecular dynamics using the accelerated weighted ensemble.

A limitation of traditional molecular dynamics (MD) is that reaction rates are difficult to compute. This is due to the rarity of observing transitions between metastable states since high energy barriers trap the system in these states. Recently the weighted ensemble (WE) family of methods have emerged which can flexibly and efficiently sample conformational space without being trapped and allow calculation of unbiased rates. However, while WE can sample correctly and efficiently, a scalable implementation applicable to interesting biomolecular systems is not available. We provide here a GPLv2 implementation called AWE-WQ of a WE algorithm using the master/worker distributed computing WorkQueue (WQ) framework. AWE-WQ is scalable to thousands of nodes and supports dynamic allocation of computer resources, heterogeneous resource usage (such as central processing units (CPU) and graphical processing units (GPUs) concurrently), seamless heterogeneous cluster usage (i.e., campus grids and cloud providers), and support for arbitrary MD codes such as GROMACS, while ensuring that all statistics are unbiased. We applied AWE-WQ to a 34 residue protein which simulated 1.5 ms over 8 months with peak aggregate performance of 1000 ns/h. Comparison was done with a 200 µs simulation collected on a GPU over a similar timespan. The folding and unfolded rates were of comparable accuracy.

Folding Proteins at 500 ns/hour with Work Queue.

Molecular modeling is a field that traditionally has large computational costs. Until recently, most simulation techniques relied on long trajectories, which inherently have poor scalability. A new class of methods is proposed that requires only a large number of short calculations, and for which minimal communication between computer nodes is required. We considered one of the more accurate variants called Accelerated Weighted Ensemble Dynamics (AWE) and for which distributed computing can be made efficient. We implemented AWE using the Work Queue framework for task management and applied it to an all atom protein model (Fip35 WW domain). We can run with excellent scalability by simultaneously utilizing heterogeneous resources from multiple computing platforms such as clouds (Amazon EC2, Microsoft Azure), dedicated clusters, grids, on multiple architectures (CPU/GPU, 32/64bit), and in a dynamic environment in which processes are regularly added or removed from the pool. This has allowed us to achieve an aggregate sampling rate of over 500 ns/hour. As a comparison, a single process typically achieves 0.1 ns/hour.

Computing generalized Langevin equations and generalized Fokker-Planck equations.

The Mori-Zwanzig formalism is an effective tool to derive differential equations describing the evolution of a small number of resolved variables. In this paper we present its application to the derivation of generalized Langevin equations and generalized non-Markovian Fokker-Planck equations. We show how long time scales rates and metastable basins can be extracted from these equations. Numerical algorithms are proposed to discretize these equations. An important aspect is the numerical solution of the orthogonal dynamics equation which is a partial differential equation in a high dimensional space. We propose efficient numerical methods to solve this orthogonal dynamics equation. In addition, we present a projection formalism of the Mori-Zwanzig type that is applicable to discrete maps. Numerical applications are presented from the field of Hamiltonian systems.

High-ionic-strength electroosmotic flows in uncharged hydrophobic nanochannels.

We report molecular dynamics simulation results of high-ionic-strength electroosmotic flows inside uncharged nanochannels. The possibility of this unusual electrokinetic phenomenon has been discussed by Dukhin et al. [A. Dukhin, S. Dukhin, P. Goetz, Langmuir 21 (2005) 9990]. Our computed velocity profiles clearly indicate the presence of a net flow with a maximum velocity around 2 m/s. We found the apparent zeta potential to be -29.7+/-6.8 mV, using the Helmholtz-Smoluchowski relation and the measured mean velocity. This value is comparable to experimentally measured values in Dukhin et al. and references therein. We also investigate the orientations of water molecules in response to an electric field by computing polarization density. Water molecules in the bulk region are oriented along the direction of the external electric field, while their near-wall orientation shows oscillations. The computation of three-dimensional density distributions of sodium and chloride ions around each individual water molecule show that chloride ions tend to concentrate near a water molecule, whereas sodium ions are diffusely distributed.

Adaptive biasing force method for scalar and vector free energy calculations.

In free energy calculations based on thermodynamic integration, it is necessary to compute the derivatives of the free energy as a function of one (scalar case) or several (vector case) order parameters. We derive in a compact way a general formulation for evaluating these derivatives as the average of a mean force acting on the order parameters, which involves first derivatives with respect to both Cartesian coordinates and time. This is in contrast with the previously derived formulas, which require first and second derivatives of the order parameter with respect to Cartesian coordinates. As illustrated in a concrete example, the main advantage of this new formulation is the simplicity of its use, especially for complicated order parameters. It is also straightforward to implement in a molecular dynamics code, as can be seen from the pseudocode given at the end. We further discuss how the approach based on time derivatives can be combined with the adaptive biasing force method, an enhanced sampling technique that rapidly yields uniform sampling of the order parameters, and by doing so greatly improves the efficiency of free energy calculations. Using the backbone dihedral angles Phi and Psi in N-acetylalanyl-N'-methylamide as a numerical example, we present a technique to reconstruct the free energy from its derivatives, a calculation that presents some difficulties in the vector case because of the statistical errors affecting the derivatives.

Molecular dynamics simulation of electro-osmotic flows in rough wall nanochannels.

We performed equilibrium and nonequilibrium molecular dynamics simulation to study electro-osmotic flows inside charged nanochannels with different types of surface roughness. We modeled surface roughness as a sequence of two-dimensional subnanoscale grooves and ridges (step function-type roughness) along the flow direction. The amplitude, spatial period, and symmetry of surface roughness were varied. The amplitude of surface roughness was on the order of the Debye length. The walls have uniform negative charges at the interface with fluids. We included only positive ions (counterions) for simplicity of computation. For the smooth wall, we compared our molecular dynamics simulation results to the well-known Poisson-Boltzmann theory. The density profiles of water molecules showed "layering" near the wall. For the rough walls, the density profiles measured from the wall are similar to those for the smooth wall except near where the steps are located. Because of the layering of water molecules and the finite size effect of ions and the walls, the ionic distribution departs from the Boltzmann distribution. To further understand the structure of water molecules and ions, we computed the polarization density. Near the wall, its z component dominates the other components, indicating the preferred orientation ("ordering") of water molecules. Especially, inside the groove for the rough walls, its maximum is 10% higher (stronger ordering) than for the smooth wall. The dielectric constant, computed with a Clausius-Mosotti-type equation, confirmed the ordering near the wall and the enhanced ordering inside the groove. The residence time and the diffusion coefficient, computed using the velocity autocorrelation function, showed that the diffusion of water and ions along the direction normal to the wall is significantly reduced near the wall and further decreases inside the groove. Along the flow direction, the diffusion of water and ions inside the groove is significantly lowered while it is similar to the bulk value elsewhere. We performed nonequilibrium molecular dynamics simulation to compute electro-osmotic velocities and flow rates. The velocity profiles correspond to those for overlapped electric double layers. For the rough walls, velocity inside the groove is close to zero, meaning that the channel height is effectively reduced. The flow rate was found to decrease as the period of surface roughness decreases or the amplitude of surface roughness increases. We defined the zeta potential as the electrostatic potential at the location of a slip plane. We computed the electrostatic potential with the ionic distribution and the dielectric constant both from our molecular dynamics simulation. We estimated the slip plane from the velocity profile. The zeta potential showed the same trend as the flow rate: it decreases with an increasing amplitude and a decreasing period of surface roughness.

A fast multipole method for Maxwell equations stable at all frequencies.

The solution of Helmholtz and Maxwell equations by integral formulations (kernel in exp(i kr)/r) leads to large dense linear systems. Using direct solvers requires large computational costs in O(N(3)). Using iterative solvers, the computational cost is reduced to large matrix-vector products. The fast multipole method provides a fast numerical way to compute convolution integrals. Its application to Maxwell and Helmholtz equations was initiated by Rokhlin, based on a multipole expansion of the interaction kernel. A second version, proposed by Chew, is based on a plane-wave expansion of the kernel. We propose a third approach, the stable-plane-wave expansion, which has a lower computational expense than the multipole expansion and does not have the accuracy and stability problems of the plane-wave expansion. The computational complexity is Nlog N as with the other methods.

Assessing the efficiency of free energy calculation methods.

The efficiencies of two recently developed methods for calculating free energy changes along a generalized coordinate in a system are discussed in the context of other, related approaches. One method is based on Jarzynski's identity [Phys. Rev. Lett. 78, 2690 (1997)]. The second method relies on thermodynamic integration of the average force and is called the adaptive biasing force method [Darve and Pohorille, J. Chem. Phys. 115, 9169 (2001)]. Both methods are designed such that the system evolves along the chosen coordinate(s) without experiencing free energy barriers and they require calculating the instantaneous, unconstrained force acting on this coordinate using the formula derived by Darve and Pohorille. Efficiencies are analyzed by comparing analytical estimates of statistical errors and by considering two numerical examples-internal rotation of hydrated 1,2-dichloroethane and transfer of fluoromethane across a water-hexane interface. The efficiencies of both methods are approximately equal in the first but not in the second case. During transfer of fluoromethane the system is easily driven away from equilibrium and, therefore, the performance of the method based on Jarzynski's identity is poor.