A quantitative quasispecies theory-based model of virus escape mutation under immune selection.

Viral infections involve a complex interplay of the immune response and escape mutation of the virus quasispecies inside a single host. Although fundamental aspects of such a balance of mutation and selection pressure have been established by the quasispecies theory decades ago, its implications have largely remained qualitative. Here, we present a quantitative approach to model the virus evolution under cytotoxic T-lymphocyte immune response. The virus quasispecies dynamics are explicitly represented by mutations in the combined sequence space of a set of epitopes within the viral genome. We stochastically simulated the growth of a viral population originating from a single wild-type founder virus and its recognition and clearance by the immune response, as well as the expansion of its genetic diversity. Applied to the immune escape of a simian immunodeficiency virus epitope, model predictions were quantitatively comparable to the experimental data. Within the model parameter space, we found two qualitatively different regimes of infectious disease pathogenesis, each representing alternative fates of the immune response: It can clear the infection in finite time or eventually be overwhelmed by viral growth and escape mutation. The latter regime exhibits the characteristic disease progression pattern of human immunodeficiency virus, while the former is bounded by maximum mutation rates that can be suppressed by the immune response. Our results demonstrate that, by explicitly representing epitope mutations and thus providing a genotype-phenotype map, the quasispecies theory can form the basis of a detailed sequence-specific model of real-world viral pathogens evolving under immune selection.

Thermodynamic basis for the emergence of genomes during prebiotic evolution.

The RNA world hypothesis views modern organisms as descendants of RNA molecules. The earliest RNA molecules must have been random sequences, from which the first genomes that coded for polymerase ribozymes emerged. The quasispecies theory by Eigen predicts the existence of an error threshold limiting genomic stability during such transitions, but does not address the spontaneity of changes. Following a recent theoretical approach, we applied the quasispecies theory combined with kinetic/thermodynamic descriptions of RNA replication to analyze the collective behavior of RNA replicators based on known experimental kinetics data. We find that, with increasing fidelity (relative rate of base-extension for Watson-Crick versus mismatched base pairs), replications without enzymes, with ribozymes, and with protein-based polymerases are above, near, and below a critical point, respectively. The prebiotic evolution therefore must have crossed this critical region. Over large regions of the phase diagram, fitness increases with increasing fidelity, biasing random drifts in sequence space toward 'crystallization.' This region encloses the experimental nonenzymatic fidelity value, favoring evolutions toward polymerase sequences with ever higher fidelity, despite error rates above the error catastrophe threshold. Our work shows that experimentally characterized kinetics and thermodynamics of RNA replication allow us to determine the physicochemical conditions required for the spontaneous crystallization of biological information. Our findings also suggest that among many potential oligomers capable of templated replication, RNAs may have evolved to form prebiotic genomes due to the value of their nonenzymatic fidelity.

Spontaneous buckling of lipid bilayer and vesicle budding induced by antimicrobial peptide magainin 2: a coarse-grained simulation study.

Molecular mechanisms of the action of antimicrobial peptides on bacterial membranes were studied by large scale coarse-grained simulations of magainin 2-dipalmitoylphosphatidylcholine/palmitoyloleoylphosphatidylglycerol (DPPC/POPG) mixed bilayer systems with spatial extents up to 0.1 µm containing up to 1600 peptides. Equilibrium simulations exhibit disordered toroidal pores stabilized by peptides. However, when a layer of peptides is placed near the lipid head groups on one side of the bilayer only, their incorporation leads to a spontaneous buckling of the bilayer. This buckling is followed by the formation of a quasi-spherical vesicular bud connected to the bilayer by a narrow neck. The mean curvature of the budding region is consistent with what is expected based on the dependence of the area per lipid on the peptide-to-lipid ratio in equilibrium simulations. Our simulations suggest that the incorporation of antimicrobial peptides on the exterior surface of a vesicle or a bacterial cell leads to buckling and vesicle budding, presumably accompanied by nucleations of giant transient pores of sizes that are much larger than indicated by equilibrium measurements and simulations.

Nonequilibrium phase transitions associated with DNA replication.

Thermodynamics governing the synthesis of DNA and RNA strands under a template is considered analytically and applied to the population dynamics of competing replicators. We find a nonequilibrium phase transition for high values of polymerase fidelity in a single replicator, where the two phases correspond to stationary states with higher elongation velocity and lower error rate than the other. At the critical point, the susceptibility linking velocity to thermodynamic force diverges. The overall behavior closely resembles the liquid-vapor phase transition in equilibrium. For a population of self-replicating macromolecules, Eigen's error catastrophe transition precedes this thermodynamic phase transition during starvation. For a given thermodynamic force, the fitness of replicators increases with increasing polymerase fidelity above a threshold.

ATP-induced conformational changes in Hsp70: molecular dynamics and experimental validation of an in silico predicted conformation.

The 70 kDa heat shock proteins (Hsp70s) play important roles in preventing the misfolding of proteins and repairing damage under stress by coupling ATP binding and hydrolysis to protein substrate release and binding, respectively. ATP binding is believed to induce closing of the Hsp70 nucleotide binding domain (NBD) around the nucleotide. We report here a combined computational-experimental study of this open-closed transition. All-atom molecular dynamics simulations were performed for isolated open state NBDs with and without bound ATP. The nucleotide-free NBD samples a wide range of open configurations exhibiting flexible rearrangements of its four subdomains (IA-IIB). In contrast, the ATP-bound Hsp70 NBD closes to a range of configurations that is substantially more closed than the conformation observed in crystals of ATP-complexed NBDs. The close approach of subdomains IB and IIB observed in the simulations results in a strong coordination of the fluorescence probe Trp90 of IB with Arg261 of IIB, a feature not seen in the crystal structures. To determine if this computationally observed conformation occurs in solution, we constructed an R261A mutant. The mutation was found to increase the K(m) and k(cat) for ATP and to significantly reduce the extent of the fluorescence quench observed upon ATP binding. Our results thus account for the previously unexplained ATP-driven change in Trp90 fluorescence seen in the isolated NBD.

Relaxation dynamics near nonequilibrium stationary states in Brownian ratchets.

A comprehensive study of the static and dynamical properties of a representative stochastic model of Brownian ratchet effects for molecular motors is reported. The model describes Brownian motions on two periodic potentials under static and time-dependent forces, where there are two distinct locations of chemical reactions coupling the levels with reversible rates within a period. Complete stationary properties have been obtained analytically for arbitrary potentials under external force. Dynamical relaxation properties near nonequilibrium stationary states were examined by considering the response function of velocity upon time-dependent external force, expressed in terms of the conditional probability density of the model. The latter is fully calculated using a systematic numerical method using matrix diagonalization, which is easily generalized to more complicated models for studying both static and dynamical properties. The behavior of the time-dependent response examined for model potentials suggests that the characteristic relaxation time near stationary states generally decreases linearly with respect to increasing velocity as one goes away from equilibrium via an increase in chemical potential of fuel species, a prediction testable in single molecule experiments.

Molecular dynamics studies of the energetics of translocation in model T7 RNA polymerase elongation complexes.

Translocation in the single subunit T7 RNA polymerase elongation complex was studied by molecular dynamics simulations using the posttranslocated crystal structure with the fingers domain open, an intermediate stable in the absence of pyrophosphate, magnesium ions, and nucleotide substrate. Unconstrained and umbrella sampling simulations were performed to examine the energetics of translocations. The extent of translocation was quantified using reaction coordinates representing the average and individual displacements of the RNA-DNA hybrid base pairs with respect to a reference structure. In addition, an unconstrained simulation was also performed for the product complex with the fingers domain closed, but with the pyrophosphate and magnesium removed, in order to examine the local stability of the pretranslocated closed state after the pyrophosphate release. The average spatial movement of the entire hybrid was found to be energetically costly in the post- to pretranslocated direction in the open state, while the pretranslocated state was stable in the closed complex, supporting the notion that the conformational state dictates the global stability of translocation states. However, spatial fluctuations of the RNA 3'-end in the open conformation were extensive, with the typical range reaching 3-4 A. Our results suggest that thermal fluctuations play more important roles in the translocation of individual nucleotides than in the movement of large sections of nucleotide strands: RNA 3'-end can move into and out of the active site within a single conformational state, while a global movement of the hybrid may be thermodynamically unfavorable without the conformational change.

Energetics of subdomain movements and fluorescence probe solvation environment change in ATP-bound myosin.

Energetics of conformational changes experienced by an ATP-bound myosin head detached from actin was studied by all-atom explicit water umbrella sampling simulations. The statistics of coupling between large scale domain movements and smaller scale structural features were examined, including the closing of the ATP binding pocket, and a number of key hydrogen bond formations shown to play roles in structural and biochemical studies. The statistics for the ATP binding pocket open/close transition show an evolution of the relative stability from the open state in the early stages of the recovery stroke to the stable closed state after the stroke. The change in solvation environment of the fluorescence probe Trp507 (scallop numbering; 501 in Dictyostelium discoideum) indicates that the probe faithfully reflects the closing of the binding pocket as previously shown in experimental studies, while being directly coupled to roughly the early half of the overall large scale conformational change of the converter domain rotation. The free energy change of this solvation environment change, in particular, is -1.3 kcal/mol, in close agreement with experimental estimates. In addition, our results provide direct molecular level data allowing for interpretations of the fluorescence experiments of myosin conformational change in terms of the de-solvation of Trp side chain.

Calculation of absolute protein-ligand binding constants with the molecular dynamics free energy perturbation method.

Reliable first-principles calculations of protein-ligand binding constants can play important roles in the study and characterization of biological recognition processes and applications to drug discovery. A detailed procedure for such a calculation is outlined in this chapter. The methodology is computationally implemented using the molecular dynamics sampling of relevant configurational spaces and free energy perturbation techniques. The procedure is illustrated with the model system of the phosphotyrosine peptide binding to the Src SH2 domain.

Exploration of the conformational space of myosin recovery stroke via molecular dynamics.

Muscle contractions are driven by cyclic conformational changes of myosin, whose molecular mechanisms of operation are being elucidated by recent advances in crystallographic studies and single molecule experiments. To complement such structural studies and consider the energetics of the conformational changes of myosin head, umbrella sampling molecular dynamics (MD) simulations were performed with the all-atom model of the scallop myosin sub-fragment 1 (S1) with a bound ATP in solution in explicit water using the crystallographic near-rigor and transition state conformations as two references. The constraints on RMSD reaction coordinates used for the umbrella sampling were found to steer the conformational changes efficiently, and relatively close correlations have been observed between the set of characteristic structural changes including the lever arm rotation and the closing of the nucleotide binding pocket. The lever arm angle and key residue interaction distances in the nucleotide binding pocket and the relay helix show gradual changes along the recovery stroke reaction coordinate, consistent with previous crystallographic and computational minimum energy studies. Thermal fluctuations, however, appear to make the switch-2 coordination of ATP more flexible than suggested by crystal structures. The local solvation environment of the fluorescence probe, Trp 507 (scallop numbering), also appears highly mobile in the presence of thermal fluctuations.

Exploration of molecular dynamics during transient sorption of fluids in mesoporous materials.

In recent years, considerable progress has been made in the development of novel porous materials with controlled architectures and pore sizes in the mesoporous range. An important feature of these materials is the phenomenon of adsorption hysteresis: for certain ranges of applied pressure, the amount of a molecular species adsorbed by the mesoporous host is higher on desorption than on adsorption, indicating a failure of the system to equilibrate. Although this phenomenon has been known for over a century, the underlying internal dynamics responsible for the hysteresis remain poorly understood. Here we present a combined experimental and theoretical study in which microscopic and macroscopic aspects of the relaxation dynamics associated with hysteresis are quantified by direct measurement and computer simulations of molecular models. Using nuclear magnetic resonance techniques and Vycor porous glass as a model mesoporous system, we have explored the relationship between molecular self-diffusion and global uptake dynamics. For states outside the hysteresis region, the relaxation process is found to be essentially diffusive in character; within the hysteresis region, the dynamics slow down dramatically and, at long times, are dominated by activated rearrangement of the adsorbate density within the host material.

Analytical theory of the nonequilibrium spatial distribution of RNA polymerase translocations.

A continuum Fokker-Planck model is considered for the RNA polymerase in the elongation phase, where the topology of a single free energy profile as a function of the translocation variable distinguishes the Brownian ratchet and power stroke mechanisms. The model yields a simple analytical stationary solution for arbitrary functional forms of the free energy. With the translocation potential of mean force estimated by the time-series data of the recent high-resolution single-molecule experiment [Abbondanzieri et al., Nature (London) 438, 460 (2005)], predictions of the model for the mechanical properties agree with experiments quantitatively with reasonable values of parameters. The evolution of the spatial distribution of translocation variable away from equilibrium with increasing nucleoside triphosphate concentration shows qualitatively different behavior in the two alternative scenarios, which could serve as an additional measurable signature of the underlying mechanism.

Protein-protein interactions in actin-myosin binding and structural effects of R405Q mutation: a molecular dynamics study.

Detailed residue-wise interactions involved in the binding of myosin to actin in the rigor conformation without nucleotides have been examined using molecular dynamics simulations of the chicken skeletal myosin head complexed with two actin monomers, based on the cryo-microscopic model of Holmes et al. (Nature 2003;425:423-427). The overall interaction is largely electrostatic in nature, because of the charged residues in the four loops surrounding the central primary binding site. The 50k/20k loop, disordered in crystal structures and in simulations of free myosin in solution, was found to be in a conformation stabilized with 1 - 2 internal salt bridges. The cardiomyopathy loop forms 2 - 3 interprotein salt bridges with actin monomers upon binding, whereas its Arg405 residue, the mutation site associated with the hypertrophic cardiomyopathy, forms a strong salt bridge with Glu605 in the neighboring helix away from actin in the actin-bound myosin. The myopathy loop of the R405Q mutant maintains a high degree of two-strand beta-sheet character when bound to actin with the corresponding salt bridges broken.

Calculation of absolute protein-ligand binding free energy from computer simulations.

A general methodology for calculating the equilibrium binding constant of a flexible ligand to a protein receptor is formulated on the basis of potentials of mean force. The overall process is decomposed into several stages that can be computed separately: the free ligand in the bulk is first restrained into the conformation it adopts in the bound state, position, and orientation by applying biasing potentials, then it is translated into the binding site, where it is released completely. The conformational restraining potential is based on the root-mean-square deviation of the peptide coordinates relative to its average conformation in the bound complex. Free energy contributions from each stage are calculated by means of free energy perturbation potential of mean force techniques by using appropriate order parameters. The present approach avoids the need to decouple the ligand from its surrounding (bulk solvent and receptor protein) as is traditionally performed in the double-decoupling scheme. It is believed that the present formulation will be particularly useful when the solvation free energy of the ligand is very large. As an application, the equilibrium binding constant of the phosphotyrosine peptide pYEEI to the Src homology 2 domain of human Lck has been calculated. The results are in good agreement with experimental values.

Grand canonical Monte Carlo simulations of water in protein environments.

The grand canonical simulation algorithm is considered as a general methodology to sample the configuration of water molecules confined within protein environments. First, the probability distribution of the number of water molecules and their configuration in a region of interest for biochemical simulations, such as the active site of a protein, is derived by considering a finite subvolume in open equilibrium with a large system serving as a bulk reservoir. It is shown that the influence of the bulk reservoir can be represented as a many-body potential of mean force acting on the atoms located inside the subvolume. The grand canonical Monte Carlo (GCMC) algorithm, augmented by a number of technical advances to increase the acceptance of insertion attempts, is implemented, and tested for simple systems. In particular, the method is illustrated in the case of a pure water box with periodic boundary conditions. In addition, finite spherical systems of pure water and containing a dialanine peptide, are simulated with GCMC while the influence of the surrounding infinite bulk is incorporated using the generalized solvent boundary potential [W. Im, S. Berneche, and B. Roux, J. Chem. Phys. 114, 2924 (2001)]. As a last illustration of water confined in the interior of a protein, the hydration of the central cavity of the KcsA potassium channel is simulated.

Modeling desorption of fluids from disordered mesoporous materials.

The desorption mechanism of fluids in disordered mesoporous glasses is studied by Monte Carlo simulations of a coarse-grained lattice model with realistic matrix configurations representative of Vycor. Two methods of simulation are considered: grand canonical ensemble Monte Carlo simulations and dynamic Monte Carlo simulations which mimic the diffusion of the fluid in and out of the material using Kawasaki dynamics. In the grand canonical simulations, cavitation via nucleation of bubbles inside the pores plays the dominant role in determining the fluid configurations along the desorption isotherm. The Kawasaki dynamics simulations indicate that such configurations are achieved dynamically via the gradual advancement of macroscopic front interfaces toward the interior. This is made possible by the bubble nucleation mechanism operating on a length scale that is determined by both the typical pore size and the strength of the solid-fluid interaction.

Phase behavior and dynamics of fluids in mesoporous glasses.

Equilibrium and dynamical relaxation behavior of fluids confined in disordered mesoporous glasses such as Vycor are studied based on a lattice model using mean field theory and Monte Carlo simulations. Preferential attractive interactions between the solid surfaces and the fluid suppresses macroscopic phase separation, while making the relaxation rate increasingly slow. The free energy landscape characterized by the presence of the many metastable minima separated by finite barriers dominates both the static and dynamic behavior of fluids at low temperature. Our results provide additional insight into the nature of hysteresis in adsorption measurements of gases in porous glasses.