Shankar Subramaniam

Analytical shape computation of macromolecules: I. molecular area and volume through alpha shape

The size and shape of macromolecules such as proteins and nucleic acids play an important role in their functions. Prior efforts to quantify these properties have been based on various discretization or tessellation procedures involving analytical or numerical computations. In this article, we present an analytically exact method for computing the metric properties of macromolecules based on the alpha shape theory. This method uses the duality between alpha complex and the weighted Voronoi decomposition of a molecule. We describe the intuitive ideas and concepts behind the alpha shape theory and the algorithm for computing areas and volumes of macromolecules. We apply our method to compute areas and volumes of a number of protein systems. We also discuss several difficulties commonly encountered in molecular shape computations and outline methods to overcome these problems. Proteins 33:1–17, 1998.

Analytical Shape Computation of Macromolecules: II. Inaccessible Cavities in Proteins

The structures of proteins are well‐packed, yet they contain numerous cavities which play key roles in accommodating small molecules, or enabling conformational changes. From high‐resolution structures it is possible to identify these cavities. We have developed a precise algorithm based on alpha shapes for measuring space‐filling‐based molecular models (such as van der Waals, solvent accessible, and molecular surface descriptions). We applied this method for accurate computation of the surface area and volume of cavities in several proteins. In addition, all of the atoms/residues lining the cavities are identified. We use this method to study the structure and the stability of proteins, as well as to locate cavities that could contain structural water molecules in the proton transport pathway in the membrane protein bacteriorhodopsin. Proteins 33:18–29, 1998.

Computation of molecular electrostatics with boundary element methods

In continuum approaches to molecular electrostatics, the boundary element method (BEM) can provide accurate solutions to the Poisson-Boltzmann equation. However, the numerical aspects of this method pose significant problems. We describe our approach, applying an alpha shape-based method to generate a high-quality mesh, which represents the shape and topology of the molecule precisely. We also describe an analytical method for mapping points from the planar mesh to their exact locations on the surface of the molecule. We demonstrate that derivative boundary integral formulation has numerical advantages over the nonderivative formulation: the well-conditioned influence matrix can be maintained without deterioration of the condition number when the number of the mesh elements scales up. Singular integrand kernels are characteristics of the BEM. Their accurate integration is an important issue. We describe variable transformations that allow accurate numerical integration. The latter is the only plausible integral evaluation method when using curve-shaped boundary elements.

Combined Monte Carlo and Molecular Dynamics Simulation of Fully Hydrated Dioleyl and Palmitoyl-oleyl Phosphatidylcholine Lipid Bilayers

We have applied a new equilibration procedure for the atomic level simulation of a hydrated lipid bilayer to hydrated bilayers of dioleyl-phosphatidylcholine (DOPC) and palmitoyl-oleyl phosphatidylcholine (POPC). The procedure consists of alternating molecular dynamics trajectory calculations in a constant surface tension and temperature ensemble with configurational bias Monte Carlo moves to different regions of the configuration space of the bilayer in a constant volume and temperature ensemble. The procedure is applied to bilayers of 128 molecules of POPC with 4628 water molecules, and 128 molecules of DOPC with 4825 water molecules. Progress toward equilibration is almost three times as fast in central processing unit (CPU) time compared with a purely molecular dynamics (MD) simulation. Equilibration is complete, as judged by the lack of energy drift in 200-ps runs of continuous MD. After the equilibrium state was reached, as determined by agreement between the simulation volume per lipid molecule with experiment, continuous MD was run in an ensemble in which the lateral area was restrained to fluctuate about a mean value and a pressure of 1 atm applied normal to the bilayer surface. Three separate continuous MD runs, 200 ps in duration each, separated by 10,000 CBMC steps, were carried out for each system. Properties of the systems were calculated and averaged over the three separate runs. Results of the simulations are presented and compared with experimental data and with other recent simulations of POPC and DOPC. Analysis of the hydration environment in the headgroups supports a mechanism by which unsaturation contributes to reduced transition temperatures. In this view, the relatively horizontal orientation of the unsaturated bond increases the area per lipid, resulting in increased water penetration between the headgroups. As a result the headgroup-headgroup interactions are attenuated and shielded, and this contributes to the lowered transition temperature.

Time-correlation analysis of simulated water motion in flexible and rigid gramicidin channels.

Molecular dynamics simulations have been done on a system consisting of the polypeptide membrane channel former gramicidin, plus water molecules in the channel and caps of waters at the two ends of the channel. In the absence of explicit simulation of the surrounding membrane, the helical form of the channel was maintained by artificial restraints on the peptide motion. The characteristic time constant of the artificial restraint was varied to assess the effect of the restraints on the channel structure and water motions. Time-correlation analysis was done on the motions of individual channel waters and on the motions of the center of mass of the channel waters. It is found that individual water molecules confined in the channel execute higher frequency motions than bulk water, for all degrees of channel peptide restraint. The center-of-mass motion of the chain of channel waters (which is the motion that is critical for transmembrane transport, due to the mandatory single filing of water in the channel) does not exhibit these higher frequency motions. The mobility of the water chain is dramatically reduced by holding the channel rigid. Thus permeation through the channel is not like flow through a rigid pipe; rather permeation is facilitated by peptide motion. For the looser restraints we used, the mobility of the water chain was not very much affected by the degree of restraint. Depending on which set of experiments is considered, the computed mobility of our water chain in the flexible channel is four to twenty times too high to account for the experimentally measured resistance of the gramicidin channel to water flow. From this result it appears likely that the peptide motions of an actual gramicidin channel embedded in a lipid membrane may be more restrained than in our flexible channel model, and that these restraints may be a significant modulator of channel permeability. For the completely rigid channel model the trapping of the water molecules in preferred positions throughout the molecular dynamics run precludes a reasonable assessment of mobility, but it seems to be quite low.

Water and polypeptide conformations in the gramicidin channel. A molecular dynamics study

Theoretical studies of ion channels address several important questions. The mechanism of ion transport, the role of water structure, the fluctuations of the protein channel itself, and the influence of structural changes are accessible from these studies. In this paper, we have carried out a 70-ps molecular dynamics simulation on a model structure of gramicidin A with channel waters. The backbone of the protein has been analyzed with respect to the orientation of the carbonyl and the amide groups. The results are in conformity with the experimental NMR data. The structure of water and the hydrogen bonding network are also investigated. It is found that the water molecules inside the channel act as a collective chain; whereas the conformation in which all the waters are oriented with the dipoles pointing along the axis of the channel is a preferred one, others are also accessed during the dynamics simulation. A collective coordinate involving the channel waters and some of the hydrogen bonding peptide partners is required to describe the transition of waters from one configuration to the other.

Explicit solvent models in protein pKa calculations.

Continuum methods for calculation of protein electrostatics treat buried and ordered water molecules by one of two approximations; either the dielectric constant of regions containing ordered water molecules is equal to the bulk solvent dielectric constant, or it is equal to the protein dielectric constant though no fixed atoms are used to represent water molecules. A method for calculating the titration behavior of individual residues in proteins has been tested on models of hen egg white lysozyme containing various numbers of explicit water molecules. Water molecules were included based on hydrogen bonding, solvent accessibility, and/or proximity to titrating groups in the protein. Inclusion of water molecules significantly alters the calculated titration behavior of individual titrating sites, shifting calculated pKa values by up to 0.5 pH unit. Our results suggest that approximately one water molecule within hydrogen-bonding distance of each charged group should be included in protein electrostatics calculations.

Computer modeling of electrostatic steering and orientational effects in antibody-antigen association

Brownian dynamics simulations are performed to investigate the role of long-range electrostatic forces in the association of the monoclonal antibody HyHEL-5 with hen egg lysozyme. The electrostatic field of the antibody is obtained from a solution of the nonlinear Poisson-Boltzmann using the x-ray crystal coordinates of this protein. The lysozyme is represented as an asymmetric dumbell consisting of two spheres of unequal size, an arrangement that allows for the modeling of the orientational requirements for docking. Calculations are done with the wild-type antibody and several point mutants at different ionic strengths. Changes in the charge distribution of the lysozyme are also considered. Results are compared with experiment and a simpler model in which the lysozyme is approximately by a single charged sphere.

Collective motion artifacts arising in long-duration molecular dynamics simulations

We tested a variety of molecular dynamics simulation strategies in long‐duration (up to several nanoseconds) constant‐temperature simulations of liquid water under periodic boundary conditions. Such long durations are necessary to achieve adequate conformational sampling in simulations of membrane assemblies and other large biomolecular systems. Under a variety of circumstances, serious artifacts arise in the form of spurious collective behavior that becomes obvious only after the simulation has gone at least several hundred picoseconds. The potential energy of the system drops and the system changes from a liquid to an icy or glassy state. The underlying cause is accumulated center‐of‐mass motion of the system, coupled with velocity rescaling associated with constant‐temperature control. The velocity rescaling in the constant‐temperature algorithm reduces the thermal velocity as the net center‐of‐mass velocity grows, effectively causing the kinetic energy of the system to drain from thermal motions into coordinated motions. We found that the incidence and magnitude of the underlying artifactual motion leading to the spurious transition is mediated by: choice of method for computing electrostatic interactions; choice of ensemble; size of the simulation cell; SHAKE tolerance; frequency of nonbonded pairlist updating; and closeness of coupling to the temperature bath. The appearance of the spurious transition can be avoided by periodically subtracting net center‐of‐mass motion during the dynamics, or by improving the accuracy of the simulation by means of tightening SHAKE tolerance and updating nonbonded pairlists every timestep.

APPLICATION OF COMBINED MONTE CARLO AND MOLECULAR DYNAMICS METHOD TO SIMULATION OF DIPALMITOYL PHOSPHATIDYLCHOLINE LIPID BILAYER

We describe a new equilibration procedure for the atomic level simulation of a hydrated lipid bilayer. The procedure consists of alternating molecular dynamics trajectory calculations in a constant surface tension and temperature ensemble with configurational bias Monte Carlo moves to different regions of the configuration space of the bilayer, in a constant volume and temperature ensemble. The procedure is described in detail and is applied to a bilayer of 100 molecules of dipalmitoyl phosphatidylcholine (DPPC) and 3205 water molecules. We find that the hybrid simulation procedure enhances the equilibration of the bilayer as measured by the convergence of the area per molecule and the segmental order parameters, as compared with a simulation using only molecular dynamics (MD). Progress toward equilibration is almost three times as fast in CPU time, compared with a purely MD simulation. Equilibration is complete, as judged by the lack of energy drift in three separate 200‐ps runs of continuous MD started from different initial states. Results of the simulation are presented and compared with experimental data and with other recent simulations of DPPC.u2003©1999 John Wiley & Sons, Inc.u2003J Comput Chem 20: 1153–1164, 1999