lea (likelihood‐based estimation of admixture): a program to estimate simultaneously admixture and time since the admixture event

We consider an admixture event, T generations in the past, where two ‘parental’ populations, P₁ and P₂, of size N₁ and N₂, contribute different proportions into the gene pool of an admixed population, H of size N_h. lea (likelihood‐based estimator of admixture) is a program which allows the user to obtain the posterior distribution of the parameters of the model. This includes p₁, the contribution of P₁, and t₁, t₂ and t_h, the time since the admixture event (scaled by the population size) for the three populations. lea allows the user to stop and restart the analyses at any time.     

Gibbs Ensemble Calculations with an Equation of State: An Application to Vapor-Liquid Equilibria

Abstract

A new modification of the Gibbs ensemble Monte Carlo computer simulation method for fluid phase equilibria is described. The modification is based on a thermodynamic model for the vapor phase, and uses an equation of state to account for the weak interactions between the vapor phase molecules. Reductions in the computational time by 30–40% as compared to the original Gibbs ensemble method are obtained. The algorithm is applied to Lennard-Jones - (12,6) fluids and their mixtures and the results are in good agreement with results obtained from simulations using the full Gibbs ensemble method.     

A comparison of alternative methods to compute conditional genotype probabilities for genetic evaluation with finite locus models

An increased availability of genotypes at marker loci has prompted the development of models that include the effect of individual genes. Selection based on these models is known as marker-assisted selection (MAS). MAS is known to be efficient especially for traits that have low heritability and non-additive gene action. BLUP methodology under non-additive gene action is not feasible for large inbred or crossbred pedigrees. It is easy to incorporate non-additive gene action in a finite locus model. Under such a model, the unobservable genotypic values can be predicted using the conditional mean of the genotypic values given the data. To compute this conditional mean, conditional genotype probabilities must be computed. In this study these probabilities were computed using iterative peeling, and three Markov chain Monte Carlo (MCMC) methods – scalar Gibbs, blocking Gibbs, and a sampler that combines the Elston Stewart algorithm with iterative peeling (ESIP). The performance of these four methods was assessed using simulated data. For pedigrees with loops, iterative peeling fails to provide accurate genotype probability estimates for some pedigree members. Also, computing time is exponentially related to the number of loci in the model. For MCMC methods, a linear relationship can be maintained by sampling genotypes one locus at a time. Out of the three MCMC methods considered, ESIP, performed the best while scalar Gibbs performed the worst.     

Reconciliation of different simulation methods in the determination of the equilibrium branch for adsorption in pores

The recently proposed mid-density scheme [Liu Z, Herrera L, Nguyen VT, Do DD, Nicholson D. A Monte Carlo scheme based on mid-density in a hysteresis loop to determine equilibrium phase transition. Mol Simul. 2011; 37(11):932–939, Liu Z, Do DD, Nicholson D. A thermodynamic study of the mid-density scheme to determine the equilibrium phase transition in cylindrical pores. Mol Simul. 2012; 38(3):189–199] is tested against a method 2V-NVT (similar to the well-established gauge cell method) and the canonical ensemble (CE) method, using argon adsorption at 87 K in graphitic slit pores of infinite and finite length. In infinitely long pores, the equilibrium transition is vertical that is expected for an infinite system to have a first-order transition and this vertical transition was found to lie at the middle of the hysteresis loop and satisfies the well-known Maxwell rule of equal area. For pores of finite length, the equilibrium transitions are steep and are close to, but not exactly identical to, the desorption branch. This lends support to the conventional view that the desorption branch is nearest to equilibrium, although both adsorption and desorption branches are strictly speaking metastable; a view proposed originally by Everett [Everett DH. Capillary condensation and adsorption hysteresis. Berichte Der Bunsen-Gesellschaft [Phys Chem Chem Phys]. 1975; 79(9):732–734]. As a consequence, the Maxwell rule of equal area does not apply to finite systems. As the widely accepted CE and gauge cell methods do not falsify the mid-density scheme, this study lends strong support to the validity of this technique for the study of equilibria.     

Differential dose contributions on total dose distribution of 125I brachytherapy source

This work provides an improvement of the approach using Monte Carlo simulation for the Amersham Model 6711 ¹²⁵I brachytherapy seed source, which is well known by many theoretical and experimental studies. The source which has simple geometry was researched with respect to criteria of AAPM Tg-43 Report. The approach offered by this study involves determination of differential dose contributions that come from virtual partitions of a massive radioactive element of the studied source to a total dose at analytical calculation point. Some brachytherapy seeds contain multi-radioactive elements so the dose at any point is a total of separate doses from each element. It is momentous to know well the angular and radial dose distributions around the source that is located in cancerous tissue for clinical treatments. Interior geometry of a source is effective on dose characteristics of a distribution. Dose information of inner geometrical structure of a brachytherapy source cannot be acquired by experimental methods because of limits of physical material and geometry in the healthy tissue, so Monte Carlo simulation is a required approach of the study. EGSnrc Monte Carlo simulation software was used. In the design of a simulation, the radioactive source was divided into 10 rings, partitioned but not separate from each other. All differential sources were simulated for dose calculation, and the shape of dose distribution was determined comparatively distribution of a single-complete source. In this work anisotropy function was examined also mathematically.     

Estimating the location and shape of hybrid zones

We propose a new model to make use of georeferenced genetic data for inferring the location and shape of a hybrid zone. The model output includes the posterior distribution of a parameter that quantifies the width of the hybrid zone. The model proposed is implemented in the GUI and command-line versions of the Geneland program versions ≥ 3.3.0. Information about the program can be found on http://www2.imm.dtu.dk/gigu/Geneland/.     

Computation of identity by descent probabilities conditional on DNA markers via a Monte Carlo Markov Chain method

The accurate estimation of the probability of identity by descent (IBD) at loci or genome positions of interest is paramount to the genetic study of quantitative and disease resistance traits. We present a Monte Carlo Markov Chain method to compute IBD probabilities between individuals conditional on DNA markers and on pedigree information. The IBDs can be obtained in a completely general pedigree at any genome position of interest, and all marker and pedigree information available is used. The method can be split into two steps at each iteration. First, phases are sampled using current genotypic configurations of relatives and second, crossover events are simulated conditional on phases. Internal track is kept of all founder origins and crossovers such that the IBD probabilities averaged over replicates are rapidly obtained. We illustrate the method with some examples. First, we show that all pedigree information should be used to obtain line origin probabilities in F2 crosses. Second, the distribution of genetic relationships between half and full sibs is analysed in both simulated data and in real data from an F2 cross in pigs.     

Simulation and measurement of microbeam dose distribution in lung tissue

Microbeam radiation therapy (MRT), a so far preclinical method in radiation oncology, modulates treatment doses on a micrometre scale. MRT uses treatment fields with a few ten micrometre wide high dose regions (peaks) separated by a few hundred micrometre wide low dose regions (valleys) and was shown to spare tissue much more effectively than conventional radiation therapy at similar tumour control rates. While preclinical research focused primarily on tumours of the central nervous system, recently also lung tumours have been suggested as a potential target for MRT.This study investigates the effect of the lung microstructure, comprising air cavities of a few hundred micrometre diameter, on the microbeam dose distribution in lung. In Monte Carlo simulations different models of heterogeneous lung tissue are compared with pure water and homogeneous air–water mixtures. Experimentally, microbeam dose distributions in porous foam material with cavity sizes similar to the size of lung alveoli were measured with film dosimetry at the European Synchrotron Radiation Facility (ESRF) in Grenoble, France.Simulations and experiments show that the microstructure of the lung has a huge impact on the local doses in the microbeam fields. Locally, material inhomogeneities may change the dose by a factor of 1.7, and also average peak and valley doses substantially differ from those in homogeneous material.Our results imply that accurate dose prediction for MRT in lung requires adequate models of the lung microstructure. Even if only average peak and valley doses are of interest, the assumption of a simple homogeneous air–water mixture is not sufficient. Since anatomic information on a micrometre scale are unavailable for clinical treatment planning, alternative methods and models have to be developed.     

Development of a dose distribution shifter to fit inside the collimator of a Boron Neutron Capture Therapy irradiation system to treat superficial tumours

The Kansai BNCT Medical Center has a cyclotron based epithermal neutron source for clinical Boron Neutron Capture Therapy. The system accelerates a proton to an energy of 30 MeV which strikes a beryllium target producing fast neutrons which are moderated down to epithermal neutrons for BNCT use. While clinical studies in the past have shown BNCT to be highly effective for malignant melanoma of the skin, to apply BNCT for superficial lesions using this system it is necessary to shift the thermal neutron distribution so that the maximum dose occurs near the surface. A dose distribution shifter was designed to fit inside the collimator to further moderate the neutrons to increase the surface dose and reduce the dose to the underlying normal tissue. Pure polyethylene was selected, and a Monte Carlo simulation was performed to determine the optimum thickness of the polyethylene slab. Compared with the original neutron beam, the shifter increased the thermal neutron flux at the skin by approximately 4 times. The measured and simulated central axis depth distribution and off axis distribution of the thermal neutron flux were found to be in good agreement. Compared with a 2 cm thick water equivalent bolus, a 26% increase in the thermal neutron flux at the surface was obtained, which would reduce the treatment time by approximately 29%. The DDS is a safe, simple and an effective tool for the treatment of superficial tumours for BNCT if an initially fast neutron beam requires moderation to maximise the thermal neutron flux at the tissue surface.     

Use of quantitative structure-property relationships to predict the folding ability of model proteins

We investigate the folding of a 125-bead heteropolymer model for proteins subject to Monte Carlo dynamics on a simple cubic lattice. Detailed study of a few sequences revealed a folding mechanism consisting of a rapid collapse followed by a slow search for a stable core that served as the transition state for folding to a near-native intermediate. Rearrangement from the intermediate to the native state slowed folding further because it required breaking native-like local structure between surface monomers so that those residues could condense onto the core. We demonstrate here the generality of this mechanism by a statistical analysis of a 200 sequence database using a method that employs a genetic algorithm to pick the sequence attributes that are most important for folding and an artificial neural network to derive the corresponding functional dependence of folding ability on the chosen sequence attributes [quantitative structure-property relationships (QSPRs)]. QSPRs that use three sequence attributes yielded substantially more accurate predictions than those that use only one. The results suggest that efficient search for the core is dependent on both the native state's overall stability and its amount of kinetically accessible, cooperative structure, whereas rearrangement from the intermediate is facilitated by destabilization of contacts between surface monomers. Implications for folding and design are discussed. Proteins 33:177–203, 1998. © 1998 Wiley-Liss, Inc.     

Analyses of genetic ancestry enable key insights for molecular ecology

Gene flow and recombination in admixed populations produce genomes that are mosaic combinations of chromosome segments inherited from different source populations, that is, chromosome segments with different genetic ancestries. The statistical problem of estimating genetic ancestry from DNA sequence data has been widely studied, and analyses of genetic ancestry have facilitated research in molecular ecology and ecological genetics. In this review, we describe and compare different model‐based statistical methods used to infer genetic ancestry. We describe the conceptual and mathematical structure of these models and highlight some of their key differences and shared features. We then discuss recent empirical studies that use estimates of genetic ancestry to analyse population histories, the nature and genetic basis of species boundaries, and the genetic architecture of traits. These diverse studies demonstrate the breadth of applications that rely on genetic ancestry estimates and typify the genomics‐enabled research that is becoming increasingly common in molecular ecology. We conclude by identifying key research areas where future studies might further advance this field.     

Calculating pH-dependent free energy of proteins by using Monte Carlo protonation probabilities of ionizable residues

Protein folding, stability, and function are usually influenced by pH. And free energy plays a fundamental role in analysis of such pH-dependent properties. Electrostatics-based theoretical framework using dielectric solvent continuum model and solving Poisson-Boltzmann equation numerically has been shown to be very successful in understanding the pH-dependent properties. However, in this approach the exact computation of pH-dependent free energy becomes impractical for proteins possessing more than several tens of ionizable sites (e.g.>30), because exact evaluation of the partition function requires a summation over a vast number of possible protonation microstates. Here we present a method which computes the free energy using the average energy and the protonation probabilities of ionizable sites obtained by the well-established Monte Carlo sampling procedure. The key feature is to calculate the entropy by using the protonation probabilities. We used this method to examine a well-studied protein (lysozyme) and produced results which agree very well with the exact calculations. Applications to the optimum pH of maximal stability of proteins and protein–DNA interactions have also resulted in good agreement with experimental data. These examples recommend our method for application to the elucidation of the pH-dependent properties of proteins.     

Leveraging eDNA to expand the study of hybrid zones

Hybrid zones are important windows into ecological and evolutionary processes. Our understanding of the significance and prevalence of hybridization in nature has expanded with the generation and analysis of genome‐spanning data sets. That said, most hybridization research still has restricted temporal and spatial resolution, which limits our ability to draw broad conclusions about evolutionary and conservation related outcomes. Here, we argue that rapidly advancing environmental DNA (eDNA) methodology could be adopted for studies of hybrid zones to increase temporal sampling (contemporary and historical), refine and geographically expand sampling density, and collect data for taxa that are difficult to directly sample. Genomic data in the environment offer the potential for near real‐time biological tracking of hybrid zones, and eDNA provides broad, but as yet untapped, potential to address eco‐evolutionary questions.     

In silico investigation of factors affecting the MV imaging performance of a novel water-equivalent EPID

PurposeA Geant4 model of a novel, water-equivalent electronic portal imaging device (EPID) prototype for radiotherapy imaging and dosimetry utilising an array of plastic scintillating fibres (PSFs) has been developed. Monte Carlo (MC) simulations were performed to quantify the PSF-EPID imaging performance and to investigate design aspects affecting performance for optimisation.MethodsUsing the Geant4 model, the PSF-EPID’s imaging performance for 6 MV photon beams was quantified in terms of its modulation transfer function (MTF), noise power spectrum (NPS) and detective quantum efficiency (DQE). Model parameters, including fibre dimensions, optical cladding reflectivity and scintillation yield, were varied to investigate impact on imaging performance.ResultsThe MC-calculated DQE(0) for the reference PSF-EPID geometry employing 30 mm fibres was approximately nine times greater than values reported for commercial EPIDs. When using 10 mm long fibres, the PSF-EPID DQE(0) was still approximately three times greater than that of a commercial EPID. Increased fibre length, cladding reflectivity and scintillation yield produced the greatest decreases in NPS and increases in DQE.ConclusionsThe potential to develop an optimised next-generation water-equivalent EPID with MV imaging performance at least comparable to commercial EPIDs has been demonstrated. Factors most important for optimising prototype design include fibre length, cladding reflectivity and scintillation yield.     

Statistical analysis of occupancy rates for overdispersed populations by redistribution procedures: Application to the european corn borer EGG masses distribution

Standard statistical analyses of distributions of individuals from contingency tables are generally invalid if the individuals are not distributed independently of each other. In this paper, we discuss a method of testing hypotheses about classification category occupancy rates for overdispersed population or for population whose individuals are distributed by groups rather than lonely. These methods are based on population redistribution simulations and provide valid, exact and powerful tests in situations for which classical methods are not appropriate. Illustrations are given from the European Corn Borer eggs data.     

Energy spectrum and dose enhancement due to the depth of the Lipiodol position using flattened and unflattened beams

Aim

Lipiodol was used for stereotactic body radiotherapy combining trans arterial chemoembolization. Lipiodol used for tumour seeking in trans arterial chemoembolization remains in stereotactic body radiation therapy. In our previous study, we reported the dose enhancement effect in Lipiodol with 10× flattening-filter-free (FFF). The objective of our study was to evaluate the dose enhancement and energy spectrum of photons and electrons due to the Lipiodol depth with flattened (FF) and FFF beams.

A Monte Carlo study of the influence of focal spot size,intensity distribution,breast thickness and magnification on spatial resolution of an a-Se digital mammography system using the generalized MTF

In this study the generalized Modulation Transfer Function (GMTF) and the geometric sharpness (S_geo) were used (i) to study the effects of various focal spot sizes (0.04 mm–0.3 mm), x-ray intensity distributions (Gaussian and double Gaussian), breast thicknesses (2–7 cm) and magnifications M (1.0–2.0) on the spatial resolution of an a-Se digital mammography system, (ii) to identify suitable focal spots for magnification mammography and (iii) derive optimum magnifications. For the calculation of GMTF the required components were: focal spot MTF, obtained from theory, detector MTF, scatter MTF and scatter fraction obtained from Monte Carlo simulations. The results showed that focal spots with sizes up to 0.18 mm are suitable for magnification mammography offering a GMTF which is >50% and >20% at the respective object frequencies of 6.5 mm⁻¹ and 9 mm⁻¹. Focal spots with sizes < 0.16 mm and Gaussian. intensity distribution, or sizes ≤ 0.1 mm and double Gaussian, offer a system resolution which improves or does not deteriorate with magnification for most object frequencies. For larger focal spots, i.e. 0.16–0.18 mm for a Gaussian and 0.12–0.18 mm for a double Gaussian. intensity distribution, optimum magnifications exist which depend on the object frequency and breast thickness. System resolution (in terms of S_geo) is maximized at M = 1.8–2.0 (all breast thicknesses) for Gaussian intensity distribution, and at M = 1.4–1.6 (breast thicknesses ≤ 4 cm) and M = 1.6–1.8 (thicker breasts) for double Gaussian.     

Test of the Inverse Monte Carlo Method for the Calculation of Interatomic Potential Energies in Atomic Liquids

Abstract

The principle purpose of this paper is to demonstrate the use of the Inverse Monte Carlo technique for calculating pair interaction energies in monoatomic liquids from a given equilibrium property. This method is based on the mathematical relation between transition probability and pair potential given by the fundamental equation of the “importance sampling” Monte Carlo method. In order to have well defined conditions for the test of the Inverse Monte Carlo method a Metropolis Monte Carlo simulation of a Lennard Jones liquid is carried out to give the equilibrium pair correlation function determined by the assumed potential. Because an equilibrium configuration is prerequisite for an Inverse Monte Carlo simulation a model system is generated reproducing the pair correlation function, which has been calculated by the Metropolis Monte Carlo simulation and therefore representing the system in thermal equilibrium. This configuration is used to simulate virtual atom displacements. The resulting changes in atom distribution for each single simulation step are inserted in a set of non-linear equations defining the transition probability for the virtual change of configuration. The solution of the set of equations for pair interaction energies yields the Lennard Jones potential by which the equilibrium configuration has been determined.