High-Dimensional Parameter Search Method to Determine Force Field Mixing Terms in Molecular Simulations.

Molecular dynamics (MD) force fields for lipids and ions are typically developed independently of one another. In simulations consisting of both lipids and ions, lipid-ion interaction energies are estimated using a predefined set of mixing rules for Lennard-Jones (LJ) interactions. This, however, does not guarantee their reliability. In fact, compared to the quantum mechanical reference data, Lorentz-Berthelot mixing rules substantially underestimate the binding energies of Na+ ions with small-molecule analogues of lipid headgroups, yielding errors on the order of 80 and 130 kJ/mol, respectively, for methyl acetate and diethyl phosphate. Previously, errors associated with mixing force fields have been reduced using approaches such as "NB-fix" in which LJ interactions are computed using explicit cross terms rather than those from mixing rules. Building on this idea, we derive explicit lipid-ion cross terms that also may implicitly include many-body cooperativity effects. Additionally, to account for the interdependency between cross terms, we optimize all cross terms simultaneously by performing high-dimensional searches using our ParOpt software. The cross terms we obtain reduce the errors due to mixing rules to below 10 kJ/mol. MD simulation of the lipid bilayer conducted using these optimized cross terms resolves the structural discrepancies between our previous simulations and small-angle X-ray and neutron scattering experiments. These results demonstrate that simulations of lipid bilayers with ions that are accurate up to structural data from scattering experiments can be performed without explicit polarization terms. However, it is worth noting that such NB-fix cross terms are not based on any physical principle; a polarizable lipid model would be more realistic and is still desired. Our approach is generic and can be applied to improve the accuracies of simulations employing mixed force fields.

Molecular basis for higher affinity of SARS-CoV-2 spike RBD for human ACE2 receptor.

Severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) has caused substantially more infections, deaths, and economic disruptions than the 2002-2003 SARS-CoV. The key to understanding SARS-CoV-2's higher infectivity lies partly in its host receptor recognition mechanism. Experiments show that the human angiotensin converting enzyme 2 (ACE2) protein, which serves as the primary receptor for both CoVs, binds to the receptor binding domain (RBD) of CoV-2's spike protein stronger than SARS-CoV's spike RBD. The molecular basis for this difference in binding affinity, however, remains unexplained from X-ray structures. To go beyond insights gained from X-ray structures and investigate the role of thermal fluctuations in structure, we employ all-atom molecular dynamics simulations. Microseconds-long simulations reveal that while CoV and CoV-2 spike-ACE2 interfaces have similar conformational binding modes, CoV-2 spike interacts with ACE2 via a larger combinatorics of polar contacts, and on average, makes 45% more polar contacts. Correlation analysis and thermodynamic calculations indicate that these differences in the density and dynamics of polar contacts arise from differences in spatial arrangements of interfacial residues, and dynamical coupling between interfacial and non-interfacial residues. These results recommend that ongoing efforts to design spike-ACE2 peptide blockers will benefit from incorporating dynamical information as well as allosteric coupling effects.

Transferable interactions of Li+ and Mg2+ ions in polarizable models.

Therapeutic implications of Li+, in many cases, stem from its ability to inhibit certain Mg2+-dependent enzymes, where it interacts with or substitutes for Mg2+. The underlying details of its action are, however, unknown. Molecular simulations can provide insights, but their reliability depends on how well they describe relative interactions of Li+ and Mg2+ with water and other biochemical groups. Here, we explore, benchmark, and recommend improvements to two simulation approaches: the one that employs an all-atom polarizable molecular mechanics (MM) model and the other that uses a hybrid quantum and MM implementation of the quasi-chemical theory (QCT). The strength of the former is that it describes thermal motions explicitly and that of the latter is that it derives local contributions from electron densities. Reference data are taken from the experiment, and also obtained systematically from CCSD(T) theory, followed by a benchmarked vdW-inclusive density functional theory. We find that the QCT model predicts relative hydration energies and structures in agreement with the experiment and without the need for additional parameterization. This implies that accurate descriptions of local interactions are essential. Consistent with this observation, recalibration of local interactions in the MM model, which reduces errors from 10.0 kcal/mol to 1.4 kcal/mol, also fixes aqueous phase properties. Finally, we show that ion-ligand transferability errors in the MM model can be reduced significantly from 10.3 kcal/mol to 1.2 kcal/mol by correcting the ligand's polarization term and by introducing Lennard-Jones cross-terms. In general, this work sets up systematic approaches to evaluate and improve molecular models of ions binding to proteins.

Interactions of Monovalent and Divalent Cations at Palmitoyl-Oleoyl-Phosphatidylcholine Interface.

Li+ is a biologically active and medically important cation. Experiments show that Li+ modulates some phospholipid bilayer properties in a manner similar to divalent cations, rather than other monovalent cations. We previously performed a comparative simulation study of the interaction of several monovalent cations with palmitoyl-oleoyl-phosphatidylcholine bilayers and reported that Li+ exhibited the highest association with lipids and formed a unique tetrahedral coordinated structure with lipid head groups. Here we extend these studies to two biologically important divalent cations, Mg2+ and Ca2+, and observe that, just like monovalent cations, Mg2+ and Ca2+ reduce bilayer areas and increase chain order. Bilayer area changes induced by cations are strongly correlated with the amount of charge inside the headgroup region; however, Mg2+ and Li+ are clear outliers. At the same time though, Mg2+ adsorption in the bilayer is the smallest among all cations, which is in contrast to Li+ that binds strongly to lipids. In fact, in contrast to all other cations, Mg2+ remains fully hydrated in the lipid headgroup region. However, Li+ and Mg2+ share high overlap between their inner-shell coordination topologies. This suggests that Li+ can structurally replace Mg2+, which is bound to other biomolecules with up to fourfold coordination, provided such replacement is energetically feasible. We compute structural topologies and compare them quantitatively using a new weighted-graphs-based method. Finally, we find that the specificity of cation interaction with lipid head groups exhibit consistent trend with the solvation shell energetics of ions in lipid headgroup and bulk water regions.

Molecular dynamics simulations of ether- and ester-linked phospholipids.

Dissimilarities in the bulk structure of bilayers composed of ether- vs ester-linked lipids are well-established; however, the atomistic interactions responsible for these differences are not well known. These differences are important in understanding of why archaea have a different bilayer composition than the other domains of life and why humans have larger concentrations of plasmalogens in specialized membranes? In this paper, we simulate two lipid bilayers, the ester linked dipalmitoylphosphatidylcholine (DPPC) and the ether lined dihexadecylphosphatidylcholine (DHPC), to study these variations. The structural analysis of the bilayers reveals that DPPC is more compressible than DHPC. A closer examination of dipole potential shows DHPC, despite having a smaller dipole potential of the bilayer, has a higher potential barrier than DPPC at the surface. Analysis of water order and dynamics suggests DHPC has a more ordered, less mobile layer of water in the headgroup. These results seem to resolve the issue as to whether the decrease in permeability of DHPC is due to of differences in minimum area per lipid (A0) or diffusion coefficient of water in the headgroup region (Dhead) (Guler et al., 2009) since we have shown significant changes in the order and mobility of water in that region.

Effects of Lithium and Other Monovalent Ions on Palmitoyl Oleoyl Phosphatidylcholine Bilayer.

Interactions of monovalent salts with lipid membranes are explored with molecular dynamics (MD) simulations. The simulations included the monovalent ions Na+ and K+, for their importance in physiology, Li+ for its small size and importance in several medical conditions including bipolar disorder, and Rb+ for its large size. All simulations included Cl- as counterions. One bilayer was simulated without salt as a control. Palmitoyl oleoyl phosphatidylcholine (POPC) bilayers experienced reductions in area per lipid with the addition of salt; the smaller the ion the smaller the area, with the exception of Li+. Li+ exhibited unique binding affinities between phosphates and sn-2 carbonyls that lowered the order of the top part of sn-2 chain, which increased the area per lipid, compared to other ionic simulations. Further, we observe that monovalent salts alter bilayer properties through structural changes and not so much through the changes in surface potential.

Atomically detailed lipid bilayer models for the interpretation of small angle neutron and X-ray scattering data.

We present a new atom density profile (ADP) model and a statistical approach for extracting structural characteristics of lipid bilayers from X-ray and neutron scattering data. Models for five lipids with varying head and tail chemical composition in the fluid phase, 1,2-dioleoyl-sn-glycero-3-phosphatidylcholine (DOPC), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine (POPC), 1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine (DPPC), 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylserine (POPS), and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylglycerol (POPG), are optimized using a simplex based method to simultaneously reproduce both neutron and X-ray scattering data. Structural properties are determined using statistical analysis of multiple optimal model structures. The method and models presented make minimal assumptions regarding the atomic configuration, while taking into account the underlying physical properties of the system. The more general model and statistical approach yield data with well defined uncertainties, indicating the precision in determining density profiles, atomic locations, and bilayer structural characteristics. Resulting bilayer structures include regions exhibiting large conformational variation. Due to the increased detail in the model, the results demonstrate the possibility of a distinct hydration layer within the interfacial (backbone) region.

Why do chimpanzee males attack the females of neighboring communities?

Our closest nonhuman primate relatives, chimpanzees, engage in potentially lethal between-group conflict; this collective aggressive behavior shows parallels with human warfare. In some communities, chimpanzee males also severely attack and even kill females of the neighboring groups. This is surprising given their system of resource defense polygyny, where males are expected to acquire potential mates. We develop a simple mathematical model based on reproductive skew among primate males to solve this puzzle. The model predicts that it is advantageous for high-ranking males but not for low-ranking males to attack females. It also predicts that more males gain a benefit from attacking females as the community's reproductive skew decreases, i.e., as mating success is more evenly distributed. Thus, fatal attacks on females should be concentrated in communities with low reproductive skew. These attacks should also concur with between-community infanticide. A review of the chimpanzee literature provides enough preliminary support for this prediction to warrant more detailed testing.

Mating strategies in primates: a game theoretical approach to infanticide.

Infanticide by newly immigrated or newly dominant males is reported among a variety of taxa, such as birds, rodents, carnivores and primates. Here we present a game theoretical model to explain the presence and prevalence of infanticide in primate groups. We have formulated a three-player game involving two males and one female and show that the strategies of infanticide on the males' part and polyandrous mating on the females' part emerge as Nash equilibria that are stable under certain conditions. Moreover, we have identified all the Nash equilibria of the game and arranged them in a novel hierarchical scheme. Only in the subspace spanned by the males are the Nash equilibria found to be strict, and hence evolutionarily stable. We have therefore proposed a selection mechanism informed by adaptive dynamics to permit the females to transition to, and remain in, optimal equilibria after successive generations. Our model concludes that polyandrous mating by females is an optimal strategy for the females that minimizes infanticide and that infanticide confers advantage to the males only in certain regions of parameter space. We have shown that infanticide occurs during turbulent changes accompanying male immigration into the group. For changes in the dominance hierarchy within the group, we have shown that infanticide occurs only in primate groups where the chance for the killer to sire the next infant is high. These conclusions are confirmed by observations in the wild. This model thus has enabled us to pinpoint the fundamental processes behind the reproductive decisions of the players involved, which was not possible using earlier theoretical studies.

Multiscale simulations of heterogeneous model membranes.

This review will focus on computer modeling aimed at providing insights into the existence, structure, size, and thermodynamic stability of localized domains in membranes of heterogeneous composition. Modeling the lateral organization within a membrane is problematic due to the relatively slow lateral diffusion rate for lipid molecules so that microsecond or longer time scales are needed to fully model the formation and stability of a raft in a membrane. Although atomistic simulations currently are not able to reach this scale, they can provide data on the intermolecular forces and correlations that are involved in lateral organization. These data can be used to define coarse grained models that are capable of predictions of lateral organization in membranes. In this paper, we review modeling efforts that use interaction data from MD simulations to construct coarse grained models for heterogeneous bilayers. In this review we will discuss MD simulations done with the aim of gaining the information needed to build accurate coarse-grained models. We will then review some of the coarse-graining work, emphasizing modeling that has resulted from or has a basis in atomistic simulations.

Lateral organization of complex lipid mixtures from multiscale modeling.

The organizational properties of complex lipid mixtures can give rise to functionally important structures in cell membranes. In model membranes, ternary lipid-cholesterol (CHOL) mixtures are often used as representative systems to investigate the formation and stabilization of localized structural domains ("rafts"). In this work, we describe a self-consistent mean-field model that builds on molecular dynamics simulations to incorporate multiple lipid components and to investigate the lateral organization of such mixtures. The model predictions reveal regions of bimodal order on ternary plots that are in good agreement with experiment. Specifically, we have applied the model to ternary mixtures composed of dioleoylphosphatidylcholine:18:0 sphingomyelin:CHOL. This work provides insight into the specific intermolecular interactions that drive the formation of localized domains in these mixtures. The model makes use of molecular dynamics simulations to extract interaction parameters and to provide chain configuration order parameter libraries.

A reactive molecular dynamics simulation of the silica-water interface.

We report our study of a silica-water interface using reactive molecular dynamics. This first-of-its-kind simulation achieves length and time scales required to investigate the detailed chemistry of the system. Our molecular dynamics approach is based on the ReaxFF force field of van Duin et al. [J. Phys. Chem. A 107, 3803 (2003)]. The specific ReaxFF implementation (SERIALREAX) and force fields are first validated on structural properties of pure silica and water systems. Chemical reactions between reactive water and dangling bonds on a freshly cut silica surface are analyzed by studying changing chemical composition at the interface. In our simulations, reactions involving silanol groups reach chemical equilibrium in approximately 250 ps. It is observed that water molecules penetrate a silica film through a proton-transfer process we call "hydrogen hopping," which is similar to the Grotthuss mechanism. In this process, hydrogen atoms pass through the film by associating and dissociating with oxygen atoms within bulk silica, as opposed to diffusion of intact water molecules. The effective diffusion constant for this process, taken to be that of hydrogen atoms within silica, is calculated to be 1.68 x 10(-6) cm(2)/s. Polarization of water molecules in proximity of the silica surface is also observed. The subsequent alignment of dipoles leads to an electric potential difference of approximately 10.5 V between the silica slab and water.

Why Class Formation Occurs in Humans but Not among Other Primates : A Primate Coalitions Model.

Most human societies exhibit a distinct class structure, with an elite, middle classes, and a bottom class, whereas animals form simple dominance hierarchies in which individuals with higher fighting ability do not appear to form coalitions to "oppress" weaker individuals. Here, we extend our model of primate coalitions and find that a division into a bottom class and an upper class is inevitable whenever fitness-enhancing resources, such as food or real estate, are exploitable or tradable and the members of the bottom class cannot easily leave the group. The model predicts that the bottom class has a near flat, low payoff and always comprises at least half the society. The upper class may subdivide into one or more middle class(es), resulting in improved payoff for the topmost members (elite). The model predicts that the bottom class on its own is incapable of mounting effective counter-coalitions against the upper class, except when receiving support from dissatisfied members of the middle class(es). Such counter-coalitions can be prevented by keeping the payoff to the lowest-ranked members of the middle classes (through concessions) well above that of the bottom class. This simple model explains why classes are also absent in nomadic hunter-gatherers and predominate in (though are not limited to) societies that produce and store food. Its results also agree well with various other known features of societies with classes.

An improved united atom force field for simulation of mixed lipid bilayers.

We introduce a new force field (43A1-S3) for simulation of membranes by the Gromacs simulation package. Construction of the force fields is by standard methods of electronic structure computations for bond parameters and charge distribution and specific volumes and heats of vaporization for small-molecule components of the larger lipid molecules for van der Waals parameters. Some parameters from the earlier 43A1 force field are found to be correct in the context of these calculations, while others are modified. The validity of the force fields is demonstrated by correct replication of X-ray form factors and NMR order parameters over a wide range of membrane compositions in semi-isotropic NTP 1 atm simulations. 43-A1-S3 compares favorably with other force fields used in conjunction with the Gromacs simulation package with respect to the breadth of phenomena that it accurately reproduces.

Interaction of salt with ether- and ester-linked phospholipid bilayers.

A distinguishing feature of Archaeal plasma membranes is that their phospholipids contain ether-links, as opposed to bacterial and eukaryotic plasma membranes where phospholipids primarily contain ester-links. Experiments show that this chemical difference in headgroup-tail linkage does produce distinct differences in model bilayer properties. Here we examine the effects of salt on bilayer structure in the case of an ether-linked lipid bilayer. We use molecular dynamics simulations and compare equilibrium properties of two model lipid bilayers in NaCl salt solution - POPC and its ether-linked analog that we refer to as HOPC. We make the following key observations. The headgroup region of HOPC "adsorbs" fewer ions compared to the headgroup region of POPC. Consistent with this, we note that the Debye screening length in the HOPC system is ∼ 10% shorter than that in the POPC system. Herein, we introduce a protocol to identify the lipid-water interfacial boundary that reproduces the bulk salt distribution consistent with Gouy-Chapman theory. We also note that the HOPC bilayer has excess solvent in the headgroup region when compared to POPC, coinciding with a trough in the electrostatic potential. Waters in this region have longer autocorrelation times and smaller lateral diffusion rates compared to the corresponding region in the POPC bilayer, suggesting that the waters in HOPC are more strongly coordinated to the lipid headgroups. Furthermore, we note that it is this region of tightly coordinated waters in the HOPC system that has a lower density of Na+ ions. Based on these observations we conclude that an ether-linked lipid bilayer has a lower binding affinity for Na+ compared to an ester-linked lipid bilayer.

Atomistic and coarse-grained computer simulations of raft-like lipid mixtures.

Computer modeling can provide insights into the existence, structure, size, and thermodynamic stability of localized raft-like regions in membranes. However, the challenges in the construction and simulation of accurate models of heterogeneous membranes are great. The primary obstacle in modeling the lateral organization within a membrane is the relatively slow lateral diffusion rate for lipid molecules. Microsecond or longer time-scales are needed to fully model the formation and stability of a raft in a membra ne. Atomistic simulations currently are not able to reach this scale, but they do provide quantitative information on the intermolecular forces and correlations that are involved in lateral organization. In this chapter, the steps needed to carry out and analyze atomistic simulations of hydrated lipid bilayers having heterogeneous composition are outlined. It is then shown how the data from a molecular dynamics simulation can be used to construct a coarse-grained model for the heterogeneous bilayer that can predict the lateral organization and stability of rafts at up to millisecond time-scales.

The Conditions Favoring Between-Community Raiding in Chimpanzees, Bonobos, and Human Foragers.

Chimpanzees, bonobos, and human foragers share a fission-fusion social system and a mating system of joint male resource defense polygyny. Within-community skew in male strength varies among and within species. In this study, we extend a mathematical model of within-group male coalition formation among primates to derive the conditions for between-community conflicts in the form of raids. We show that the main factor affecting the presence of successful raiding is the likelihood of major discrepancies in party strength, which are set by party size distributions (and thus community size) and the skew in strength. This study confirms the functional similarities between the raiding of chimpanzees and human foragers, and it supports the "imbalance of power" hypothesis for raiding. However, it also proposes two amendments to this model. First, the absence of raiding in bonobos may be attributable more to potential female involvement in defense against raids, which increases the size of defensive coalitions. Second, the model attributes some of the raiding in humans to major contrasts in instantaneous fighting ability created by surprise raids on unarmed victims; it also draws attention to the distinction between minor raids and major raids that involve multiple bands of the same community.

DCMS: A data analytics and management system for molecular simulation.

Molecular Simulation (MS) is a powerful tool for studying physical/chemical features of large systems and has seen applications in many scientific and engineering domains. During the simulation process, the experiments generate a very large number of atoms and intend to observe their spatial and temporal relationships for scientific analysis. The sheer data volumes and their intensive interactions impose significant challenges for data accessing, managing, and analysis. To date, existing MS software systems fall short on storage and handling of MS data, mainly because of the missing of a platform to support applications that involve intensive data access and analytical process. In this paper, we present the database-centric molecular simulation (DCMS) system our team developed in the past few years. The main idea behind DCMS is to store MS data in a relational database management system (DBMS) to take advantage of the declarative query interface (i.e., SQL), data access methods, query processing, and optimization mechanisms of modern DBMSs. A unique challenge is to handle the analytical queries that are often compute-intensive. For that, we developed novel indexing and query processing strategies (including algorithms running on modern co-processors) as integrated components of the DBMS. As a result, researchers can upload and analyze their data using efficient functions implemented inside the DBMS. Index structures are generated to store analysis results that may be interesting to other users, so that the results are readily available without duplicating the analysis. We have developed a prototype of DCMS based on the PostgreSQL system and experiments using real MS data and workload show that DCMS significantly outperforms existing MS software systems. We also used it as a platform to test other data management issues such as security and compression.

Automated optimization of water-water interaction parameters for a coarse-grained model.

We have developed an automated parameter optimization software framework (ParOpt) that implements the Nelder-Mead simplex algorithm and applied it to a coarse-grained polarizable water model. The model employs a tabulated, modified Morse potential with decoupled short- and long-range interactions incorporating four water molecules per interaction site. Polarizability is introduced by the addition of a harmonic angle term defined among three charged points within each bead. The target function for parameter optimization was based on the experimental density, surface tension, electric field permittivity, and diffusion coefficient. The model was validated by comparison of statistical quantities with experimental observation. We found very good performance of the optimization procedure and good agreement of the model with experiment.

The effect of cross-link distributions in axially-ordered, cross-linked networks.

Cross-linking between the constituent chains of biopolymers has a marked effect on their materials' properties. In certain of these materials, such as fibrillar collagen, increases in cross-linking lead to an increase in the melting temperature. Fibrillar collagen is an axially-ordered network of cross-linked polymer chains exhibiting a broadened denaturation transition, which has been explained in terms of the successive denaturation with temperature of multiple species. We model axially-ordered, cross-linked materials as stiff chains with distinct arrangements of cross-link-forming sites. Simulations suggest that systems composed of chains with identical arrangements of cross-link-forming sites exhibit critical behavior. In contrast, systems composed of non-identical chains undergo a crossover. This model suggests that the arrangement of cross-link-forming sites may contribute to the broadening of the denaturation transition in fibrillar collagen.