Main page

Research Interests

Coarse Grained RNA models

Coming Soon.

Molecular Crowding of Nucleic Acid and Protein structures

Coming Soon.

Theory of Random and Block Copolymers

My interest in copolymers is related to an important technological issue of how to strengthen the interaction of molecules across an interface between two immiscible liquids. For example, a mixture of water and oil would generally be found in a phase separated state, with the two phases being almost pure water and oil. Traditionally nonionic surfactants which consist of a hydrophobic (oil soluble) tail and a hydrophilic (water soluble) head are used to stabilize water/oil mixtures. Molecules of such surfactants tend to localize at a water/oil interface and reduce its interfacial tension. In some cases the interfacial tension, defined as the free energy associated with the interface, can become negative, leading to the formation of water/oil microemulsions. In addition to nonionic surfactants, long copolymer chains that consist of both hydrophobic and hydrophilic monomers may be included in the mixture to facilitate mixing. Large molecular weight of polymer surfactants leads to a large energy gain vs. corresponding entropy loss when the surfactant molecules localize at an interface between two liquids. In particular, very long copolymer chains that on average show no preference for either water or oil will remain adsorbed at a water/oil interface over a broad range of temperatures. Diblock copolymers, which are generally used for compatibilizing purposes, can be quite laborious to synthesize. In my work (N. A. Denesyuk, 2001) I have used statistical physics methods to compare adsorption characteristics of diblock copolymers and of various random copolymers (which would be fairly easy to synthesize). I have shown that random copolymers can be effective compatibilizers, although their tendency to adsorb at interfaces is weaker than that of diblock copolymers. Another interesting result is that copolymers with randomly alternating sequences of hydrophobic and hydrophilic monomers adsorb at water/oil interfaces more easily than copolymers with regularly periodic sequences. This result is of entropic origin and is due to the fact that random copolymers, when localized at an interface, can adopt a greater number of low energy configurations than regular polyblock chains.

Another related problem is that of stabilizing homopolymer blends. In contrast to mixtures of low molecular weight liquids, a binary blend of homopolymers A and B will tend to segregate even at high temperatures. Such phase behavior is due to a large molecular weight of polymers: even if on the monomer level the incompatibility between chemical species A and B is fairly weak, it may add up to a very strong repulsion between complete polymer molecules. However many technological applications require mixing of immiscible homopolymers, which can be achieved by adding a "surfactant". Generally AB diblock copolymers are used as surfactants since their blocks A and B show affinity respectively for the A and B components of the homopolymer mixture. In particular, polymeric microemulsions akin to water/oil microemulsions have been successfully stabilized by a relatively low copolymer content in mixtures of weakly incompatible homopolymers.

Figure 1 A generic compositional phase diagram for a ternary blend of weakly incompatible homopolymers A and B and symmetric diblock copolymers AB. At high copolymer content, an A/B/AB polymer blend generally self-esembles into lamellar, hexagonal or BCC periodic structures. In addition, a small region of stability of the gyroid structures is found between the lamellar and hexagonal regions. The symmetry of an ordered structure depends on the relative amounts of homopolymers A and B. At low copolymer content, an A/B/AB polymer blend separates into A-rich and B-rich homogeneous phases. At intermediate (10-20%) copolymer content, polymeric microemulsions are formed, that represent macroscopically homogeneous phases with a spatially even distribution of homopolymers A and B.

Figure 2 This picture shows an interface between A-rich and B-rich microphase regions in a gyroid structure. Gyroid is a bicontinuous phase and, in this respect, it is somewhat similar to microemulsions. In contrast to microemulsions, however, gyroid is a spatially periodic structure belonging to the Ia3d space group.

Until very recently polymeric microemulsions could not be observed in mixtures of strongly incompatible homopolymers because the corresponding diblock copolymers AB tend to segregate in a separate phase from which homopolymers are excluded, rather than at the A/B interface. There has been some experimental evidence (J. Lee et al., 2003) that strongly immiscible homopolymers A and B can still be blended if the AB copolymer is substituted with its AC analogue, where the chemical interactions of the C-block with the A and B monomers is adjusted so as to provide a better compatibilizer. I have developed a numerical code which employs self-consistent field theory to study A/B/AB and A/B/AC polymeric alloys (N. A. Denesyuk and G. Gompper, 2006). I have shown that by appropriately adjusting the C block of the AC diblock copolymer, the concentration of surfactant required to stabilize an A/B polymer blend can be reduced by about 10%. In addition, I have shown that AB diblock copolymers may turn out to be much more efficient compatibilizers if their constituent blocks are larger than the homopolymers that need to be blended. In other words, the compatibilizing efficiency of diblock copolymers depends not only on their chemical nature but also on their degree of polymerization (which is not always recognized in experimental studies). These theoretical studies should facilitate the design of more effective polymer surfactants.

Theoretical Studies of Surface Wetting

If a small amount of liquid is placed on a solid (or liquid) substrate, in many cases it will adopt the shape of a droplet due to the liquid's surface tension. This situation, known as partial wetting, is described by a finite value of the contact angle between the liquid and the substrate surfaces. In the case of partial wetting, the substrate around the liquid droplet will be covered with a microscopically (a few angstroms) thin wetting film. However, if the liquid molecules show affinity for the underlying substrate, the liquid droplet may spread over the substrate despite its surface tension - the so-called complete wetting. Complete wetting can also be observed in the absence of a liquid droplet, for instance, when a substrate is placed in contact with a saturated vapor. In the case of strong attractive interactions between the vapor molecules and the substrate, a macroscopically thick film of dense liquid may condense onto the substrate surface. It is also possible for a liquid film to condense from an undersaturated vapor, even though liquid is not an equilibrium phase in the bulk. In the latter case, known as prewetting, the liquid film will adopt some mesoscopic thickness since the formation of a macroscopic non-equilibrium liquid phase is thermodynamically unfavorable. Apart from liquids on solid substrates, wetting is commonly studied in binary mixtures of liquids A and B which can have two, A-rich and B-rich, equilibrium phases. In such systems, a film of the B-rich phase may precipitate on the surface of the A-rich phase, even if B represents the heavier component.

Wetting transitions are the transitions that occur between the situations of partial and complete wetting (or prewetting). Wetting transitions are defined by a vanishing value of the contact angle or by a thin-to-thick transition in the thickness of the wetting film. Most experimentally observed wetting transitions are first order in which case, at some wetting temperature Tw, a microscopically thin wetting film changes discontinuously into a mesoscopic film (prewetting) or a macroscopic film (complete wetting). The complete wetting transitions can also be second order, i.e. continuous, resulting in a continuous divergence of the thickness of the wetting film as T®Tw (see review by D. Bonn and D. Ross, 2001). A wetting temperature Tw is generally found in the vicinity of the critical temperature Tc of the wetting liquid, where the liquid's surface tension is relatively small.

I have been interested in the influence of long-range Coulomb forces on the nature of wetting transitions. In particular, I have looked at the wetting transitions of ionic solutions near charged solid substrates, combining the classical Cahn-Landau theory for the solvent film with the nonliner Poisson-Boltzmann treatment of the salt solute (N. A. Denesyuk and J. P. Hansen, 2004). The result of this study has been the discovery of a large variety of wetting scenarios in salt solutions.

In very dilute salt solutions Coulomb interactions are mostly unscreened which leads to a discontinuous wetting transition, independently of the nature of the transition in the pure solvent. Electrostatic forces generally favor the formation of a macroscopically thick wetting film so that the wetting transition occurs further from the critical point than in the pure solvent. In concentrated ionic solutions electrostatic forces are largely screened and the wetting behavior is essentially the same as in the pure solvent. However, moderately concentrated ionic solutions show a very rich wetting behavior which includes i) a continuous or critical wetting transition, even if the transition for the pure solvent is first order, and ii) a sequence of two wetting transitions. In the latter case a discontinuous jump between two finite film thicknesses occurs first and is followed, upon further approach to the critical point, by a continuous divergence of the film thickness. This sequence of two wetting transitions is somewhat analogous to the prewetting transition from an undersaturated vapor followed by complete wetting on the liquid-vapor coexistence line (although in the present case both liquid and vapor are the equilibrium phases).

A similar sequence of wetting transitions has also been observed in systems dominated by long-range van der Waals forces. In these systems a negative sign of the Hamaker constant, describing the effective interaction between the two surfaces of the wetting film, corresponds to the net attraction between these surfaces and therefore disfavors complete wetting. When such a system is brought sufficiently close to its critical point it will undergo a discontinuous wetting transition into another partially wet state, characterized by a larger thickness of the wetting film. There are systems in which, upon further approach to the critical point, the Hamaker constant will in fact change sign from negative to positive, thus favoring complete wetting. In such systems the initial first-order transition will be followed by a continuous divergence of the film thickness at the point when the Hamaker constant changes its sign. Although apparently similar, the underlying physics of wetting by ionic solutions is quite different. In contrast with the van der Waals forces, the Coulomb forces always favor complete wetting so that a wetting transition will occur further away from the critical point. However, at those temperatures (or mixture compositions) complete wetting may still be prohibited by a large free energy cost associated with the liquid-vapor interface, so that instead of complete wetting the system will undergo a first-order transition to another partially wet state. Thus, if in the case of van der Waals interactions it is the short-range forces that favor complete wetting whereas the long-range forces disfavor it, in salt solutions the short-range and long-range forces switch these roles.

Generalized Debye Theory of Screening

This part of my research is concerned with the development and testing of a new method for the simulation of inhomogeneous molecular systems with long-range electrostatic interactions. With regard to soft matter and biological systems in particular, two methods have been commonly employed to account for electrostatic effects. The Debye theory of screening is usually invoked in analytical calculations, whose advantages are its simple mathematical form and clear physical meaning. However, Debye theory is a mean-field theory and, therefore, not suitable for the description of important correlations that may arise between charged objects at short distances. In addition, Debye theory is quantitatively accurate only in the case of extremely dilute ionic solutions with concentrations below typical physiological concentrations of 0.15M. In computer simulations, the Ewald image method is considered to be the most universal numerical method for treating electrostatic interactions. The disadvantage of the Ewald image method is that it approximates a large system with periodic images of a fundamental simulation box, and hence every inhomogeneity arising in the system is multiplied infinite number of times. This has in some cases produced artifacts in computer simulations of biological systems, such as false stabilization of more compact folded states of a protein.

I propose an alternative method to treat electrostatic interactions in simulation studies - generalized Debye theory of screening - which is based on Local Molecular Field theory and combines direct simulations of interactions between charges at short distances and the Debye theory of screening at large distances. The method contains an adjustable length parameter that is related to the minimal separation between charges beyond which the Debye approximation is expected to be quantitatively accurate; this length parameter can be optimized for the specific system under study. This new method has many advantages. First of all, the required computational time increases linearly with the system's size, which is particularly important for biological simulations in which typical system sizes are of order of 10000 - 100000 particles. Secondly, the proposed method does not rely on multiple periodic images of the simulation box and hence is free of any artifacts related to such periodicity. Thirdly, and probably most importantly, it departs from existing methods in that it can be applied to simulation boxes that are not strictly neutral but whose charge can fluctuate. This makes it particularly suited for studies of biological macromolecules with substantial net charge, such as DNA, that can strongly perturb the distribution of ions in the surrounding solution. In contrast, any method based on the infinite number of periodic images of the simulation box must impose the condition that the box is strictly neutral in order to ensure that the net charge density of the system of all box images is vanishing. In addition, generalized Debye theory of screening significantly broadens the class of inhomogeneous systems which are accessible to simulations. For instance, it is easily implemented for confined systems, for which the Ewald method is very involved and computationally demanding.

Simulation of a Model Protein

Recently I have been looking at the effects of electrostatics on the collapse of charged heterogeneous polymers. To this end I have developed a model of a heterogeneous polyelectrolyte, which may also be viewed as a very simple model of a protein. Although the folding of a protein molecule is known to be governed by the specific primary structure, there are still some advantages to considering simplified models. Such models generally keep simulation times short, thereby enabling the simulation of an ensemble of heterogeneous polyelectrolytes with different sequences (as opposed to the simulation of a single polyelectrolyte with a specific sequence). This in turn allows the study of electrostatic or other effects that are generic, rather than specific to one molecule.

I simulate polyelectrolyte chains which consist of hydrophobic and hydrophilic monomers, modeled as spherical beads. The hydrophobic species, which do not carry any electrostatic charge, interact via the attractive 12-6 Lennard-Jones potential to ensure the formation of a hydrophobic core at lower temperatures. The hydrophilic species are modeled by the purely repulsive 12-6 Lennard-Jones potential and they carry point-like charges near their surfaces. All polymer chains contain the same numbers of hydrophobic and hydrophilic beads, however these beads are distributed randomly to form a chain's sequence. The number of hydrophilic beads in a sequence is calculated based on the requirement that when a polymer chain is fully collapsed it should form a compact core comprised of all the hydrophobic beads, whereas all the hydrophilic ones should cover the polymer surface. This criterion yields 44 hydrophilic monomers for a chain whose total length is 100 monomers, which are the numbers that I use in my simulations. It is worth noting that the average content of hydrophilic aminoacids in globular proteins is 43.6%. I perform averages over an ensemble of sequences in order to identify those general properties that are sequence independent.

I would like to emphasize that my model is free, to a large extent, from a very common error arising in the simulations of polyelectrolytes. A well known effect of the counterion condensation takes place when the Coulomb interactions are strong and the counterions form ion pairs with the polymer charges. In this case, the polyelectrolyte chains collapse essentially as if they were neutral since there is no chain stretching due to electrostatic repulsions. However, when a polyelectrolyte chain is fully collapsed, the dielectric constant inside its core is approximately 2, as opposed to the dielectric constant of 80 for the water solvent. This fact is often neglected in simulations, where the electrostatic charges end up being buried inside the polymer interior, whereas their Coulomb interactions are assumed to be of the same magnitude as in the solvent. Since also, in most simulations, the buried counterions are of the same size as the polymer beads themselves, their excluded volume may strongly perturb the equilibrium polymer conformations. At least this part of the problem is circumvented in my model since, as a result of my sequence design, all the charged hydrophilic segments will tend to stay on the polymer surface. Of course some numerical errors remain since, when considering the Coulomb interactions of ions, I still neglect the presence of an interface between the regions of high and low dielectric constants. However I do not expect these numerical errors to be large because the surface of a polymer is not as dense as its core and this allows the solvent to envelop the hydrophilic beads as if they were almost entirely immersed in it.

Animation 1: In this animation a model protein consisting of 56 neutral hydrophobic monomers (gray spheres) and 44 charged hydrophilic monomers (green spheres) is simulated. As a first approximation, all 44 polymer charges (shown as red spheres) are assumed to have the same (nominal) negative sign and same magnitude. Positive and negative ions in solution are represented by blue and pink spheres. The density of ions is sustained by a grand canonical Monte Carlo, which causes ions to be added or removed from the simulation box based on the chemical potential. During the initial part of the animation only enough ions to neutralize the polymer are present and their number is kept fixed. In the second half of the animation the grand canonical routine is turned on and the density of ions equilibrates to a predetermined value. (Click image to view 14 Mb animated gif file)
Animation 2: This animation is similar to the one above, except that the density of ions is controlled by a grand canonical Monte Carlo during the entire animation. Halfway through the chemical potential is reduced by a factor of e2, causing the ions to re-equilibriate at a approximately half of the initial density. (Click image to view 9 Mb animated gif file)

To study the effects of electrostatics I place each polymer chain in a simulation box filled with a symmetric, low molecular weight, ionic solution. The diameter of the ions in this solution is taken to be significantly smaller than that of a polymer bead. Although both the polymer and low molecular weight ions (salt) are simulated by solving the Langevin dynamics equations of motion, it is not known a priori how much salt a simulation box of a given size may contain, given its inhomogeneous distribution around the polymer. I therefore control the salt concentration via the grand canonical MC insertions and deletions of ions (see animations above). I find that for typical physiological concentrations of ions (cM= 0.1-0.15M) the electrostatic forces are not efficiently screened and need to be taken into account to correctly predict the polymer's structure. In addition, truncation of these forces results in an incomplete screening of the polymer charge by the salt ions. In my simulations I treat electrostatics within the generalized Debye theory of screening which goes beyond the simple truncation of Coulomb interactions but, at the same time, is free from any artifacts which may result from taking into account multiple periodic images of charges (as in the Ewald sum method).

Dynamics of Helices

Another topic I am interested in is the formation of liquid crystalline phases of helical biopolymers. Virtually all biopolymers have helical structure: DNA, a-helices of globular proteins, tobacco mosaic virus (TMV), fd-viruses, actin, etc. In dense suspensions they form a variety of complex structures, the most well known among which are the nematic (achiral) and cholesteric (chiral) liquid crystalline phases. In addition, a liquid crystalline phase with a novel chiral symmetry, conical phase, has very recently been discovered in concentrated suspensions of helical flagella (E. Barry et al., 2006). However, very few theoretical tools are available to rationalize experimental data. In particular, the quantitative connection between microscopic chirality of individual molecules and macroscopic chirality of macromolecular aggregates is not yet understood. Another open question is why some helical biomolecules (TMV) form nematic phases, whereas others (fd-virus, DNA) form cholesteric ones.

The main problem that any theory would face here is the complexity of the helical structure. Mean-field theories, such as the self-consistent field theory of polymer melts commonly used to describe copolymer mesophases, neglect intermolecular correlations, which are crucial for the formation of chiral liquid crystalline phases. A chiral molecule that can rotate independently of other molecules will on average appear achiral, preventing formation of a chiral phase. Another common approximation that many theories rely on is the assumption of pairwise additive interactions between molecules. This assumption allows a system of many molecules to be described in terms of effective interactions between a pair of molecules, and thus reduces the task to the calculation of an effective interaction between two helices. Although justified in some cases, this is generally not true in charged systems where the long-range electrostatic forces lead to multimolecular correlations. The unfortunate circumstance is that liquid crystalline phases of biological helices are virtually always charged.

Animation 3: When an ensemble of 2000 attractive helices are simulated they soon begin to form a web of condensed quasi-crystalline formations. In their final configuration the helices are expected to organize into a liquid crystalline phase of defined spatial symmetry. However, observing the formation of such phases in simulations would require very long simulation times. The animation shows the view as the virtual camera moves through the simulation system in its intermediate, web-like configuration. (Click image to view 130 Mb animated gif file)
Animation 4: Animation of two helices which are subject to thermal fluctuations and a Lennard-Jones attraction. Dynamics of helices include important phenomena not captured when modeled as cylinders. These include the tendency of helices to approach each other at a non-zero twist angle, as well as rigid coupling of rotational and translational motion. (Click image to view 10 Mb animated gif file)

I am currently studying liquid crystalline phases of helical macromolecules in computer simulations. In contrast with analytical theory, computer simulations make it possible to treat the complex correlations between helical macromolecules. Furthermore, with currently available computational resources, it is possible to consider large ensembles of helices. Animation 3 shows an animated snapshot of a Langevin dynamics simulation of a system of 2000 model helical molecules. Each molecule is represented by a set of rigidly connected spherical beads whose centers of mass follow a helical path. The diameter and pitch of these helices have been chosen to correspond to those of the a-helices of globular proteins, and the bead-bead interactions have been modeled by the attractive Lennard-Jones potential. A close-up view of a pair of such helices in animation 4 indicates that they prefer to approach each other at some twist angle f not equal to 0. This non-vanishing value of the preferential twist angle is directly connected to the formation of macroscopically twisted (chiral) structures. The helix parameters in simulations can easily be adjusted to correspond to different biological systems. For instance, in contrast with the a-helices of proteins, the pitch of helical flagella is approximately an order of magnitude larger than their diameter. Consequently two helical flagella can lower their excluded volume by placing their helices in phase, which leads to the formation of a recently discovered conical phase (E. Barry et al., 2006). Apart from excluded volume and dispersion forces modeled by the Lennard-Jones potential, I also plan to include in my simulations electrostatic forces and explicit solvent (water) molecules.