fastsimcoal: a continuous-time coalescent simulator of genomic diversity under arbitrarily complex evolutionary scenarios

MOTIVATION: Genetic studies focus on increasingly larger genomic regions of both extant and ancient DNA, and there is a need for simulation software to match these technological advances. We present here a new coalescent-based simulation program fastsimcoal, which is able to quickly simulate a variety of genetic markers scattered over very long genomic regions with arbitrary recombination patterns under complex evolutionary scenarios. Availability and Implementation: fastsimcoal is a C++ program compiled for Windows, MacOsX and Linux platforms. It is freely available at cmpg.unibe.ch/software/fastsimcoal/, together with its detailed user manual and example input files.     

A Software Architecture for Multi-Cellular System Simulations on Graphics Processing Units

The first aim of simulation in virtual environment is to help biologists to have a better understanding of the simulated system. The cost of such simulation is significantly reduced compared to that of in vivo simulation. However, the inherent complexity of biological system makes it hard to simulate these systems on non-parallel architectures: models might be made of sub-models and take several scales into account; the number of simulated entities may be quite large. Today, graphics cards are used for general purpose computing which has been made easier thanks to frameworks like CUDA or OpenCL. Parallelization of models may however not be easy: parallel computer programing skills are often required; several hardware architectures may be used to execute models. In this paper, we present the software architecture we built in order to implement various models able to simulate multi-cellular system. This architecture is modular and it implements data structures adapted for graphics processing units architectures. It allows efficient simulation of biological mechanisms.     

Automated molecular library generation of proteic fragments by virtual proteolysis for molecular modelling studies

This paper presents a computer aided design method useful for simulation of a set of proteolytic cleavages upon target proteins obtained from the Brookhaven Data Bank. The method was developed by using algorithms that are able to interface themselves with other software environments, in order to assist computer analyses in the molecular modelling field, and allowing the generation of molecular libraries containing protein fragments produced by simulated proteolysis. These libraries include structures that differ for several amino acid deletions upon specified regions of the primary sequence. Target residues chosen for the simulation are compatible with enzymatic proteolysis methods used in conventional laboratory procedures. Furthermore, algorithms were able to identify a set of chemical-physical properties of the starting proteins, leading the simulation to find out the most suitable residues for proteolysis. The goal of these strategies is to generate fragments that are leaded to maintain the native-like condition of starting molecules, avoiding loss of conformational characteristics of the original tertiary structure. Proteins chosen for generating proteolytic libraries were represented by naphthalene 1,2 dioxygenase and Rigidoporus lignosus laccase.     

Diffraction-based density restraints for membrane and membrane-peptide molecular dynamics simulations

We have recently shown that current molecular dynamics (MD) atomic force fields are not yet able to produce lipid bilayer structures that agree with experimentally-determined structures within experimental errors. Because of the many advantages offered by experimentally validated simulations, we have developed a novel restraint method for membrane MD simulations that uses experimental diffraction data. The restraints, introduced into the MD force field, act upon specified groups of atoms to restrain their mean positions and widths to values determined experimentally. The method was first tested using a simple liquid argon system, and then applied to a neat dioleoylphosphatidylcholine (DOPC) bilayer at 66% relative humidity and to the same bilayer containing the peptide melittin. Application of experiment-based restraints to the transbilayer double-bond and water distributions of neat DOPC bilayers led to distributions that agreed with the experimental values. Based upon the experimental structure, the restraints improved the simulated structure in some regions while introducing larger differences in others, as might be expected from imperfect force fields. For the DOPC-melittin system, the experimental transbilayer distribution of melittin was used as a restraint. The addition of the peptide caused perturbations of the simulated bilayer structure, but which were larger than observed experimentally. The melittin distribution of the simulation could be fit accurately to a Gaussian with parameters close to the observed ones, indicating that the restraints can be used to produce an ensemble of membrane-bound peptide conformations that are consistent with experiments. Such ensembles pave the way for understanding peptide-bilayer interactions at the atomic level.     

Additional Efficient Computation of Branched Nerve Equations: Adaptive Time Step and Ideal Voltage Clamp

Various improvements are described for the simulation of biophysically and anatomically detailed compartmental models of single neurons and networks of neurons. These include adaptive time-step integration and a reordering of the circuit matrix to allow ideal voltage clamp of arbitrary nodes. We demonstrate how the adaptive time-step method can give equivalent accuracy as a fixed time-step method for typical current clamp simulation protocols, with about a 2.5 reduction in runtime. The ideal voltage clamp method is shown to be more stable than the nonideal case, in particular when used with the adaptive time-step method. Simulation results are presented using the Surf-Hippo Neuron Simulation System, a public domain object-oriented simulator written in Lisp.     

An unbiased adaptive sampling algorithm for the exploration of RNA mutational landscapes under evolutionary pressure

The analysis of the relationship between sequences and structures (i.e., how mutations affect structures and reciprocally how structures influence mutations) is essential to decipher the principles driving molecular evolution, to infer the origins of genetic diseases, and to develop bioengineering applications such as the design of artificial molecules. Because their structures can be predicted from the sequence data only, RNA molecules provide a good framework to study this sequence-structure relationship. We recently introduced a suite of algorithms called RNAmutants which allows a complete exploration of RNA sequence-structure maps in polynomial time and space. Formally, RNAmutants takes an input sequence (or seed) to compute the Boltzmann-weighted ensembles of mutants with exactly k mutations, and sample mutations from these ensembles. However, this approach suffers from major limitations. Indeed, since the Boltzmann probabilities of the mutations depend of the free energy of the structures, RNAmutants has difficulties to sample mutant sequences with low G+C-contents. In this article, we introduce an unbiased adaptive sampling algorithm that enables RNAmutants to sample regions of the mutational landscape poorly covered by classical algorithms. We applied these methods to sample mutations with low G+C-contents. These adaptive sampling techniques can be easily adapted to explore other regions of the sequence and structural landscapes which are difficult to sample. Importantly, these algorithms come at a minimal computational cost. We demonstrate the insights offered by these techniques on studies of complete RNA sequence structures maps of sizes up to 40 nucleotides. Our results indicate that the G+C-content has a strong influence on the size and shape of the evolutionary accessible sequence and structural spaces. In particular, we show that low G+C-contents favor the apparition of internal loops and thus possibly the synthesis of tertiary structure motifs. On the other hand, high G+C-contents significantly reduce the size of the evolutionary accessible mutational landscapes.     

不同气象数据空间内插对区域玉米生长模拟结果的影响

作物模型区域模拟已成为作物模型应用的一个新方向。运用作物模型进行区域研究时,遇到的问题之一就是输入模型的空间数据质量问题,研究不同空间内插法获得的气象数据对作物模型区域模拟结果的影响,可以为区域模拟对输入数据的敏感性研究提供一定的参考。利用区域校准的CERES-Maize模型,将3类内插方法(几何内插、统计内插、动力模型内插)产生的网格化天气数据分别输入到CERES-Maize模型中,模拟了50km×50km网格水平下1961—1990年我国玉米生产状况,并选取1980—1990年模拟的平均产量与同期农调队调查产量进行比较,以了解区域模拟中,不同空间内插方法所得的逐日气象数据对区域模拟结果的影响。结果表明:(1)作物模型区域应用时,所采用的3种内插方法都能满足作物模型区域模拟对网格化天气数据的要求,采用3种天气数据的区域模拟结果都能反映出玉米平均产量的空间变化特征,与网格调查平均产量之间具有极显著的相关关系,但采用不同内插天气数据对模拟结果造成了8%以内的偏差。(2)采用不同内插天气数据,在进行作物区域模拟时,各方法的模拟结果之间呈极显著的相关关系,但这些模拟结果之间,在全国大部分地区是差异显著。     

Compartmental neural simulations with spatial adaptivity

Since their inception, computational models have become increasingly complex and useful counterparts to laboratory experiments within the field of neuroscience. Today several software programs exist to solve the underlying mathematical system of equations, but such programs typically solve these equations in all parts of a cell (or network of cells) simultaneously, regardless of whether or not all of the cell is active. This approach can be inefficient if only part of the cell is active and many simulations must be performed. We have previously developed a numerical method that provides a framework for spatial adaptivity by making the computations local to individual branches rather than entire cells (Rempe and Chopp, SIAM Journal on Scientific Computing, 28: 2139–2161, 2006). Once the computation is reduced to the level of branches instead of cells, spatial adaptivity is straightforward: the active regions of the cell are detected and computational effort is focused there, while saving computations in other regions of the cell that are at or near rest. Here we apply the adaptive method to four realistic neuronal simulation scenarios and demonstrate its improved efficiency over non-adaptive methods. We find that the computational cost of the method scales with the amount of activity present in the simulation, rather than the physical size of the system being simulated. For certain problems spatial adaptivity reduces the computation time by up to 80%.     

REJECTOR: software for population history inference from genetic data via a rejection algorithm

We introduce REJECTOR, software for parameter estimation and comparison of alternate models of population history from genetic data via a rejection algorithm. Through coalescent simulation, REJECTOR generates numerous gene genealogies, and hence simulated data, under a model of population history specified by the user. Summary statistics derived from such simulated data are compared with observed statistics, leading to acceptance or rejection of a given set of parameter values. We performed tests of the software using known parameter values in order to assess the inferential power provided by each summary statistic. The tests demonstrate the precision and accuracy of estimation made possible using this approach. AVAILABILITY: http://www.rejector.org     

eXPatGen: generating dynamic expression patterns for the systematic evaluation of analytical methods

MOTIVATION: Experimental gene expression data sets, such as those generated by microarray or gene chip experiments, typically have significant noise and complicated interconnectivities that make understanding even simple regulatory patterns difficult. Given these complications, characterizing the effectiveness of different analysis techniques to uncover network groups and structures remains a challenge. Generating simulated expression patterns with known biological features of expression complexity, diversity and interconnectivities provides a more controlled means of investigating the appropriateness of different analysis methods. A simulation-based approach can systematically evaluate different gene expression analysis techniques and provide a basis for improved methods in dynamic metabolic network reconstruction. RESULTS: We have developed an on-line simulator, called eXPatGen, to generate dynamic gene expression patterns typical of microarray experiments. eXPatGen provides a quantitative network structure to represent key biological features, including the induction, repression, and cascade regulation of messenger RNA (mRNA). The simulation is modular such that the expression model can be replaced with other representations, depending on the level of biological detail required by the user. Two example gene networks, of 25 and 100 genes respectively, were simulated. Two standard analysis techniques, clustering and PCA analysis, were performed on the resulting expression patterns in order to demonstrate how the simulator might be used to evaluate different analysis methods and provide experimental guidance for biological studies of gene expression. AVAILABILITY: http://www.che.udel.edu/eXPatGen/     

indel-Seq-Gen: a new protein family simulator incorporating domains, motifs, and indels

Reconstructing the evolutionary history of protein sequences will provide a better understanding of divergence mechanisms of protein superfamilies and their functions. Long-term protein evolution often includes dynamic changes such as insertion, deletion, and domain shuffling. Such dynamic changes make reconstructing protein sequence evolution difficult and affect the accuracy of molecular evolutionary methods, such as multiple alignments and phylogenetic methods. Unfortunately, currently available simulation methods are not sufficiently flexible and do not allow biologically realistic dynamic protein sequence evolution. We introduce a new method, indel-Seq-Gen (iSG), that can simulate realistic evolutionary processes of protein sequences with insertions and deletions (indels). Unlike other simulation methods, iSG allows the user to simulate multiple subsequences according to different evolutionary parameters, which is necessary for generating realistic protein families with multiple domains. iSG tracks all evolutionary events including indels and outputs the "true" multiple alignment of the simulated sequences. iSG can also generate a larger sequence space by allowing the use of multiple related root sequences. With all these functions, iSG can be used to test the accuracy of, for example, multiple alignment methods, phylogenetic methods, evolutionary hypotheses, ancestral protein reconstruction methods, and protein family classification methods. We empirically evaluated the performance of iSG against currently available methods by simulating the evolution of the G protein-coupled receptor and lipocalin protein families. We examined their true multiple alignments, reconstruction of the transmembrane regions and beta-strands, and the results of similarity search against a protein database using the simulated sequences. We also presented an example of using iSG for examining how phylogenetic reconstruction is affected by high indel rates.     

Molecular dynamics in the endgame of protein structure prediction.

In order adequately to sample conformational space, methods for protein structure prediction make necessary simplifications that also prevent them from being as accurate as desired. Thus, the idea of feeding them, hierarchically, into a more accurate method that samples less effectively was introduced a decade ago but has not met with more than limited success in a few isolated instances. Ideally, the final stages should be able to identify the native state, show a good correlation with native similarity in order to add value to the selection process, and refine the structures even further. In this work, we explore the possibility of using state-of-the-art explicit solvent molecular dynamics and implicit solvent free energy calculations to accomplish all three of those objectives on 12 small, single-domain proteins, four each of alpha, beta and mixed topologies. We find that this approach is very successful in ranking the native and also enhances the structure selection of predictions generated from the Rosetta method.     

Simulation of bone adaptive remodeling using a stochastic process as loading history

In this simulation study for bone adaptive remodeling, loading conditions are described as stochastic processes to catch the unpredictable characteristics of daily physical activities, which are observed to be closely related with bone adaptive remodeling. This will not only eliminate the necessity of arbitrary choices for loading conditions, but also generate greater flexibility for simulations of bone adaptive remodeling. The sensitivity of simulation outcomes to the parameters in the simulation algorithm was examined by applying stochastic loading conditions on finite element models of simplified spine structures. In this way, the limitations induced by simplifying loading conditions into constant or cyclic loads can be avoided and, potentially, more clinical observations could be accommodated when more comprehensive finite element models are available.     

Evolution of taxis responses in virtual bacteria: non-adaptive dynamics

Bacteria are able to sense and respond to a variety of external stimuli, with responses that vary from stimuli to stimuli and from species to species. The best-understood is chemotaxis in the model organism Escherichia coli, where the dynamics and the structure of the underlying pathway are well characterised. It is not clear, however, how well this detailed knowledge applies to mechanisms mediating responses to other stimuli or to pathways in other species. Furthermore, there is increasing experimental evidence that bacteria integrate responses from different stimuli to generate a coherent taxis response. We currently lack a full understanding of the different pathway structures and dynamics and how this integration is achieved. In order to explore different pathway structures and dynamics that can underlie taxis responses in bacteria, we perform a computational simulation of the evolution of taxis. This approach starts with a population of virtual bacteria that move in a virtual environment based on the dynamics of the simple biochemical pathways they harbour. As mutations lead to changes in pathway structure and dynamics, bacteria better able to localise with favourable conditions gain a selective advantage. We find that a certain dynamics evolves consistently under different model assumptions and environments. These dynamics, which we call non-adaptive dynamics, directly couple tumbling probability of the cell to increasing stimuli. Dynamics that are adaptive under a wide range of conditions, as seen in the chemotaxis pathway of E. coli, do not evolve in these evolutionary simulations. However, we find that stimulus scarcity and fluctuations during evolution results in complex pathway dynamics that result both in adaptive and non-adaptive dynamics depending on basal stimuli levels. Further analyses of evolved pathway structures show that effective taxis dynamics can be mediated with as few as two components. The non-adaptive dynamics mediating taxis responses provide an explanation for experimental observations made in mutant strains of E. coli and in wild-type Rhodobacter sphaeroides that could not be explained with standard models. We speculate that such dynamics exist in other bacteria as well and play a role linking the metabolic state of the cell and the taxis response. The simplicity of mechanisms mediating such dynamics makes them a candidate precursor of more complex taxis responses involving adaptation. This study suggests a strong link between stimulus conditions during evolution and evolved pathway dynamics. When evolution was simulated under conditions of scarce and fluctuating stimulus conditions, the evolved pathway contained features of both adaptive and non-adaptive dynamics, suggesting that these two types of dynamics can have different advantages under distinct environmental circumstances.     

Controllers for imposing continuum-to-molecular boundary conditions in arbitrary fluid flow geometries

We present a new parallelised controller for steering an arbitrary geometric region of a molecular dynamics (MD) simulation towards a desired thermodynamic and hydrodynamic state. We show that the controllers may be applied anywhere in the domain to set accurately an initial MD state, or solely at boundary regions to prescribe non-periodic boundary conditions (PBCs) in MD simulations. The mean molecular structure and velocity autocorrelation function remain unchanged (when sampled a few molecular diameters away from the constrained region) when compared with those distributions measured using PBCs. To demonstrate the capability of our new controllers, we apply them as non-PBCs in parallel to a complex MD mixing nano-channel and in a hybrid MD continuum simulation with a complex coupling region. The controller methodology is easily extendable to polyatomic MD fluids.     

Construction and electrophoretic migration of single-stranded DNA knots and catenanes

In recent years there has been growing interest in the question of how the particular topology of polymeric chains affects their overall dimensions and physical behavior. The majority of relevant studies are based on numerical simulation methods or analytical treatment; however, both these approaches depend on various assumptions and simplifications. Experimental verification is clearly needed but was hampered by practical difficulties in obtaining preparative amounts of knotted or catenated polymers with predefined topology and precisely set chain length. We introduce here an efficient method of production of various single-stranded DNA knots and catenanes that have the same global chain length. We also characterize electrophoretic migration of the produced single-stranded DNA knots and catenanes with increasing complexity.     

Adaptive tests for detecting gene-gene and gene-environment interactions

There has been an increasing interest in detecting gene-gene and gene-environment interactions in genetic association studies. A major statistical challenge is how to deal with a large number of parameters measuring possible interaction effects, which leads to reduced power of any statistical test due to a large number of degrees of freedom or high cost of adjustment for multiple testing. Hence, a popular idea is to first apply some dimension reduction techniques before testing, while another is to apply only statistical tests that are developed for and robust to high-dimensional data. To combine both ideas, we propose applying an adaptive sum of squared score (SSU) test and several other adaptive tests. These adaptive tests are extensions of the adaptive Neyman test [Fan, 1996], which was originally proposed for high-dimensional data, providing a simple and effective way for dimension reduction. On the other hand, the original SSU test coincides with a version of a test specifically developed for high-dimensional data. We apply these adaptive tests and their original nonadaptive versions to simulated data to detect interactions between two groups of SNPs (e.g. multiple SNPs in two candidate regions). We found that for sparse models (i.e. with only few non-zero interaction parameters), the adaptive SSU test and its close variant, an adaptive version of the weighted sum of squared score (SSUw) test, improved the power over their non-adaptive versions, and performed consistently well across various scenarios. The proposed adaptive tests are built in the general framework of regression analysis, and can thus be applied to various types of traits in the presence of covariates.     

Object Command Language: a formalism to build molecular models and to analyze structural parameters in macromolecules,with applications to nucleic acids

We have written a programming language OCL (Object Command Language) to solve, in a general way, two recurring problems that arise during the construction of molecular models and during the geometrical characterization of macromolecules: how to move precisely and reproducibly any part of a molecular model in any user-defined local reference axes; and how to calculate standard or userdefined structural parameters that characterize the complex geometries of any macromolecule. OCL endows the user with three main capabilities: the definition of subsets of the macromolecule, called objects in OCL, with a formalism from elementary set theory or lexical analysis; the definition of sequences of elementary geometrical operations, called procedures in OCL, enabling one to build arbitrary three-dimensional (3D) orthonormal reference frames, to be associated with previously defined objects; and the transmission of these definitions to programs that allow one to display, to modify and to analyze interactively the molecular structure, or to programs that perform energy minimizations or molecular dynamics. Several applications to nucleic acids are presented.     

AClAP, Autonomous hierarchical agglomerative Cluster Analysis based protocol to partition conformational datasets

MOTIVATION: Sampling the conformational space is a fundamental step for both ligand- and structure-based drug design. However, the rational organization of different molecular conformations still remains a challenge. In fact, for drug design applications, the sampling process provides a redundant conformation set whose thorough analysis can be intensive, or even prohibitive. We propose a statistical approach based on cluster analysis aimed at rationalizing the output of methods such as Monte Carlo, genetic, and reconstruction algorithms. Although some software already implements clustering procedures, at present, a universally accepted protocol is still missing. RESULTS: We integrated hierarchical agglomerative cluster analysis with a clusterability assessment method and a user independent cutting rule, to form a global protocol that we implemented in a MATLAB metalanguage program (AClAP). We tested it on the conformational space of a quite diverse set of drugs generated via Metropolis Monte Carlo simulation, and on the poses we obtained by reiterated docking runs performed by four widespread programs. In our tests, AClAP proved to remarkably reduce the dimensionality of the original datasets at a negligible computational cost. Moreover, when applied to the outcomes of many docking programs together, it was able to point to the crystallographic pose. AVAILABILITY: AClAP is available at the "AClAP" section of the website http://www.scfarm.unibo.it.     

On adaptive accuracy and precision in natural populations

Adaptation is usually conceived as the fit of a population mean to a fitness optimum. Natural selection, however, does not act only to optimize the population mean. Rather, selection normally acts on the fitness of individual organisms in the population. Furthermore, individual genotypes do not produce invariant phenotypes, and their fitness depends on how precisely they are able to realize their target phenotypes. For these reasons we suggest that it is better to conceptualize adaptation as accuracy rather than as optimality. The adaptive inaccuracy of a genotype can be measured as a function of the expected distance of its associated phenotype from a fitness optimum. The less the distance, the more accurate is the adaptation. Adaptive accuracy has two components: the deviance of the genotypically set target phenotype from the optimum and the precision with which this target phenotype can be realized. The second component, the adaptive precision, has rarely been quantified as such. We survey the literature to quantify how much of the phenotypic variation in wild populations is due to imprecise development. We find that this component is often substantial and highly variable across traits. We suggest that selection for improved precision may be important for many traits.