Brute-force Simulations Research Articles

Equation-free Model Reduction in Agent-based Computations: Coarse-grained Bifurcation and Variable-free Rare Event Analysis

We study the coarse-grained, reduced dynamics of an agent-based market model due to Omurtag and Sirovich []. We first describe the large agent number, deterministic limit of the system dynamics by performing numerical bifurcation calculations on a continuum approximation of their model. By exploring a broad parameter space, we observe several interesting phenomena including turning points leading to unstable stationary agent density distributions as well as a type of “termination point.” Close to these deterministic turning points we expect the stochastic underlying model to exhibit rare event transitions. We then proceed to discuss a coarse-grained approach to the quantitative study of these rare events. The basic assumption is that the dynamics of the system can be decomposed into fast (noise) and slow (single reaction coordinate) dynamics, so that the system can be described by an effective, coarse-grained Fokker-Planck(FP) equation. An explicit form of this effective FP equation is not available; in our computations we bypass the lack of a closed form equation by numerically estimating its components - the drift and diffusion coefficients - from ensembles of short bursts of microscopic simulations with judiciously chosen initial conditions. The reaction coordinate is first constructed based on our understanding of the continuum model close to the turning points, and it gives results reasonably close to those from brute-force direct simulations. When no guidelines for the selection of a good reaction coordinate are available, data-mining tools, in particular Diffusion Maps, can be used to determine a suitable reaction coordinate. In the third part of this work we demonstrate this “variable-free” approach by constructing a reaction coordinate simply based on the data from the simulation itself. This Diffusion Map based, empirical coordinate gives results consistent with the direct simulation.

Risk Analysis for Large Pools of Loans

Financial institutions, government-sponsored enterprises, and asset-backed security investors are often exposed to delinquency and prepayment risk from large numbers of loans. Examples include mortgages, credit cards, and auto, student, and business loans. Because of the size of the pools, the measurement of these exposures is computationally expensive. This paper develops and tests efficient numerical methods for the analysis of large pools of loans as well as asset-backed securities backed by such pools. For a broad class of statistical and machine learning models of loan-level risk in discrete time, we prove a law of large numbers and a central limit theorem for the pool-level risk. The asymptotics are then used to construct computationally efficient approximations for a large pool. The approximations aggregate the full loan-level dynamics, making it possible to take advantage of the detailed loan-level data often available in practice. We furthermore prove a convergence rate for the approximation and a uniform integrability result. The latter allows the approximation to be applied to a wide class of pool-level risk statistics. To demonstrate the effectiveness of our approach, we implement it on a data set of over 25 million mortgages. The results show the accuracy and speed of the approximation in comparison to brute-force simulation of a pool, for a variety of pools with different risk profiles. The online appendix is available at https://doi.org/10.1287/mnsc.2017.2947 . This paper was accepted by Noah Gans, stochastic models.

Nonconservative dynamics in long atomic wires

The effect of nonconservative current-induced forces on the ions in a defect-free metallic nanowire is investigated using both steady-state calculations and dynamical simulations. Nonconservative forces were found to have a major influence on the ion dynamics in these systems, but their role in increasing the kinetic energy of the ions decreases with increasing system length. The results illustrate the importance of nonconservative effects in short nanowires and the scaling of these effects with system size. The dependence on bias and ion mass can be understood with the help of a simple pen and paper model. This material highlights the benefit of simple preliminary steady-state calculations in anticipating aspects of brute-force dynamical simulations, and provides rule of thumb criteria for the design of stable quantum wires.

A reduced model for domain walls in soft ferromagnetic films at the cross-over from symmetric to asymmetric wall types

We study the Landau-Lifshitz model for the energy of multi-scale transition layers – called "domain walls'' – in soft ferromagnetic films. Domain walls separate domains of constant magnetization vectors m^\pm_\alpha\in\mathbb{S}^2 that differ by an angle 2\alpha . Assuming translation invariance tangential to the wall, our main result is the rigorous derivation of a reduced model for the energy of the optimal transition layer, which in a certain parameter regime confirms the experimental, numerical and physical predictions: The minimal energy splits into a contribution from an asymmetric, divergence-free core which performs a partial rotation in \mathbb{S}^2 by an angle 2\theta , and a contribution from two symmetric, logarithmically decaying tails, each of which completes the rotation from angle \theta to \alpha in \mathbb{S}^1 . The angle \theta is chosen such that the total energy is minimal. The contribution from the symmetric tails is known explicitly, while the contribution from the asymmetric core is analyzed in [7]. Our reduced model is the starting point for the analysis of a bifurcation phenomenon from symmetric to asymmetric domain walls. Moreover, it allows for capturing asymmetric domain walls including their extended tails (which were previously inaccessible to brute-force numerical simulation).

Effects of surface interactions on heterogeneous ice nucleation for a monatomic water model.

Despite its importance in atmospheric science, much remains unknown about the microscopic mechanism of heterogeneous ice nucleation. In this work, we perform hybrid Monte Carlo simulations of the heterogeneous nucleation of ice on a range of generic surfaces, both flat and structured, in order to probe the underlying factors affecting the nucleation process. The structured surfaces we study comprise one basal plane bilayer of ice with varying lattice parameters and interaction strengths. We show that what determines the propensity for nucleation is not just the surface attraction, but also the orientational ordering imposed on liquid water near a surface. In particular, varying the ratio of the surface's attraction and orientational ordering can change the mechanism by which nucleation occurs: ice can nucleate on the structured surface even when the orientational ordering imposed by the surface is weak, as the water molecules that interact strongly with the surface are themselves a good template for further growth. We also show that lattice matching is important for heterogeneous nucleation on the structured surface we study. We rationalise these brute-force simulation results by explicitly calculating the interfacial free energies of ice and liquid water in contact with the nucleating surface and their variation with surface interaction parameters.

Stem cell differentiation as a many-body problem.

Stem cell differentiation has been viewed as coming from transitions between attractors on an epigenetic landscape that governs the dynamics of a regulatory network involving many genes. Rigorous definition of such a landscape is made possible by the realization that gene regulation is stochastic, owing to the small copy number of the transcription factors that regulate gene expression and because of the single-molecule nature of the gene itself. We develop an approximation that allows the quantitative construction of the epigenetic landscape for large realistic model networks. Applying this approach to the network for embryonic stem cell development explains many experimental observations, including the heterogeneous distribution of the transcription factor Nanog and its role in safeguarding the stem cell pluripotency, which can be understood by finding stable steady-state attractors and the most probable transition paths between those attractors. We also demonstrate that the switching rate between attractors can be significantly influenced by the gene expression noise arising from the fluctuations of DNA occupancy when binding to a specific DNA site is slow.

Selectivity trend of gas separation through nanoporous graphene

By means of molecular dynamics (MD) simulations, we demonstrate that porous graphene can efficiently separate gases according to their molecular sizes. The flux sequence from the classical MD simulation is H2>CO2≫N2>Ar>CH4, which generally follows the trend in the kinetic diameters. This trend is also confirmed from the fluxes based on the computed free energy barriers for gas permeation using the umbrella sampling method and kinetic theory of gases. Both brute-force MD simulations and free-energy calcualtions lead to the flux trend consistent with experiments. Case studies of two compositions of CO2/N2 mixtures further demonstrate the separation capability of nanoporous graphene.

Fast Approximation of Loan Portfolio Loss

Economic capital and Value at Risk are key measures in risk management and risk regulations. The industry practice is to calculate these risk measures based on brute force simulations using market observed index returns as risk factors and some assumed correlation structure. In this paper, we illustrate how the underlying risk factors can be constructed from time series of default probabilities and how Value at Risk and Expected Shortfall can be calculated using a closed form approximation for a multi-factor model. Furthermore, economic capital figures based on risk factors extracted from different sources are compared and analyzed using a large portfolio of over 500 obligors and a small portfolio of 26 obligors. Our findings indicate that observed index returns may not explain default variations well and could lead to an underestimation of risk. Risk neutral default probabilities capture market conditions and changes in market sentiment well, but could lead to an over estimation of risk due to the embedded risk premium.

Conformational Transitions in ATP-Driven Calcium Pump SERCA

Sarco/endoplasmic reticulum Ca2+-ATPase (SERCA) is an integral membrane protein that uses ATP hydrolysis as a source of free energy to pump two calcium ions per ATP molecule from calcium poor cytoplasm of the muscle cell to the calcium rich lumen of the sarcoplasmic reticulum, thereby maintaining a ten thousand fold concentration gradient. Detailed structural studies of the pump under different conditions provided analogues of various intermediates in the reaction cycle and revealed important changes in the tertiary structure of the protein both in the cytoplasmic and in the transmembrane parts. Two major outstanding issues are the pathways of the ions to and from the transmembrane binding sites and a detailed understanding of the large scale conformational changes among various functionally relevant states. We have applied all-atom molecular dynamics (MD) and string method with swarms-of-trajectories to study transition pathways among various experimental structures. The size of the system (291,000) and the time-scale involved in the large scale motions prohibit the use of brute-force unbiased simulations to obtain statistically meaningful information. To circumvent this challenge, the string method with swarms-of-trajectories is used to discover the optimal minimum free energy path between the two end-states. The path, which is represented as a chain of states or images in the space of relevant collective variables, is optimized by iterating two steps: moving each image along the drift calculated from an ensemble (swarm) of short unbiased MD trajectories initiated from the image and a reparametrization procedure that keeps all the images equidistant. We hope to understand the molecular details of the complex motions involved in the conformational transitions by studying the pathways produced by the string method and correlating these observations to available structural and biochemical experiments.

Efficient Risk Analysis of Mortgage Pools

Financial institutions, government-sponsored enterprises such as Fannie Mae, and asset-backed security investors are often exposed to delinquency and prepayment risk from large numbers of loans. Examples include mortgages, credit cards, auto, commercial, real estate, student, and small business loans. Due to the size of such loan pools and the potentially long maturities of their cashflows, the measurement and management of these exposures is computationally expensive. This paper develops and tests efficient numerical methods for the risk analysis of large pools of loans. For a broad class of dynamic loan-level models of delinquency and prepayment, we develop a law of large numbers and a central limit theorem for the loss and prepayment levels in the pool. The asymptotics are then used to construct efficient Monte Carlo approximations of the loss and prepayment distributions for a large pool. The approximations aggregate the full loan-level dynamics, making it possible to take advantage of the detailed loan-level data often available. The Monte Carlo approximation allows for efficient risk analysis of loan portfolios as well as MBSs, CMOs, and other ABSs backed by pools of loans. To demonstrate the effectiveness of our approach, we implement it on a data set of over 25 million actual subprime and agency mortgages. The results show the accuracy and speed of the approximation in comparison to brute-force simulation of a pool, for a variety of pools with different risk characteristics. Computational cost is often several orders of magnitude less than brute-force simulation of the actual pool with a similar level of accuracy. Furthermore, the computational expense of our efficient Monte Carlo approximation is constant no matter the dimension of the loan-level features; this is key since loan-level feature data is often high-dimensional.

Rare-event sampling for driven inertial systems via the nonequilibrium distribution function

Recently substantial generalizations have been made to the Yamada-Kawasaki type nonlinear response theory formalism for deterministic dynamics. These provide opportunities to develop new importance sampling techniques against known distribution functions. We exploit this to develop a method of calculating nonequilibrium rate constants across a large free energy barrier, allowing molecular simulations to access far greater time scales. In contrast to existing stochastic methods, we do not merely follow the dynamics forward in time, allowing systems governed by inertial equations of motion to be usefully addressed. We further generalize the Yamada-Kawasaki type formalism to the case of stochastic equations of motion and the resulting nonequilibrium importance sampling method follows straight forwardly. We quantitatively test the method on stochastic and deterministic models, specifically chosen to allow a comparison with reliable, independently obtained data. To do this we consider oscillators in a one-dimensional double well potential. Either side of the well is separated by an energy barrier and the oscillator is driven away from equilibrium by a position dependent temperature profile. The energy barrier is chosen to be low enough that we can still solve the system computationally by brute force. We demonstrate that the results of our method are consistent with those from the brute force simulations to very high accuracy. We also show how our method is more than two orders of magnitude more efficient for the given example than forward flux sampling.

A method to represent ozone response to large changes in precursor emissions using high-order sensitivity analysis in photochemical models

Abstract. Photochemical grid models (PGMs) are used to simulate tropospheric ozone and quantify its response to emission changes. PGMs are often applied for annual simulations to provide both maximum concentrations for assessing compliance with air quality standards and frequency distributions for assessing human exposure. Efficient methods for computing ozone at different emission levels can improve the quality of ozone air quality management efforts. This study demonstrates the feasibility of using the decoupled direct method (DDM) to calculate first- and second-order sensitivity of ozone to anthropogenic NOx and VOC emissions in annual PGM simulations at continental scale. Algebraic models are developed that use Taylor series to produce complete annual frequency distributions of hourly ozone at any location and any anthropogenic emission level between zero and 100%, adjusted independently for NOx and VOC. We recommend computing the sensitivity coefficients at the midpoint of the emissions range over which they are intended to be applied, in this case with 50% anthropogenic emissions. The algebraic model predictions can be improved by combining sensitivity coefficients computed at 10 and 50% anthropogenic emissions. Compared to brute force simulations, algebraic model predictions tend to be more accurate in summer than winter, at rural than urban locations, and with 100% than zero anthropogenic emissions. Equations developed to combine sensitivity coefficients computed with 10 and 50% anthropogenic emissions are able to reproduce brute force simulation results with zero and 100% anthropogenic emissions with a mean bias of less than 2 ppb and mean error of less than 3 ppb averaged over 22 US cities.

Computation of transit times using the milestoning method with applications to polymer translocation

Milestoning is an efficient approximation for computing long-time kinetics and thermodynamics of large molecular systems, which are inaccessible to brute-force molecular dynamics simulations. A common use of milestoning is to compute the mean first passage time (MFPT) for a conformational transition of interest. However, the MFPT is not always the experimentally observed timescale. In particular, the duration of the transition path, or the mean transit time, can be measured in single-molecule experiments, such as studies of polymers translocating through pores and fluorescence resonance energy transfer studies of protein folding. Here we show how to use milestoning to compute transit times and illustrate our approach by applying it to the translocation of a polymer through a narrow pore.

A multiscale algorithm for simulating spatially-extended Langmuir circulation dynamics

A multiscale algorithm is proposed for simulations of the spatially-extended dynamics of Langmuir circulation, a wind- and surface-wave-driven convective flow that dominates vertical transport and mixing on the scale of the O(100)-m deep ocean surface boundary layer. The algorithm is motivated by multiple scale asymptotic analysis of the master partial differential (Navier–Stokes-like) equations and comprises coupled equation sets for the fine- and coarse-scale dynamics. A primary virtue of the resulting multiscale formulation is that dynamical teleconnections between the fine and coarse scales are explicitly identified and computed, enabling fewer fine-scale domains than coarse-scale grid points to be used. Quantitative comparisons between the multiscale simulations and direct numerical simulations of the master partial differential equations are made. Good agreement is demonstrated in the appropriate parameter regime at a cost that is orders of magnitude less than that required for brute-force simulations using the single-scale algorithm. Several novel multiscale physical interactions are identified, highlighting the utility of the algorithm and the potential importance of spatiotemporal modulation of the convective flow on the mean dynamics.

A Site Density Functional Theory for Water: Application to Solvation of Amino Acid Side Chains.

We report a site density functional theory (SDFT) based on the conventional atomistic models of water and the universality ansatz of the bridge functional. The excess Helmholtz energy functional is formulated in terms of a quadratic expansion with respect to the local density deviation from that of a uniform system and a universal functional for all higher-order terms approximated by that of a reference hard-sphere system. With the atomistic pair direct correlation functions of the uniform system calculated from MD simulation and an analytical expression for the bridge functional from the modified fundamental measure theory, the SDFT can be used to predict the structure and thermodynamic properties of water under inhomogeneous conditions with a computational cost negligible in comparison to that of brute-force simulations. The numerical performance of the SDFT has been demonstrated with the predictions of the solvation free energies of 15 molecular analogs of amino acid side chains in water represented by SPC/E, SPC, and TIP3P models. For theTIP3P model, a comparison of the theoretical predictions with MD simulation and experimental data shows agreement within 0.64 and 1.09 kcal/mol on average, respectively.

Smoking epidemic eradication in a eco-epidemiological dynamical model

Smoking is perceived as a major epidemic with regard to mortality. Modelling is a major tool used to obtain insight in the dynamics and possible solutions to decrease or even eradicate this epidemic. Most models on smoking consider the epidemiological context explicitly, in which smoking is regarded as an ‘infectious disease’, in which individuals ‘infect’ each other. However, the population dynamics are often ignored, while these occur at roughly the same timescale as smoking, and hence should explicitly be considered in the modelling of smoking. We present a simple but dynamical eco-epidemiological model. The model formulation consists of a resource-population dynamic part coupled to an epidemiological part resembling a SIR type model for the three compartments: non-smokers, smokers and ex-smokers. The coupling is via birth of non-smokers and death of the three classes with different death rates. The final four-dimensional system of ordinary differential equations are studied using brute force simulations for the short term dynamics and bifurcation analysis for the long-term dynamics. Due to a feed-back mechanism of the two coupling terms there is a codim-two tangent-transcritical bifurcation. This leads to bi-stability of one smoker endemic interior equilibrium and a smoker free boundary equilibrium. Changing parameters beyond the emerging tangent bifurcation leads on the short term to eradicating smoking. We consider the Netherlands in this paper for parametrization, but the modelling approach may be generally applicable.

Statistical Physics Analysis of Maximum a Posteriori Estimation for Multi-channel Hidden Markov Models

The performance of Maximum a posteriori (MAP) estimation is studied analytically for binary symmetric multi-channel Hidden Markov processes. We reduce the estimation problem to a 1D Ising spin model and define order parameters that correspond to different characteristics of the MAP-estimated sequence. The solution to the MAP estimation problem has different operational regimes separated by first order phase transitions. The transition points for $L$-channel system with identical noise levels, are uniquely determined by $L$ being odd or even, irrespective of the actual number of channels. We demonstrate that for lower noise intensities, the number of solutions is uniquely determined for odd $L$, whereas for even $L$ there are exponentially many solutions. We also develop a semi analytical approach to calculate the estimation error without resorting to brute force simulations. Finally, we examine the tradeoff between a system with single low-noise channel and one with multiple noisy channels.

Simulating rare events using a weighted ensemble-based string method.

We introduce an extension to the weighted ensemble (WE) path sampling method to restrict sampling to a one-dimensional path through a high dimensional phase space. Our method, which is based on the finite-temperature string method, permits efficient sampling of both equilibrium and non-equilibrium systems. Sampling obtained from the WE method guides the adaptive refinement of a Voronoi tessellation of order parameter space, whose generating points, upon convergence, coincide with the principle reaction pathway. We demonstrate the application of this method to several simple, two-dimensional models of driven Brownian motion and to the conformational change of the nitrogen regulatory protein C receiver domain using an elastic network model. The simplicity of the two-dimensional models allows us to directly compare the efficiency of the WE method to conventional brute force simulations and other path sampling algorithms, while the example of protein conformational change demonstrates how the method can be used to efficiently study transitions in the space of many collective variables.

Conformation-Switching in Adenylate Kinase Revisited with a Path-Ensemble Simulation

Adenylate kinase (ADK), which reversibly converts ATP and AMP to two ADP molecules, has two conformational states, inactive (open) and active (closed); crystal structures of both states were solved in mid 1990s. Numerous studies using experiments as well as computer simulations have aimed to elucidate the relationship between conformational rearrangement and ADK function. However, the detailed mechanism of the switching between the states has not been fully understood. Computationally, sampling the conformation space is the key limitation. We therefore employ the “weighted ensemble” (WE) simulation technique, a parallel, enhanced sampling method for characterizing non-equilibrium and equilibrium processes that would be impractical to observe using brute force simulation. Because it yields an ensemble of transition events, WE provides a unique lens for examining the mechanism (i.e., the time sequence of intermediate states), including possible heterogeneity, of conformational switching in all-atom ADK.

Simultaneous Calculation of Dynamical and Equilibrium Properties of Atomistic Systems using an Ensemble of Weighted Trajectories

Equilibrium can be represented as an ensemble of trajectories undergoing dynamics in which detailed balance is maintained. Here, we present results for a simulation protocol based on the “weighted ensemble” (WE) approach of Huber and Kim in which weak coupling between parallel simulations can capture both dynamical and equilibrium properties without statistical bias. The approach is demonstrated in two molecular systems: alanine di-peptide and chignolin, a 10 residue mini-protein. We present results for the mean first passage time between arbitrary states and demonstrate the rapid convergence towards equilibrium for chignolin in comparison to a brute force simulation. Importantly, the analysis of weighted ensemble trajectories does not require states to be defined in advance, and therefore permits the calculation of rates between arbitrary states selected after the simulation is run. We also describe a history-dependent Markov-like formulation to correct for a potentially significant bias in kinetic properties measured from equilibrium WE simulations.