Temperature-and-density-scaling Monte Carlo: methodology and the canonical thermodynamics of Lennard-Jonesium

Umbrella sampling was used to study, within a Monte-Carlo run, substantial ranges of both temperature and density: i.e. “temperature-and-density-scaling Monte-Carlo” (TDSMC). The paper reviews in detail some of the theory and the practicalities of using the TDSMC technique. For example, we discuss how to generate an appropriate sampling distribution, and the question of appropriate error analysis for TDSMC data. In order to test whether the sampling efficiency might be improved by “stratifying”, the entire investigation was carried out in two ways: by covering the target region of thermodynamic space in a single TDSMC run, and by independently segmenting the same target region into three overlapping subregions, investigated by separate TDSMC runs: the results and relative efficiencies are compared.

These matters, and the potential of the method, are illustrated by application to a Lennard-Jones fluid. The target region covered a substantial region of the thermodynamic space, including the upper part of the gas–liquid coexistence curve, and the TDSMC run provides direct and accurate estimates of relative free energy throughout the region, along with other properties such as the pressure, internal energy, chemical potential and heat capacity. The results are also used to characterise the thermodynamics of the liquid–gas transition in the canonical ensemble. The precision of all these results is substantially better than in other simulation methods of studying such properties.     

Importance of coulombic end effects on cation accumulation near oligoelectrolyte B-DNA: a demonstration using 23Na NMR.

The local cation concentration at the surface of oligomeric or polymeric B-DNA is expected, on the basis of MC simulations (Olmsted, M. C., C. F. Anderson, and M. T. Record, Jr. 1989. Proc. Natl. Acad. Sci. USA. 86:7766-7770), to decrease sharply as either end of the molecule is approached. In this paper we report 23Na NMR measurements indicating the importance of this "coulombic" end effect on the average extent of association of Na+ with oligomeric duplex DNA. In solutions containing either 20-bp synthetic DNA or 160-bp mononucleosomal calf thymus DNA at phosphate monomer concentrations [P] of 4-10 mM, measurements were made over the range of ratios 1 < or = [Na]/[LP] < or = 20, corresponding to Na+ concentrations of 4-200 nM. The longitudinal 23Na NMR relaxation rates measured in these NaDNA solutions, Robs, are interpreted as population-weighted averages of contributions from "bound" (RB) and "free" (RF) 23Na relaxation rates. The observed enhancements of Robs indicate that RB significantly exceeds RF, which is approximately equal to the 23Na relaxation rate in an aqueous solution containing only NaCl. Under salt-fre-tconditions ([Na]/[P] = 1), where the enhancement in Robs is maximal, we find that Robs--RF in the solution containing 160-bp DNA is approximately 1.8 times that observed for the 20-bp DNA. For the 160-bp oligomer (which theoretical calculations predict to be effectively polyion-like), we find that a plot of Robs v. [P]/[Na] is linear, as observed previously for sonicated (approximately 700 bp) DNA samples. For the 20-bp oligonucleotide this plot exhibits a marked departure from linearity that can be fitted to a quadratic function of [P]/[Na]. Monte Carlo simulations based on a simplified model are capable of reproducing the qualitative trends in the 23Na NMR measurements analyzed here. In particular, the dependences of Robs--RF on DNA charge magnitude of Z(320 vs. 38 phosphates) and (for the 20-bp oligomer) on [Na]/[P] are well correlated with the calculated average surface concentration of Na+. Thus, effects of sodium concentration on RB appear to be of secondary importance. We conclude that 23Na NMR relaxation measurements are a sensitive probe of the effects of oligomer charge on the extent of ion accumulation near B-DNA oligonucleotides, as a function of [Na] and [P].     

Hydrogen adsorption capacities of multi-walled boron nitride nanotubes and nanotube arrays: a grand canonical Monte Carlo study

Hydrogen adsorption in multi-walled boron nitride nanotubes and their arrays was studied using grand canonical Monte Carlo simulation. The results show that hydrogen storage increases with tube diameter and the distance between the tubes in multi-walled boron nitride nanotube arrays. Also, triple-walled boron nitride nanotubes present the lowest level of hydrogen physisorption, double-walled boron nitride nanotubes adsorb hydrogen better when the diameter of the inner tube diameter is sufficiently large, and single-walled boron nitride nanotubes adsorb hydrogen well when the tube diameter is small enough. Boron nitride nanotube arrays adsorb hydrogen, but the percentage of adsorbed hydrogen (by weight) in boron nitride nanotube arrays is rather similar to that found in multi-walled boron nitride nanotubes. Also, when the Langmuir and Langmuir-Freundlich equations were fitted to the simulated data, it was found that multi-layer adsorptivity occurs more prominently as the number of walls and the tube diameter increase. However, in single-walled boron nitride nanotubes with a small diameter, the dominant mechanism is monolayer adsorptivity.     

Confinement effects on the thermodynamics of protein folding: Monte Carlo simulations

The effects of chaperonin-like cage-induced confinement on protein stability have been studied for molecules of varying sizes and topologies. Minimalist models based on Gō-like interactions are employed for the proteins, and density-of-states-based Monte Carlo simulations are performed to accurately characterize the thermodynamic transitions. This method permits efficient sampling of conformational space and yields precise estimates of free energy and entropic changes associated with protein folding. We find that confinement-driven stabilization is not only dependent on protein size and cage radius, but also on the specific topology. The choice of the confining potential is also shown to have an effect on the observed stabilization and the scaling behavior of the stabilization with respect to the cage size.     

Monte Carlo simulation of hydration of the nucleic acid fragments

Systems containing a base or a base pair and 25 water molecules, as well as a helical stack and 30 water molecules per base pair, have been simulated. Changes in the base hydration shell structure, after the bases have been included into the pair and then into the base pair stack, are discussed. Hydration shells of several configurations of the base pair stacks are discussed. Probabilities of formation of the hydrogen-bonded bridges of 1, 2 and 3 water molecules between hydrophilic centres have been estimated. The hydration shell structure was shown to depend on the nature of the base pair and on the stack configuration, while dependence of the global hydration shell characteristics on the stack configuration has been proved to be rather slight. The most typical structural elements of hydration shells, in the glycosidic (minor in B-like conformation) and non-glycosidic (major) grooves, for different configurations of AU and GC stacks, have been found and discussed. The number of hydrogen bonds between water molecules and bases per water molecule was shown to change upon transformation of the stack from A to B configuration. This result is discussed in connection with the reasons for B to A conformational transition and the concept of "water economy". Hydration shell patterns of NH2-groups of AU and GC helical stacks differ significantly.     

A novel algorithm to accelerate the convergence of grand canonical Monte Carlo simulation of non-uniform fluids

A novel Monte Carlo scheme, MTZ-GCMC, utilising the ‘mass transfer zone’ (MTZ) between the fluid phase and the adsorbed phase, is proposed as an effective method for simulating the approach of a non-uniform fluid to equilibrium. We have applied this procedure to study the adsorption of gases in a closed-end graphitic slit pore, paying special attention to the region where there is a very sharp (but not vertical) condensation transition within the pore. In this region, conventional Monte Carlo (MC) requires a much greater computational effort to achieve convergence and may lead to incorrect results. In the MTZ-GCMC scheme, the insertion/deletion trials during the course of a simulation are restricted to the MTZ, which is identified and characterised by the percentage of successful insertions (PSI). Since the MTZ changes during the equilibration stage, the PSI must be updated on the fly until equilibrium has been reached. The distinct advantage of this scheme is that unnecessary insertion and deletion trials in the dense adsorbed phase, and the gas-like region are avoided; since most of these trials would be rejected, this means that convergence to the true equilibrium can be achieved an order of magnitude faster.     

Multivariate curve resolution: a powerful tool for the analysis of conformational transitions in nucleic acids

A successful application is reported of the multivariate curve resolution alternating least-squares method (MCR-ALS) for the analysis of nucleic acid melting and salt-induced transitions. Under conditions where several structures co-exist in a conformational equilibrium, MCR-ALS analysis of the UV and circular dichroism (CD) spectra at different temperatures, ionic strength and oligonucleotide concentration allows for the resolution of concentration profiles and pure spectra of the different species. The methodology is illustrated by the case of the cyclic oligonucleotide d<pTGCTCGCT>. The melting transition of this molecule at different oligonucleotide concentrations was studied at 0, 2 and 10 mM MgCl₂ by UV and CD spectroscopy. In addition, salt titration experiments were carried out at 21.0 and 54.0°C. The MCR-ALS analysis indicates that three different conformations of this molecule co-exist in solution. In agreement with previous NMR studies, these conformations were assigned to a monomeric dumbbell-like structure, a dimeric four-stranded conformation and a disordered (random coil) structure. The MCR-ALS methodology allows for a detailed analysis of how this equilibrium is affected by temperature, salt and oligonucleotide concentration.     

Nature of the stacking of nucleic acid bases in water: a Monte Carlo simulation

The results of a Monte Carlo simulation of the hydration of uracil and thymine molecules, their stacked dimers and hydrogen-bonded base pairs are presented. Simulations have been performed in a cluster approximation. The semiempirical atom-atom potential functions have been used (cluster consisting of 200 water molecules). It has been shown that the stacking interactions of uracil and thymine molecules in water arise mainly due to the increase in the water-water interaction during the transition from monomers to dimer. It has been found out that stacked base associates are more preferable than base pairs in water. This preference is mainly due to the energetically more favourable structure of water around the stack.     

Solvation studies of DMP323 and A76928 bound to HIV protease: analysis of water sites using grand canonical Monte Carlo simulations.

We examine the water solvation of the complex of the inhibitors DMP323 and A76928 bound to HIV-1 protease through grand canonical Monte Carlo simulations, and demonstrate the ability of this method to reproduce crystal waters and effectively predict water positions not seen in the DMP323 or A76928 structures. The simulation method is useful for identifying structurally important waters that may not be resolved in the crystal structures. It can also be used to identify water positions around a putative drug candidate docked into a binding pocket. Knowledge of these water positions may be useful in designing drugs to utilize them as bridging groups or displace them in the binding pocket. In addition, the method should be useful in finding water sites in homology models of enzymes for which crystal structures are unavailable.     

Methane storage in homogeneous armchair open-ended single-walled boron nitride nanotube triangular arrays: a grand canonical Monte Carlo simulation study

The physisorption of methane in homogeneous armchair open-ended SWBNNT triangular arrays was evaluated using grand canonical ensemble Monte Carlo simulation for tubes 11.08, 13.85, 16.62, and 19.41 ? [(8,8), (10,10), (12,12), and (14,14), respectively] in diameter, at temperatures of 273, 298, 323, and 373 K, and at fugacities of 0.5-9.0 Mpa. The intermolecular forces were modeled using the Lennard-Jones potential model. The absolute, excess, and delivery adsorption isotherms of methane were calculated for the various boron nitride nanotube arrays. The specific surface areas and the isosteric heats of adsorption, Q(st), were also studied, different isotherm models were fitted to the simulated adsorption data, and the model parameters were correlated. According to the results, it is possible to reach 108% and 140% of the US Department of Energy's target for CH(4) storage (180 v/v at 298 K and 35 bar) using the SWBNNT array with nanotubes 16.62 and 19.41 ? in diameter, respectively, as adsorbent. The results show that for a van der Waals gap of 3.4 ?, there is no interstitial adsorption except for arrays containing nanotubes with diameters of >15.8 ?. Multilayer adsorption starts to occur in arrays containing nanotubes with diameters of >16.62 ?, and the minimum pressure required for multilayer adsorption is 1.0 MPa. A brief comparison of the methane adsorption capacities of single-walled carbon and boron nitride nanotube arrays was also performed.     

Aqueous hydration of nucleic acid constituents: Monte Carlo computer simulation studies

Monte Carlo computer simulations were performed on dilute aqueous solutions of thymine, cytosine, uracil, adenine, guanine, the dimethyl phosphate anion in the gauche-gauche conformation and a ribose and deoxyribose derivative. The aqueous hydration of each molecule was analysed in terms of quasi-component distribution functions based on the Proximity Criterion, and partitioned into hydrophobic, hydrophilic and ionic contributions. Color stereo views of selected hydration complexes are also presented. A preliminary discussion of the transferability of functional group coordination numbers is given. The results enable to comment on two current problems related to the hydration of nucleic acids: a) the theory of Dickerson and coworkers on the role of water in the relative stability of the A and B form of DNA and b) the idea of water bridges and filaments emerging from the computer simulation results on the hydration of DNA fragments by Clementi.     

Development of a grand canonical-kinetic Monte Carlo scheme for simulation of mixtures

A rejection-free methodology-based kinetic Monte Carlo (kMC) method has been developed in the grand canonical ensemble to simulate fluid mixtures. It comprises two different moves: entropic displacement of a selected molecule (based on the Rosenbluth algorithm) in the volume space of the system, and exchange of molecules with the surroundings (insertion or deletion). These two moves are made sequentially with M displacement moves followed by one exchange. The displacement moves are treated as sub-NVT sequences within a grand canonical ensemble. The procedure for deletion or insertion of a molecule is either, based on the Rosenbluth algorithm, or on a direct comparison, in which the average activity of one component is compared with its specified activity. The components are chosen either with equal probability or with a probability proportional to their density. The implementation of rejection-free kMC is much simpler than the Metropolis importance sampling MC procedure, which requires three different types of move, all of which must be tested for acceptance or rejection. The new scheme has been evaluated by applying it to fluid argon and to an equimolar mixture of methane, ethane and propane.     

Spermine: an "invisible" component in the crystals of B-DNA. A grand canonical Monte Carlo and molecular dynamics simulation study

The association of spermine(4+) (Spm(4+)), Mg(2+) and monovalent (M(+)) ions with DNA in crystal form, have been studied using grand canonical Monte Carlo (GCMC) and molecular dynamics (MD) computer simulations. GCMC calculations were used to calculate the distribution of Spm(4+), Mg(2+), and M(+) between the equilibrating solvent and the DNA crystal under conditions mimicking the crystal-growing protocols reported in a number of recent X-ray diffraction studies of DNA oligomers. The GCMC simulations show that the composition of ions neutralizing the negative charge of DNA can vary in a broad range. The GCMC simulations were used to provide appropriate conditions for subsequent 6 ns constant pressure and temperature MD simulations of DNA in a typical crystalline environment consisting of three DNA double helix decamers in a periodic hexagonal cell, containing 1200 water molecules, eight Spm(4+), 32 Na(+) and four Cl(-) ions. Based on the simulation results, it seems possible to give an explanation why spermine molecules are usually not detected in X-ray studies in spite of their high concentration in the preparatory samples used as the crystallizing agent. It appears that this flexible polyamine molecule has several binding modes, interacting in fairly irregular manner with different sites on DNA and showing no regular ordering in the DNA crystals. Ions of Na(+) and Spm(4+) compete with each other and with water molecules in binding to bases in the minor groove and they influence the structure of the DNA hydration shell in different ways.     

Grand canonical Monte Carlo molecular and thermodynamic predictions of ion effects on binding of an oligocation (L8+) to the center of DNA oligomers.

Grand canonical Monte Carlo (GCMC) simulations are reported for aqueous solutions containing excess univalent salt (activities a +/- = 1.76-12.3 mM) and one of the following species: an octacationic rod-like ligand, L8+; a B-DNA oligomer with N phosphate charges (8 < or = N < or = 100); or a complex resulting from the binding of L8+ at the center of an N-mer (24 < or = N < or = 250). Simplified models of these multiply charged species are used in the GCMC simulations to predict the fundamental coulombic contributions to the following experimentally relevant properties: 1) the axial distance over which ligand binding affects local counterion concentrations at the surface of the N-mer; 2) the dependence on N of GCMC preferential interaction coefficients, gamma 32MC identical to delta C3/delta C2l a +/-, T, where C3 and C2 are, respectively, the molar concentrations of salt and the multiply charged species (ligand, N-mer or complex); and 3) the dependence on N of SaKobs identical to d in Kobs/d in a +/- = delta (magnitude of ZJ + 2 gamma 32J), where Kobs is the equilibrium concentration quotient for the binding of L8+ to the center of an N-mer and delta denotes the stoichiometric combination of terms, each of which pertains to a reactant or product J having magnitude of ZJ charges. The participation of electrolyte ions in the ligand binding interaction is quantified by the magnitude of SaKobs, which reflects the net (stoichiometrically weighted) difference in the extent of thermodynamic binding of salt ions to the products and reactants. Results obtained here from GCMC simulations yield a picture of the salient molecular consequences of binding a cationic ligand, as well as thermodynamic predictions whose applicability can be tested experimentally. Formation of the central complex is predicted to cause a dramatic reduction in the surface counterion (e.g., Na+) concentration over a region including but extending well beyond the location of the ligand binding site. For binding a cationic ligand, SaKobs is predicted to be negative, indicating net electrolyte ion release in the binding process. At small enough N, -SaKobs is predicted to decrease strongly toward zero with decreasing N. At intermediate N, -SaKobs appears to exceed its limiting value as N-->infinity.     

The electrostatically driven Monte Carlo method: application to conformational analysis of decaglycine

The Electrostatically Driven Monte Carlo (EDMC) method was applied in a study of a decamer of glycine whose conformational behavior is described by the Empirical Conformational Energy Program for Peptides (ECEPP/2) potential energy model. When free neutral end groups were used, it was found that conformations that were not alpha-helical had significantly lower potential energies than fully alpha-helical ones. However, when the N- and C-termini were blocked by acetyl and methyl amide groups, respectively, the number of unsatisfied hydrogen-bond donors and acceptors at the helix termini was diminished from 8 to 6; in this case, the possibility of forming two additional alpha-helical hydrogen bonds was an important enough factor in making the alpha-helical conformation the one with the lowest energy. The EDMC method was used as a global energy optimizer since it does not often become trapped in high-energy local minima.     

Modeling of hydration of incorrect nucleic acid base pairs by the Monte Carlo method

The hydration of water bridged base pairs of nucleic acids have been simulated via the Monte Carlo method. The simulation have shown that water molecules forming H-bonds with both bases preserve this H-bonding with large probability in the water surrounding. This fact supports the supposition about the important role of water molecules in wrong base pair formation and about the role of these base pairs in the structure and functioning of nucleic acids.     

3D-Lattice Monte Carlo simulations of model proteins. Size effects on folding thermodynamics and kinetics

Recently, we devised an energy scale to vary systematically amino-acid residue-solvent interactions for Monte Carlo simulations of lattice-model proteins in water. For 27-mer proteins, the folding behavior varies appreciably with the choice of interaction parameters. We now perform similar simulations with 64-mers to study the size dependence of the optimal energy parameter set for representing realistic behavior typical of many real proteins (i.e. fast folding and high cooperativity for single chains). We find that 64-mers are considerably more stable and more cooperative compared to 27-mers. The optimal interfacial-interaction-energy parameter set, however, is relatively size independent.     

Ligand-induced conformational changes in tissue transglutaminase: Monte Carlo analysis of small-angle scattering data

Small-angle neutron and x-ray scattering experiments have been performed on type 2 tissular transglutaminase to characterize the conformational changes that bring about Ca(2+) activation and guanosine triphosphate (GTP) inhibition. The native and a proteolyzed form of the enzyme, in the presence and in the absence of the two effectors, were considered. To describe the shape of transglutaminase in the different conformations, a Monte Carlo method for calculating small-angle neutron scattering profiles was developed by taking into account the computer-designed structure of the native transglutaminase, the results of the Guinier analysis, and the essential role played by the solvent-exposed peptide loop for the conformational changes of the protein after activation. Although the range of the neutron scattering data is rather limited, by using the Monte Carlo analysis, and because the structure of the native protein is available, the distribution of the protein conformations after ligand interaction was obtained. Calcium activation promotes a rotation of the C-terminal with respect to the N-terminal domain around the solvent-exposed peptide loop that connects the two regions. The psi angle between the longest axes of the two pairs of domains is found to be above 50 degrees, larger than the psi value of 35 degrees calculated for the native transglutaminase. On the other hand, the addition of GTP makes possible conformations characterized by psi angles lower than 34 degrees. These results are in good agreement with the proposed enzyme activity regulation: in the presence of GTP, the catalytic site is shielded by the more compact protein structure, while the conformational changes induced by Ca(2+) make the active site accessible to the substrate.     

Ab initio conformational analysis of nucleic acid components: intrinsic energetic contributions to nucleic acid structure and dynamics

In recent years, the use of high-level ab initio calculations has allowed for the intrinsic conformational properties of nucleic acid building blocks to be revisited. This has provided new insights into the intrinsic conformational energetics of these compounds and its relationship to nucleic acids structure and dynamics. In this article we review recent developments and present new results. New data include comparison of various levels of theory on conformational properties of nucleic acid building blocks, calculations on the abasic sugar, known to occur in vivo in DNA, on the TA conformation of DNA observed in the complex with the TATA box binding protein, and on inosine. Tests of the Hartree-Fock (HF), second-order M?ller-Plesset (MP2), and Density Functional Theory/Becke3, Lee, Yang and Par (DFT/B3LYP) levels of theory show the overall shape of backbone torsional energy profiles (for gamma, epsilon, and chi) to be similar for the different levels, though some systematic differences are identified between the MP2 and DFT/B3LYP profiles. The east pseudorotation energy barrier in deoxyribonucleosides is also sensitive to the level of theory, with the HF and DFT/B3LYP east barriers being significantly lower (approximately 2.5 kcal/mol) than the MP2 counterpart (approximately 4.0 kcal/mol). Additional calculations at various levels of theory suggest that the east barrier in deoxyribonucleosides is between 3.0 and 4.0 kcal/mol. In the abasic sugar, the west pseudorotation energy barrier is found to be slightly lower than the east barrier and the south pucker is favored more than in standard nucleosides. Results on the TA conformation suggest that, at the nucleoside level, this conformation is significantly destabilized relative to the global energy minimum, or relative to the A- and B-DNA conformations. Deoxyribocytosine would destabilize the TA conformation more than other bases relative to the A-DNA conformation, but not relative to the B-DNA conformation.     

Grand canonical ensemble Monte Carlo simulation of the dCpG/proflavine crystal hydrate.

The grand canonical ensemble Monte Carlo molecular simulation method is used to investigate hydration patterns in the crystal hydrate structure of the dCpG/proflavine intercalated complex. The objective of this study is to show by example that the recently advocated grand canonical ensemble simulation is a computationally efficient method for determining the positions of the hydrating water molecules in protein and nucleic acid structures. A detailed molecular simulation convergence analysis and an analogous comparison of the theoretical results with experiments clearly show that the grand ensemble simulations can be far more advantageous than the comparable canonical ensemble simulations.