Dynamics and mechanisms of CRISPR-Cas9 through the lens of computational methods

The clustered regularly interspaced short palindromic repeat (CRISPR) genome-editing revolution established the beginning of a new era in life sciences. Here, we review the role of state-of-the-art computations in the CRISPR-Cas9 revolution, from the early refinement of cryo-EM data to enhanced simulations of large-scale conformational transitions. Molecular simulations reported a mechanism for RNA binding and the formation of a catalytically competent Cas9 enzyme, in agreement with subsequent structural studies. Inspired by single-molecule experiments, molecular dynamics offered a rationale for the onset of off-target effects, while graph theory unveiled the allosteric regulation. Finally, the use of a mixed quantum-classical approach established the catalytic mechanism of DNA cleavage. Overall, molecular simulations have been instrumental in understanding the dynamics and mechanism of CRISPR-Cas9, contributing to understanding function, catalysis, allostery, and specificity.     

Speeding up Monte Carlo molecular simulation by a non-conservative early rejection scheme

Monte Carlo (MC) molecular simulation describes fluid systems with rich information, and it is capable of predicting many fluid properties of engineering interest. In general, it is more accurate and representative than equations of state. On the other hand, it requires much more computational effort and simulation time. For that purpose, several techniques have been developed in order to speed up MC molecular simulations while preserving their precision. In particular, early rejection schemes are capable of reducing computational cost by reaching the rejection decision for the undesired MC trials at an earlier stage in comparison to the conventional scheme. In a recent work, we have introduced a ‘conservative’ early rejection scheme as a method to accelerate MC simulations while producing exactly the same results as the conventional algorithm. In this paper, we introduce a ‘non-conservative’ early rejection scheme, which is much faster than the conservative scheme, yet it preserves the precision of the method. The proposed scheme is tested for systems of structureless Lennard-Jones particles in both canonical and NVT-Gibbs ensembles. Numerical experiments were conducted at several thermodynamic conditions for different number of particles. Results show that at certain thermodynamic conditions, the non-conservative method is capable of doubling the speed of the MC molecular simulations in both canonical and NVT-Gibbs ensembles.     

Effective energy function for proteins in solution

A Gaussian solvent-exclusion model for the solvation free energy is developed. It is based on theoretical considerations and parametrized with experimental data. When combined with the CHARMM 19 polar hydrogen energy function, it provides an effective energy function (EEF1) for proteins in solution. The solvation model assumes that the solvation free energy of a protein molecule is a sum of group contributions, which are determined from values for small model compounds. For charged groups, the self-energy contribution is accounted for primarily by the exclusion model. Ionic side-chains are neutralized, and a distance-dependent dielectric constant is used to approximate the charge-charge interactions in solution. The resulting EEF1 is subjected to a number of tests. Molecular dynamics simulations at room temperature of several proteins in their native conformation are performed, and stable trajectories are obtained. The deviations from the experimental structures are similar to those observed in explicit water simulations. The calculated enthalpy of unfolding of a polyalanine helix is found to be in good agreement with experimental data. Results reported elsewhere show that EEF1 clearly distinguishes correctly from incorrectly folded proteins, both in static energy evaluations and in molecular dynamics simulations and that unfolding pathways obtained by high-temperature molecular dynamics simulations agree with those obtained by explicit water simulations. Thus, this energy function appears to provide a realistic first approximation to the effective energy hypersurface of proteins.     

Estimating chikungunya virus transmission parameters and vector control effectiveness highlights key factors to mitigate arboviral disease outbreaks

BackgroundViruses transmitted by Aedes mosquitoes have greatly expanded their geographic range in recent decades. They are considered emerging public health threats throughout the world, including Europe. Therefore, public health authorities must be prepared by quantifying the potential magnitude of virus transmission and the effectiveness of interventions.MethodologyWe developed a mathematical model with a vector-host structure for chikungunya virus transmission and estimated model parameters from epidemiological data of the two main autochthonous chikungunya virus transmission events that occurred in Southern France, in Montpellier (2014) and in Le Cannet-des-Maures (2017). We then performed simulations of the model using these estimates to forecast the magnitude of the foci of transmission as a function of the response delay and the moment of virus introduction.ConclusionsThe results of the different simulations underline the relative importance of each variable and can be useful to stakeholders when designing context-based intervention strategies. The findings emphasize the importance of, and advocate for early detection of imported cases and timely biological confirmation of autochthonous cases to ensure timely vector control measures, supporting the implementation and the maintenance of sustainable surveillance systems.     

Peptide Partitioning Properties from Direct Insertion Studies

Partitioning properties of polypeptides are at the heart of biological membrane phenomena and their precise quantification is vital for ab-initio structure prediction and the accurate simulation of membrane protein folding and function. Recently the cellular translocon machinery has been employed to determine membrane insertion propensities and transfer energetics for a series of polyleucine segments embedded in a carrier sequence. We show here that the insertion propensity, pathway, and transfer energetics into synthetic POPC bilayers can be fully described by direct atomistic peptide partitioning simulations. The insertion probability as a function of peptide length follows two-state Boltzmann statistics, in agreement with the experiments. The simulations expose a systematic offset between translocon-mediated and direct insertion free energies. Compared to the experiment the insertion threshold is shifted toward shorter peptides by ∼2 leucine residues. The simulations reveal many hitherto unknown atomic-resolution details about the partitioning process and promise to provide a powerful tool for urgently needed calibration of lipid parameters to match experimentally observed peptide transfer energies.     

Computer simulations of cell-target encounter including biased cell motion toward targets: Single and multiple cell-target simulations in two dimensions

In order for immune cells to carry out many of their functions, including clearance of infectious agents from tissue, they must first encounter their targets in the tissue. This encounter process is often the rate-limiting step in the overall function. Most immune cells exhibit chemotactic ability, and previous continuum models for encounter rates and dynamics have shown that chemotaxis can be a great advantage to cells by greatly increasing encounter rates relative to those for randomly moving cells. This paper describes computer simulations of discrete cell-target encounter events in two dimensions, for the two cases considered by the continuum models: where only a single cell and a single target are present, and where many cells and targets are present. The results of these simulations verify our previous model predictions that a small amount of chemotactic bias dramatically decreases the encounter time, while further increases in the amount of bias have a much smaller effect. Chemotactic ability is shown to be an important determinant of the kinetics of target clearance, and its effects depend on the initial cell-target ratio and the initial distributions of cells and targets. To the best of our knowledge, this work provides the first computer simulations of particle-target encounter in which there is biased motion of particles toward their targets, and is therefore of general interest beyond specific application to immune cell function.     

Molecular Dynamics of Channelrhodopsin at the Early Stages of Channel Opening

Channelrhodopsin (ChR) is a light-gated cation channel that responds to blue light. Since ChR can be readily expressed in specific neurons to precisely control their activities by light, it has become a powerful tool in neuroscience. Although the recently solved crystal structure of a chimeric ChR, C1C2, provided the structural basis for ChR, our understanding of the molecular mechanism of ChR still remains limited. Here we performed electrophysiological analyses and all-atom molecular dynamics (MD) simulations, to investigate the importance of the intracellular and central constrictions of the ion conducting pore observed in the crystal structure of C1C2. Our electrophysiological analysis revealed that two glutamate residues, Glu122 and Glu129, in the intracellular and central constrictions, respectively, should be deprotonated in the photocycle. The simulation results suggested that the deprotonation of Glu129 in the central constriction leads to ion leakage in the ground state, and implied that the protonation of Glu129 is important for preventing ion leakage in the ground state. Moreover, we modeled the 13-cis retinal bound; i.e., activated C1C2, and performed MD simulations to investigate the conformational changes in the early stage of the photocycle. Our simulations suggested that retinal photoisomerization induces the conformational change toward channel opening, including the movements of TM6, TM7 and TM2. These insights into the dynamics of the ground states and the early photocycle stages enhance our understanding of the channel function of ChR.     

Modeling the Hydration Layer around Proteins: HyPred

Protein hydration plays an integral role in determining protein function and stability. We develop a simple method with atomic level precision for predicting the solvent density near the surface of a protein. A set of proximal radial distribution functions are defined and calculated for a series of different atom types in proteins using all-atom, explicit solvent molecular dynamic simulations for three globular proteins. A major improvement in predicting the hydration layer is found when the protein is held immobile during the simulations. The distribution functions are used to develop a model for predicting the hydration layer with sub-1-Ångstrom resolution without the need for additional simulations. The model and the distribution functions for a given protein are tested in their ability to reproduce the hydration layer from the simulations for that protein, as well as those for other proteins and for simulations in which the protein atoms are mobile. Predictions for the density of water in the hydration shells are then compared with high occupancy sites observed in crystal structures. The accuracy of both tests demonstrates that the solvation model provides a basis for quantitatively understanding protein solvation and thereby predicting the hydration layer without additional simulations.     

Brownian dynamics simulations of probe and self-diffusion in concentrated protein and DNA solutions.

We have developed a Brownian dynamics algorithm for simulating probe and self-diffusion in concentrated solutions of DNA and protein. In these simulations, proteins are represented as spheres with radii given by their hydrodynamic radii, while DNA is modeled as a wormlike chain of hydrodynamically equivalent spherical frictional elements. The molecular interaction potentials employed by the program allow for intramolecular stretching and bending motions of the DNA chains, short-range Lennard-Jones interactions, and long-range electrostatic interactions. To test the program, we have carried out simulations of bovine serum albumin (BSA) probe diffusion and DNA self-diffusion in solutions of short-chain DNA as a function of both DNA concentration and solution ionic strength. In addition, we report on simulations of BSA self-diffusion as a function of BSA concentration and ionic strength. Based on a comparison to available experimental data, we find that our simulations accurately predict these transport properties under conditions of physiological salt concentration and predict the stronger concentration dependence observed at lower salt concentrations. These results are discussed in light of the nature of the intermolecular interactions in such systems and the approximations and limitations of the simulation algorithm.     

A model of stem cell population dynamics: in silico analysis and in vivo validation

The proper renewal and maintenance of tissues by stem cell populations is simultaneously influenced by anatomical constraints, cell proliferation dynamics and cell fate specification. However, their relative influence is difficult to examine in vivo. To address this difficulty we built, as a test case, a cell-centered state-based computational model of key behaviors that govern germline development in C. elegans, and used it to drive simulations of cell population dynamics under a variety of perturbations. Our analysis provided unexpected possible explanations for laboratory observations, including certain 'all-or-none' phenotypes and complex differentiation patterns. The simulations also offered insights into niche-association dynamics and the interplay between cell cycle and cell fate. Subsequent experiments validated several predictions generated by the simulations. Notably, we found that early cell cycle defects influence later maintenance of the progenitor cell population. This general modeling approach is potentially applicable to other stem cell systems.     

Visualization of early events in acetic acid denaturation of HIV-1 protease: a molecular dynamics study

Protein denaturation plays a crucial role in cellular processes. In this study, denaturation of HIV-1 Protease (PR) was investigated by all-atom MD simulations in explicit solvent. The PR dimer and monomer were simulated separately in 9 M acetic acid (9 M AcOH) solution and water to study the denaturation process of PR in acetic acid environment. Direct visualization of the denaturation dynamics that is readily available from such simulations has been presented. Our simulations in 9 M AcOH reveal that the PR denaturation begins by separation of dimer into intact monomers and it is only after this separation that the monomer units start denaturing. The denaturation of the monomers is flagged off by the loss of crucial interactions between the α-helix at C-terminal and surrounding β-strands. This causes the structure to transit from the equilibrium dynamics to random non-equilibrating dynamics. Residence time calculations indicate that denaturation occurs via direct interaction of the acetic acid molecules with certain regions of the protein in 9 M AcOH. All these observations have helped to decipher a picture of the early events in acetic acid denaturation of PR and have illustrated that the α-helix and the β-sheet at the C-terminus of a native and functional PR dimer should maintain both the stability and the function of the enzyme and thus present newer targets for blocking PR function.     

The role of oligomerization and cooperative regulation in protein function: the case of tryptophan synthase

The oligomerization/co-localization of protein complexes and their cooperative regulation in protein function is a key feature in many biological systems. The synergistic regulation in different subunits often enhances the functional properties of the multi-enzyme complex. The present study used molecular dynamics and Brownian dynamics simulations to study the effects of allostery, oligomerization and intermediate channeling on enhancing the protein function of tryptophan synthase (TRPS). TRPS uses a set of α/β-dimeric units to catalyze the last two steps of L-tryptophan biosynthesis, and the rate is remarkably slower in the isolated monomers. Our work shows that without their binding partner, the isolated monomers are stable and more rigid. The substrates can form fairly stable interactions with the protein in both forms when the protein reaches the final ligand-bound conformations. Our simulations also revealed that the α/β-dimeric unit stabilizes the substrate-protein conformation in the ligand binding process, which lowers the conformation transition barrier and helps the protein conformations shift from an open/inactive form to a closed/active form. Brownian dynamics simulations with a coarse-grained model illustrate how protein conformations affect substrate channeling. The results highlight the complex roles of protein oligomerization and the fine balance between rigidity and dynamics in protein function.     

Folding simulations of small proteins

Understanding how a protein folds is a long-standing challenge in modern science. We have used an optimized atomistic model (united-residue force field) to simulate folding of small proteins of various structures: HP-36 (alpha protein), protein A (beta), 1fsd (alpha+beta), and betanova (beta). Extensive Monte Carlo folding simulations (ten independent runs with 10(9) Monte Carlo steps at a temperature) starting from non-native conformations are carried out for each protein. In all cases, proteins fold into their native-like conformations at appropriate temperatures, and glassy transitions occur at low temperatures. To investigate early folding trajectories, 200 independent runs with 10(6) Monte Carlo steps are also performed at a fixed temperature for a protein. There are a variety of possible pathways during non-equilibrium early processes (fast process, approximately 10(4) Monte Carlo steps). Finally, these pathways converge to the point unique for each protein. The convergence point of the early folding pathways can be determined only by direct folding simulations. The free energy surface, an equilibrium thermodynamic property, dictates the rest of the folding (slow process, approximately 10(8) Monte Carlo steps).     

Quantitative Analysis of the Timing of the Origin and Diversification of Extant Placental Orders

Fossil evidence is consistent with origination and diversification of extant placental orders in the early Tertiary (Explosive Model), and with the possibility of some orders having stem taxa extending into the Cretaceous (Long Fuse Model). Fossil evidence that 15 of 18 extant placental orders appeared and began diversification in the first 16 m.y. of the Cenozoic is, however, at odds with molecular studies arguing some orders diversified up to 40 m.y. earlier in the Early Cretaceous (Short Fuse Model). The quality of the fossil record was assessed by tabulating localities of all mammals in the last 105 m.y. Global locality data (except Africa) for 105 m.y. of eutherian evolution indicate discernible biogeographic patterns by the last 15 m.y. of the Cretaceous. Eutherian genera increase from 11 in latest Cretaceous to 139 in earliest Tertiary, although both are represented by about 50 localities. Yet even in the Late Cretaceous of North America and Asia where eutherians are abundant, none of the 18 extant orders are definitely known. A series of Monte Carlo simulations test whether the rapid appearance of most mammalian orders is statistically significant, and if so, whether it is a radiation event or an artifact of a limited fossil record. Monte Carlo tests affirm that the clustering of appearances in the early Cenozoic is statistically significant. Quantitative analysis of the locality data suggests that the number of genera described is a function of the number of localities sampled. In contrast, the number of orders is not a simple function of localities and thus does not appear to be limited by localities. A second set of Monte Carlo simulations confirms that the increase in orders cannot be explained by the limited number of localities sampled. Even for best-fit simulations, the observed pattern of ordinal appearances is steeper than expected under a variety of null models. These quantitative analyses of the fossil record demonstrate that the rapid ordinal appearances cannot be ascribed to limited Late Cretaceous sample sizes; thus, early Tertiary ordinal diversification is real. Although the fossil record is incomplete, it appears adequate to reject the hypothesis that orders of placentals began to diversify before the K/T boundary.     

Some theoretical considerations on cytosolic redox balance during anaerobiosis in marine invertebrates

Studies on anaerobiosis in marine invertebrates have shown that many rely on malate, octopine, or alanopine dehydrogenases, rather than lactate dehydrogenase, for cytosolic redox balance. These systems were studied by computer simulations with the assumption that these dehydrogenases maintain their substrates and products at instantaneous equilibrium. The simulations permit a study of the redox ratio (NADH/NAD⁺) as a function of the concentration of lactate, malate, or octopine. The redox ratio was found to increase as these products accumulated. It was substantially less in the simulations of malate and octopine dehydrogenases when compared to those of lactate dehydrogenase. This factor may be important for maintaining glycolysis in these organisms, and suggests an advantage for the use of octopine dehydrogenase rather than its analogue lactate dehydrogenase.     

Muscle and prosthesis contributions to amputee walking mechanics: A modeling study

Unilateral, below-knee amputees have altered gait mechanics, which can significantly affect their mobility. Below-knee amputees lose the functional use of the ankle muscles, which are critical during walking to provide body support, forward propulsion, leg-swing initiation and mediolateral balance. Thus, either muscles must compensate or the prosthesis must provide the functional tasks normally provided by the ankle muscles. Three-dimensional (3D) forward dynamics simulations of amputee and non-amputee walking were generated to identify muscle and prosthesis contributions to amputee walking mechanics, including the subtasks of body support, forward propulsion, leg-swing initiation and mediolateral balance. Results showed that the prosthesis provided body support in the absence of the ankle muscles. The prosthesis contributed to braking from early to mid-stance and propulsion in late stance. The prosthesis also functioned like the uniarticular soleus muscle by transferring energy from the residual leg to the trunk to provide trunk propulsion. The residual-leg vasti and rectus femoris reduced their contributions to braking in early stance, which mitigated braking from the prosthesis during this period. The prosthesis did not replace the function of the gastrocnemius, which normally generates energy to the leg to initiate swing. As a result, lower overall energy was delivered to the residual leg. The prosthesis also acted to accelerate the body laterally in the absence of the ankle muscles. These results provide further insight into muscle and prosthesis function in below-knee amputee walking and can help guide rehabilitation methods and device designs to improve amputee mobility.     

The molecular origin of the MMR-dependent apoptosis pathway from dynamics analysis of MutSα-DNA complexes

The cellular response to DNA damage signaling by mismatch-repair (MMR) proteins is incompletely understood. It is generally accepted that MMR-dependent apoptosis pathway in response to DNA damage detection is independent of MMR's DNA repair function. In this study, we investigate correlated motions in response to the binding of mismatched and platinum cross-linked DNA fragments by MutSα, as derived from 50 ns molecular dynamics simulations. The protein dynamics in response to the mismatched and damaged DNA recognition suggests that MutSα signals their recognition through independent pathways providing evidence for the molecular origin of the MMR-dependent apoptosis. MSH2 subunit is indicated to play a key role in signaling both mismatched and damaged DNA recognition; localized and collective motions within the protein allow identifying sites on the MSH2 surface possible involved in recruiting proteins responsible for downstream events. Unlike in the mismatch complex, predicted key communication sites specific for the damage recognition are on the list of known cancer-causing mutations or deletions. This confirms MSH2's role in signaling DNA damage-induced apoptosis and suggests that defects in MMR alone is sufficient to trigger tumorigenesis, supporting the experimental evidence that MMR-damage response function could protect from the early occurrence of tumors. Identifying these particular communication sites may have implications for the treatment of cancers that are not defective for MMR, but are unable to function optimally for MMR-dependent responses following DNA damage such as the case of resistance to cisplatin.     

Native structure-based modeling and simulation of biomolecular systems per mouse click

Background

Molecular dynamics (MD) simulations provide valuable insight into biomolecular systems at the atomic level. Notwithstanding the ever-increasing power of high performance computers current MD simulations face several challenges: the fastest atomic movements require time steps of a few femtoseconds which are small compared to biomolecular relevant timescales of milliseconds or even seconds for large conformational motions. At the same time, scalability to a large number of cores is limited mostly due to long-range interactions. An appealing alternative to atomic-level simulations is coarse-graining the resolution of the system or reducing the complexity of the Hamiltonian to improve sampling while decreasing computational costs. Native structure-based models, also called Gō-type models, are based on energy landscape theory and the principle of minimal frustration. They have been tremendously successful in explaining fundamental questions of, e.g., protein folding, RNA folding or protein function. At the same time, they are computationally sufficiently inexpensive to run complex simulations on smaller computing systems or even commodity hardware. Still, their setup and evaluation is quite complex even though sophisticated software packages support their realization.

Results

Here, we establish an efficient infrastructure for native structure-based models to support the community and enable high-throughput simulations on remote computing resources via GridBeans and UNICORE middleware. This infrastructure organizes the setup of such simulations resulting in increased comparability of simulation results. At the same time, complete workflows for advanced simulation protocols can be established and managed on remote resources by a graphical interface which increases reusability of protocols and additionally lowers the entry barrier into such simulations for, e.g., experimental scientists who want to compare their results against simulations. We demonstrate the power of this approach by illustrating it for protein folding simulations for a range of proteins.

Conclusions

We present software enhancing the entire workflow for native structure-based simulations including exception-handling and evaluations. Extending the capability and improving the accessibility of existing simulation packages the software goes beyond the state of the art in the domain of biomolecular simulations. Thus we expect that it will stimulate more individuals from the community to employ more confidently modeling in their research.

Electronic supplementary material

The online version of this article (doi:10.1186/1471-2105-15-292) contains supplementary material, which is available to authorized users.     

DNA damage in synchronized HeLa cells irradiated with ultraviolet.

The lethal effect of UV radiation of HeLa cells is least in mitosis and greatest in late G1-early S. Photochemical damage to HeLa DNA, as measured by thymine-containing dimer formation and by alkaline sucrose sedimentation, also increases from mitosis towards early S phase. Computer simulations of UV absorption by an idealized HeLa cell at different stages of the cell cycle indicate that changes in damage could be due solely to changes in chromatin geometry. But survival is not exclusively a function of damage.     

Distinct allosteric pathways in imidazole glycerol phosphate synthase from yeast and bacteria

Understanding the relationship between protein structures and their function is still an open question that becomes very challenging when allostery plays an important functional role. Allosteric proteins, in fact, exploit different ranges of motions (from sidechain local fluctuations to long-range collective motions) to effectively couple distant binding sites, and of particular interest is whether allosteric proteins of the same families with similar functions and structures also necessarily share the same allosteric mechanisms. Here, we compared the early dynamics initiating the allosteric communication of a prototypical allosteric enzyme from two different organisms, i.e., the imidazole glycerol phosphate synthase (IGPS) enzymes from the thermophilic bacteria and the yeast, working at high and room temperatures, respectively. By combining molecular dynamics simulations and network models derived from graph theory, we found rather distinct early allosteric dynamics in the IGPS from the two organisms, involving significatively different allosteric pathways in terms of both local and collective motions. Given the successful prediction of key allosteric residues in the bacterial IGPS, whose mutation disrupts its allosteric communication, the outcome of this study paves the way for future experimental studies on the yeast IGPS that could foster therapeutic applications by exploiting the control of IGPS enzyme allostery.