Multi-level Monte Carlo for continuous time Markov chains, with applications in biochemical kinetics

We show how to extend a recently proposed multi-level Monte Carlo approach to the continuous time Markov chain setting, thereby greatly lowering the computational complexity needed to compute expected values of functions of the state of the system to a specified accuracy. The extension is non-trivial, exploiting a coupling of the requisite processes that is easy to simulate while providing a small variance for the estimator. Further, and in a stark departure from other implementations of multi-level Monte Carlo, we show how to produce an unbiased estimator that is significantly less computationally expensive than the usual unbiased estimator arising from exact algorithms in conjunction with crude Monte Carlo. We thereby dramatically improve, in a quantifiable manner, the basic computational complexity of current approaches that have many names and variants across the scientific literature, including the Bortz-Kalos-Lebowitz algorithm, discrete event simulation, dynamic Monte Carlo, kinetic Monte Carlo, the n-fold way, the next reaction method,the residence-time algorithm, the stochastic simulation algorithm, Gillespie's algorithm, and tau-leaping. The new algorithm applies generically, but we also give an example where the coupling idea alone, even without a multi-level discretization, can be used to improve efficiency by exploiting system structure. Stochastically modeled chemical reaction networks provide a very important application for this work. Hence, we use this context for our notation, terminology, natural scalings, and computational examples.Comment: Improved description of the constants in statement of Theorem

Bifurcation analysis informs Bayesian inference in the Hes1 feedback loop

<p>Background Ordinary differential equations (ODEs) are an important tool for describing the dynamics of biological systems. However, for ODE models to be useful, their parameters must first be calibrated. Parameter estimation, that is, finding parameter values given experimental data, is an inference problem that can be treated systematically through a Bayesian framework.</p> <p>A Markov chain Monte Carlo approach can then be used to sample from the appropriate posterior probability distributions, provided that suitable prior distributions can be found for the unknown parameter values. Choosing these priors is therefore a vital first step in the inference process. We study here a negative feedback loop in gene regulation where an ODE incorporating a time delay has been proposed as a realistic model and where experimental data is available. Our aim is to show that a priori mathematical analysis can be exploited in the choice of priors.</p> <p>Results By focussing on the onset of oscillatory behaviour through a Hopf Bifurcation, we derive a range of analytical expressions and constraints that link the model parameters to the observed dynamics of the system. Computational tests on both simulated and experimental data emphasise the usefulness of this analysis.</p> <p>Conclusion Mathematical analysis not only gives insights into the possible dynamical behaviour of gene expression models, but can also be used to inform the choice of priors when parameters are inferred from experimental data in a Bayesian setting.</p&gt

THE CHRONOLOGY AND STATUS OF NON NOK THA, NORTHEAST THAILAND

  Excavations at Non Nok Tha, in Northeast Thailand in 1965-1968 revealed for the first time in Southeast Asia, a stratigraphic transition from the Neolithic into the Bronze Age. Based on conventional charcoal radiocarbon determinations, early reports identified fourth millennium bronze casting. The proposed length of the prehistoric sequence, and the division of the Neolithic to Bronze age mortuary sequence into at least 11 phases, has stimulated a series of social interpretations all of which have in common, a social order based on ascriptive ranking into at least two groups which saw increased hierarchical divisions emerge with the initial Bronze Age. This paper presents the results of a new dating initiative, based on the ultrafiltration of human bones. The results indicate that the initial Neolithic occupation took place during the 14th century BC. The earliest Bronze Age has been placed in the 10th centuries BC. These dates are virtually identical with those obtained for the sites of Ban Chiang and Ban Non Wat. Compared with the elite early Bronze Age graves of Ban Non Wat, Non Nok Tha burials display little evidence for significant divisions in society

Radiocarbon dates from the Oxford AMS system: archaeometry datelist 35

This is the 35th list of AMS radiocarbon determinations measured at the Oxford Radiocarbon Accelerator Unit (ORAU). Amongst some of the sites included here are the latest series of determinations from the key sites of Abydos, El Mirón, Ban Chiang, Grotte de Pigeons (Taforalt), Alepotrypa and Oberkassel, as well as others dating to the Palaeolithic, Mesolithic and later periods. Comments on the significance of the results are provided by the submitters of the material

Hydropyrolysis: implications for radiocarbon pre-treatment and characterization of Black Carbon

Charcoal is the result of natural and anthropogenic burning events, when biomass is exposed to elevated temperatures under conditions of restricted oxygen. This process produces a range of materials, collectively known as pyrogenic carbon, the most inert fraction of which is known as Black Carbon (BC). BC degrades extremely slowly, and is resistant to diagenetic alteration involving the addition of exogenous carbon making it a useful target substance for radiocarbon dating particularly of more ancient samples, where contamination issues are critical. We present results of tests using a new method for the quantification and isolation of BC, known as hydropyrolysis (hypy). Results show controlled reductive removal of non-BC organic components in charcoal samples, including lignocellulosic and humic material. The process is reproducible and rapid, making hypy a promising new approach not only for isolation of purified BC for 14C measurement but also in quantification of different labile and resistant sample C fractions

Non-backtracking Walk Centrality for Directed Networks

The theory of zeta functions provides an expression for the generating function of nonbacktracking walk counts on a directed network. We show how this expression can be used to produce a centrality measure that eliminates backtracking walks at no cost. We also show that the radius of convergence of the generating function is related to the spectrum of a three-by-three block matrix involving the original adjacency matrix. This gives a means to choose appropriate values of the attenuation parameter. We find that three important additional benefits arise when we use this technique to eliminate traversals around the network that are unlikely to be of relevance. First, we obtain a larger range of choices for the attenuation parameter. Second, a natural approach for determining a suitable parameter range is invariant under the removal of certain types of nodes, we can gain computational efficiencies through reducing the dimension of the resulting eigenvalue problem. Third, the dimension of the linear system defining the centrality measures may be reduced in the same manner. We show that the new centrality measure may be interpreted as standard Katz on a modified network, where self loops are added, and where nonreciprocal edges are augmented with negative weights. We also give a multilayer interpretation, where negatively weighted walks between layers compensate for backtracking walks on the only non-empty layer. Studying the limit as the attenuation parameter approaches its upper bound also allows us to propose a generalization of eigenvector-based nonbacktracking centrality measure to this directed network setting. In this context, we find that the two-by-two block matrix arising in previous studies focused on undirected networks must be extended to a new three-by-three block structure to allow for directed edges. We illustrate the centrality measure on a synthetic network, where it is shown to eliminate a localization effect present in standard Katz centrality. Finally, we give results for real networks

Deep learnability: using neural networks to quantify language similarity and learnability

Learning a second language (L2) usually progresses faster if a learner's L2 is similar to their first language (L1). Yet global similarity between languages is difficult to quantify, obscuring its precise effect on learnability. Further, the combinatorial explosion of possible L1 and L2 language pairs, combined with the difficulty of controlling for idiosyncratic differences across language pairs and language learners, limits the generalisability of the experimental approach. In this study, we present a different approach, employing artificial languages and artificial learners. We built a set of five artificial languages whose underlying grammars and vocabulary were manipulated to ensure a known degree of similarity between each pair of languages. We next built a series of neural network models for each language, and sequentially trained them on pairs of languages. These models thus represented L1 speakers learning L2s. By observing the change in activity of the cells between the L1-speaker model and the L2-learner model, we estimated how much change was needed for the model to learn the new language. We then compared the change for each L1/L2 bilingual model to the underlying similarity across each language pair. The results showed that this approach can not only recover the facilitative effect of similarity on L2 acquisition, but can also offer new insights into the differential effects across different domains of similarity. These findings serve as a proof of concept for a generalisable approach that can be applied to natural languages

Discrete Razumikhin-type technique and stability of the Euler-Maruyama method to stochastic functional differential equations

A discrete stochastic Razumikhin-type theorem is established to investigate whether the Euler--Maruyama (EM) scheme can reproduce the moment exponential stability of exact solutions of stochastic functional differential equations (SFDEs). In addition, the Chebyshev inequality and the Borel-Cantelli lemma are applied to show the almost sure stability of the EM approximate solutions of SFDEs. To show our idea clearly, these results are used to discuss stability of numerical solutions of two classes of special SFDEs, including stochastic delay differential equations (SDDEs) with variable delay and stochastically perturbed equations

An Ensemble Kalman-Particle Predictor-Corrector Filter for Non-Gaussian Data Assimilation

An Ensemble Kalman Filter (EnKF, the predictor) is used make a large change in the state, followed by a Particle Filer (PF, the corrector) which assigns importance weights to describe non-Gaussian distribution. The weights are obtained by nonparametric density estimation. It is demonstrated on several numerical examples that the new predictor-corrector filter combines the advantages of the EnKF and the PF and that it is suitable for high dimensional states which are discretizations of solutions of partial differential equations.Comment: ICCS 2009, to appear; 9 pages; minor edit

Dynamic DNA and human disease: mathematical modelling and statistical inference for myotonic dystrophy type 1 and Huntington disease

Several human genetic diseases, including myotonic dystrophy type 1 (DM1) and Huntington disease (HD), are associated with inheriting an abnormally large unstable DNA simple sequence tandem repeat. These sequences mutate, by changing the number of repeats, many times during the lifetime of those affected, with a bias towards expansion. High repeat numbers are associated with early onset and disease severity. The presence of somatic instability compromises attempts to measure intergenerational repeat dynamics and infer genotype-phenotype relationships. Modelling the progression of repeat length throughout the lifetime of individuals has potential for improving prognostic information as well as providing a deeper understanding of the underlying biological process. Dr Fernando Morales, Dr Anneli Cooper and others from the Monckton lab have characterised more than 25,000 de novo somatic mutations from a large cohort of DM1 patients using single-molecule polymerase chain reaction (SM-PCR). This rich dataset enables us to fully quantify levels of somatic instability across a representative DM1 population for the first time. We establish the relationship between inherited or progenitor allele length, age at sampling and levels of somatic instability using linear regression analysis. We show that the estimated progenitor allele length genotype is significantly better than modal repeat length (the current clinical standard) at predicting age of onset and this novel genotype is the major modifier of the age of onset phenotype. Further we show that somatic variation (adjusted for estimated progenitor allele length and age at sampling) is also a modifier of the age of onset phenotype. Several families form the large cohort, and we find that the level of somatic instability is highly heritable, implying a role for individual-specific trans-acting genetic modifiers. We develop new mathematical models, the main focus of this thesis, by modifying a previously proposed stochastic birth process to incorporate possible contraction. A Bayesian likelihood approach is used as the basis for inference and parameter estimation. We use model comparison analysis to reveal, for the first time, that the expansion bias observed in the distributions of repeat lengths is likely to be the cumulative effect of many expansion and contraction events. We predict that mutation events can occur as frequently as every other day, which matches the timing of regular cell activities such as DNA repair and transcription, but not DNA replication. Mutation rates estimated under the models described above are lower than expected among individuals with inherited repeat lengths less than 100 CTGs, suggesting that these rates may be suppressed at the lower end of the disease causing range. We propose that a length-specific effect may be operating within this range and test this hypothesis by introducing such an effect into the model. To calibrate this extended model, we use blood DNA data from DM1 individuals with small alleles (inherited repeat lengths less than 100 CTGs) and buccal DNA from HD individuals who almost always have inherited repeat lengths less than 100 CAGs. These datasets comprise single DNA molecules sized using SM-PCR. We find statistical support for a general length-specific effect which suppresses mutational rates among the smaller alleles and gives rise to a distinctive pattern in the repeat length distributions. In a novel application of this new model, fitted to a large cohort of DM1 individuals, we also show that this distinctive pattern may help identify individuals whose effective repeat length, with regards to somatic instability, is less than their actual repeat length. A plausible explanation for this distinction is that the expanded repeat tract is compromised by interruptions or other unusual features. For these individuals, we estimate the effective repeat length of their expanded repeat tracts and contribute to the on-going discussion about the effect of interruptions on phenotype. The interpretation of the levels of somatic instability in many of the affected tissues in the triplet repeat diseases is hindered by complex cell compositions. We extend our model to two cell populations whose repeat lengths have different rates of mutation (fast and slow). Swami et al. have recently characterised repeat length distributions in end stage HD brain. Applying our model, we infer for each frontal cortex HD dataset the likely relative weight of these cell populations and their corresponding contribution towards somatic variation. By comparison with data from laser captured single cells we conclude that the neuronal repeat lengths most likely mutate at a higher rate than glial repeat lengths, explaining the characteristic skewed distributions observed in mixed cell tissue from the brain. We confirm that individual-specific mutation rates in neurons are, in addition to the inherited repeat length, a modifier of age of onset. Our results support a model of disease progression where individuals with the same inherited repeat length may reach age of onset, as much as 30 years earlier, because of greater somatic expansions underpinned by higher mutational rates. Therapies aimed at reducing somatic expansions would therefore have considerable benefits with regard to extending the age of onset. Currently clinical diagnosis of DM1 is based on a measure of repeat length from blood cells, but variance in modal length only accounts for between 20 - 40% of the variance in age of onset and, therefore, is not a an accurate predictive tool. We show that in principle progenitor allele length improves the inverse correlation with age of onset over the traditional model length measure. We make use of second blood samples that are now available from 40 DM1 individuals. We show that inherited repeat length and the mutation rates underlying repeat length instability in blood, inferred from samples at two time points rather than one, are better predictors of age of onset than the traditional modal length measure. Our results are a step towards providing better prognostic information for DM1 individuals and their families. They should also lead to better predictions for drug/therapy response, which is emerging as key to successful clinical trials. Microsatellites are another type of tandem repeat found in the genome with high levels of intergenerational and somatic mutation. Differences between individuals make microsatellites very useful biomarkers and they have many applications in forensics and medicine. As well as a general application to other expanded repeat diseases, the mathematical models developed here could be used to better understand instability at other mutational hotspots such as microsatellites