Iterated birth and death process as a model of radiation cell survival

The iterated birth and death process is defined as an n-fold iteration of a stochastic process consisting of the combination of instantaneous random killing of individuals in a certain population with a given survival probability s with a Markov birth and death process describing subsequent population dynamics. A long standing problem of computing the distribution of the number of clonogenic tumor cells surviving a fractionated radiation schedule consisting of n equal doses separated by equal time intervals tau is solved within the framework of iterated birth and death processes. For any initial tumor size i, an explicit formula for the distribution of the number M of surviving clonogens at moment tau after the end of treatment is found. It is shown that if i-->infinity and s-->0 so that is(n) tends to a finite positive limit, the distribution of random variable M converges to a probability distribution, and a formula for the latter is obtained. This result generalizes the classical theorem about the Poisson limit of a sequence of binomial distributions. The exact and limiting distributions are also found for the number of surviving clonogens immediately after the nth exposure. In this case, the limiting distribution turns out to be a Poisson distribution.     

Dynamics of human T-cell lymphotropic virus I (HTLV-I) infection of CD4+ T-cells

Stilianakis and Seydel (Bull. Math. Biol., 1999) proposed an ODE model that describes the T-cell dynamics of human T-cell lymphotropic virus I (HTLV-I) infection and the development of adult T-cell leukemia (ATL). Their model consists of four components: uninfected healthy CD4+ T-cells, latently infected CD4+ T-cells, actively infected CD4+ T-cells, and ATL cells. Mathematical analysis that completely determines the global dynamics of this model has been done by Wang et al. (Math. Biosci., 2002). In this note, we first modify the parameters of the model to distinguish between contact and infectivity rates. Then we introduce a discrete time delay to the model to describe the time between emission of contagious particles by active CD4+ T-cells and infection of pure cells. Using the results in Culshaw and Ruan (Math. Biosci., 2000) in the analysis of time delay with respect to cell-free viral spread of HIV, we study the effect of time delay on the stability of the endemically infected equilibrium. Numerical simulations are presented to illustrate the results.     

A novel analytical method for evolutionary graph theory problems

Evolutionary graph theory studies the evolutionary dynamics of populations structured on graphs. A central problem is determining the probability that a small number of mutants overtake a population. Currently, Monte Carlo simulations are used for estimating such fixation probabilities on general directed graphs, since no good analytical methods exist. In this paper, we introduce a novel deterministic framework for computing fixation probabilities for strongly connected, directed, weighted evolutionary graphs under neutral drift. We show how this framework can also be used to calculate the expected number of mutants at a given time step (even if we relax the assumption that the graph is strongly connected), how it can extend to other related models (e.g. voter model), how our framework can provide non-trivial bounds for fixation probability in the case of an advantageous mutant, and how it can be used to find a non-trivial lower bound on the mean time to fixation. We provide various experimental results determining fixation probabilities and expected number of mutants on different graphs. Among these, we show that our method consistently outperforms Monte Carlo simulations in speed by several orders of magnitude. Finally we show how our approach can provide insight into synaptic competition in neurology.     

A compound poisson process for relaxing the molecular clock

The molecular clock hypothesis remains an important conceptual and analytical tool in evolutionary biology despite the repeated observation that the clock hypothesis does not perfectly explain observed DNA sequence variation. We introduce a parametric model that relaxes the molecular clock by allowing rates to vary across lineages according to a compound Poisson process. Events of substitution rate change are placed onto a phylogenetic tree according to a Poisson process. When an event of substitution rate change occurs, the current rate of substitution is modified by a gamma-distributed random variable. Parameters of the model can be estimated using Bayesian inference. We use Markov chain Monte Carlo integration to evaluate the posterior probability distribution because the posterior probability involves high dimensional integrals and summations. Specifically, we use the Metropolis-Hastings-Green algorithm with 11 different move types to evaluate the posterior distribution. We demonstrate the method by analyzing a complete mtDNA sequence data set from 23 mammals. The model presented here has several potential advantages over other models that have been proposed to relax the clock because it is parametric and does not assume that rates change only at speciation events. This model should prove useful for estimating divergence times when substitution rates vary across lineages.     

An alternative approach to modelling relapse in cancer with an application to adenocarcinoma of the prostate

This paper proposes an alternative approach to modelling relapse in cancer. In particular, the dynamic model for the tumor or biomarker will be subjected to a lower elastic boundary at which the process either will be absorbed or reflected. The likelihood of reflection then can be interpreted as the probability of relapse. This framework will be exemplified for prostatic cancer by extending the recently proposed stochastic model of Dayananda et al. [P.W.A. Dayananda, J.T. Kemper, M.M. Shvartsman, A stochastic model for prostate-specific antigen levels, Math. Biosci. 190 (2004) 113] that focussed on the dynamics of the prostate-specific antigen (PSA) biomarker. Analytical results for the conditional density function, given a non-negative lower boundary, are obtained for the extreme cases of certain cure and of certain relapse. Simulations illustrate the relevance of the relapse probability and of the normal value of the biomarker for the design of treatment strategies. The paper thus points to two additional (patient-specific) characteristics that might enter treatment design and monitoring of progress in therapy.     

How well is the probability of tumor cure after fractionated irradiation described by Poisson statistics?

The probability of tumor cure in a homogeneous population of tumors exposed to fractionated radiotherapy was modeled using numerical simulations and compared with the predictions of Poisson statistics, assuming exact knowledge of the relevant tumor parameters (clonogen number, radiosensitivity, and growth kinetics). The results show that although Poisson statistics (based on exact knowledge of all parameters) accurately describes the probability of tumor cure when no proliferation occurs during treatment, it underestimates the cure rate when proliferation does occur. In practice, however, the inaccuracy is not likely to be more than about 10%. When the tumor parameters are unknown and are estimated by fitting an empirical Poisson model to tumor-cure data from a homogeneous population of proliferative tumors, the resulting estimates of tumor growth rate and radiosensitivity accurately reflect the true values, but the estimate of initial clonogen number is biased downward. A new formula that is more accurate than Poisson statistics in predicting the probability of tumor cure when proliferation occurs during treatment is discussed.     

Mixed Poisson regression models with covariate dependent rates

This paper studies a class of Poisson mixture models that includes covariates in rates. This model contains Poisson regression and independent Poisson mixtures as special cases. Estimation methods based on the EM and quasi-Newton algorithms, properties of these estimates, a model selection procedure, residual analysis, and goodness-of-fit test are discussed. A Monte Carlo study investigates implementation and model choice issues. This methodology is used to analyze seizure frequency and Ames salmonella assay data.     

Calculation of the cumulative distribution function of the time to a small observable tumor

Multistage mathematical models of carcinogenesis (when applied to tumor incidence data) have historically assumed that the growth kinetics of cells in the malignant state are disregarded and the formation of a single malignant cell is equated with the emergence of a detectable tumor. The justification of this simplification is, from a mathematical point of view, to make the estimation of tumor incidence rates tractable. However, analytical forms are not mandatory in the estimation of tumor incidence rates. Portier et al. (1996b, Math. Biosci. 135, 129–146) have demonstrated the utility of the Kolmogorov backward equations in numerically calculating tumor incidence. By extending their results, the cumulative distribution function of the time to a small observable tumor may be numerically obtained.     

Simulation methods for protein structure fluctuations

Three numerical techniques for generating thermally accessible configurations of globular proteins are considered; these techniques are the molecular dynamics method, the Metropolis Monte Carlo method, and a modified Monte Carlo method which takes account of the forces acting on the protein atoms. The molecular dynamics method is shown to be more efficient than either of the Monte Carlo methods. Because it may be necessary to use Monte Carlo methods in certain important types of sampling problems, the behavior of these methods is examined in some detail. It is found that an acceptance ratio close to 1/6 yields optimum efficiency for the Metropolis method, in contrast to what is often assumed. This result, together with the overall inefficiency of the Monte Carlo methods, appears to arise from the anisotropic forces acting on the protein atoms due to their covalent bonding. Possible ways of improving the Monte Carlo methods are suggested.     

A biological interpretation of transient anomalous subdiffusion. II. Reaction kinetics

Reaction kinetics in a cell or cell membrane is modeled in terms of the first passage time for a random walker at a random initial position to reach an immobile target site in the presence of a hierarchy of nonreactive binding sites. Monte Carlo calculations are carried out for the triangular, square, and cubic lattices. The mean capture time is expressed as the product of three factors: the analytical expression of Montroll for the capture time in a system with a single target and no binding sites; an exact expression for the mean escape time from the set of lattice points; and a correction factor for the number of targets present. The correction factor, obtained from Monte Carlo calculations, is between one and two. Trapping may contribute significantly to noise in reaction rates. The statistical distribution of capture times is obtained from Monte Carlo calculations and shows a crossover from power-law to exponential behavior. The distribution is analyzed using probability generating functions; this analysis resolves the contributions of the different sources of randomness to the distribution of capture times. This analysis predicts the distribution function for a lattice with perfect mixing; deviations reflect imperfect mixing in an ordinary random walk.     

A Scan Statistic for Binary Outcome Based on Hypergeometric Probability Model,with an Application to Detecting Spatial Clusters of Japanese Encephalitis

As a useful tool for geographical cluster detection of events, the spatial scan statistic is widely applied in many fields and plays an increasingly important role. The classic version of the spatial scan statistic for the binary outcome is developed by Kulldorff, based on the Bernoulli or the Poisson probability model. In this paper, we apply the Hypergeometric probability model to construct the likelihood function under the null hypothesis. Compared with existing methods, the likelihood function under the null hypothesis is an alternative and indirect method to identify the potential cluster, and the test statistic is the extreme value of the likelihood function. Similar with Kulldorff’s methods, we adopt Monte Carlo test for the test of significance. Both methods are applied for detecting spatial clusters of Japanese encephalitis in Sichuan province, China, in 2009, and the detected clusters are identical. Through a simulation to independent benchmark data, it is indicated that the test statistic based on the Hypergeometric model outweighs Kulldorff’s statistics for clusters of high population density or large size; otherwise Kulldorff’s statistics are superior.     

A Monte Carlo approach to calculate dose distribution around the lineal brachytherapy sources.

Monte Carlo means the statistical methods used to solve stochastic or deterministic physical or mathematical problems. The use of this technique for solving deterministic problems is indicated whenever there is no analytical solution for the problem or when this solution is very difficult or when the use of numerical methods requires an excessive amount of CPU-time. This paper describes the application of Monte Carlo methods using a computer program, called ISODOSE, which calculates isodose curves around linear radioactive sources to be used in brachytherapy treatment-planning. Brachytherapy is a special form of cancer treatment in which the radioactive source is placed near or inside the tumor, in order to produce cancer tissue necrosis. The program is written in C language, and the results will be compared to similar data obtained from a commercially available program at Hospital do Cancer de Recife (Brazil), radiation therapy institute. The use of Monte Carlo techniques in the ISODOSE program allows for plotting isodose curves around linear sources, and it is especially more precise near the borders of the source (singularities), taking less time than it would be possible by using deterministic methods. The ISODOSE program can be used for brachytherapy planning in small clinics in Recife and adjacent areas.     

Poisson Sequential Sampling Modified Towards Maximal Safety in Adverse Event Monitoring

For applying Poisson sequential sampling in safety monitoring concerning new medical conditions (as well as in general control of safety) the one-sample Poisson sequential probability ratio test (SPRT) is modified. The modification results in an adaptive procedure with respect to the smallest detectable relative risk at any time point in the course of the monitoring process. Although this quality would imply maximal safety, the result must be seen as a modification towards maximal safety since the Wald constants of the Poisson SPRT, basic for the derivation of the proposed procedure, are approximate only. The procedure may detect high relative risks even by one observed event right after start of monitoring. The alternative hypothesis must not be specified and the chosen error probabilities α and β completely specify the sampling scheme which can be used with any value of the null hypothesis. The power function and run lengths are established by Monte Carlo runs for an exemplifying choice of α = 10% and β = 5%. Three specified SPRTs are evaluated, by Monte Carlo runs also, to compare their results with those of the modified procedure.     

Ring closure probabilities for DNA fragments by Monte Carlo simulation

The rate of ligation of DNA molecules into circular forms depends on the ring closure probability, commonly called the j-factor, which is a sensitive measure of the extent to which thermal fluctuations contribute to bending and twisting of DNA molecules in solution. We present a theoretical treatment of the cyclization equilibria of DNA that employs a special Monte Carlo method for generating large ensembles of model DNA chains. Using this method, the chain length dependence of the j-factor was calculated for molecules. in the size range 250 to 2000 base-pairs. The Monte Carlo results are compared with recent analytical theory and experimental data. We show that a value of 475 A for the persistence length of DNA, close to values measured by a number of other methods, is in excellent agreement with the cyclization results. Preliminary applications of the Monte Carlo method to the problem of systematically bent DNA molecules are presented. The calculated j-factor is shown to be very sensitive to the amount of bending in these fragments. This fact suggests that ligase closure measurements of systematically bent DNA molecules should be a useful method for studying sequence-directed bending in DNA.     

Determination of order parameters and correlation times in proteins: A comparison between Bayesian, Monte Carlo and simple graphical methods

We describe a novel approach to deducing order parameters and correlation times in proteins using a Bayesian statistical method, and show how likelihood contours, P(,S), and confidence levels can be obtained. These results are then compared with those obtained from a simple graphical method, as well as those from Monte Carlo simulations. The Bayes approach has the advantage that it is simple and accurate. Unlike Monte Carlo methods, it gives useful contour plots of probability (also not provided by the simple graphical method), and provides likelihood/confidence information. In addition, the Bayesian approach gives results in very good agreement with those obtained from Monte Carlo simulations, and as such use of Bayesian statistical methods appears to have a promising future for studies of order and dynamics in macromolecules.     

On the number of Monte Carlo runs in comparative probabilistic LCA

Introduction

The Monte Carlo technique is widely used and recommended for including uncertainties LCA. Typically, 1000 or 10,000 runs are done, but a clear argument for that number is not available, and with the growing size of LCA databases, an excessively high number of runs may be a time-consuming thing. We therefore investigate if a large number of runs are useful, or if it might be unnecessary or even harmful.

Probability theory

We review the standard theory or probability distributions for describing stochastic variables, including the combination of different stochastic variables into a calculation. We also review the standard theory of inferential statistics for estimating a probability distribution, given a sample of values. For estimating the distribution of a function of probability distributions, two major techniques are available, analytical, applying probability theory and numerical, using Monte Carlo simulation. Because the analytical technique is often unavailable, the obvious way-out is Monte Carlo. However, we demonstrate and illustrate that it leads to overly precise conclusions on the values of estimated parameters, and to incorrect hypothesis tests.

Numerical illustration

We demonstrate the effect for two simple cases: one system in a stand-alone analysis and a comparative analysis of two alternative systems. Both cases illustrate that statistical hypotheses that should not be rejected in fact are rejected in a highly convincing way, thus pointing out a fundamental flaw.

Discussion and conclusions

Apart form the obvious recommendation to use larger samples for estimating input distributions, we suggest to restrict the number of Monte Carlo runs to a number not greater than the sample sizes used for the input parameters. As a final note, when the input parameters are not estimated using samples, but through a procedure, such as the popular pedigree approach, the Monte Carlo approach should not be used at all.

     

Rotational dynamics of surface probes in lipid vesicles

Translational and rotational diffusion of fluorescent molecules on the surface of small biological systems such as vesicles, proteins and micelles depolarize the fluorescence. A recent study has treated the case of the translational dynamics of surface probes (M.M.G. Krishna, R. Das, N. Periasamy and R. Nityananda, J. Chem. Phys., 112 (2000) 8502-8514) using Monte Carlo and theoretical methods. Here we extend the application of the methodologies to apply the case of rotational dynamics of surface probes. The corresponding fluorescence anisotropy decays were obtained using the Monte Carlo simulation methods for the two cases: surface probes undergoing rotational dynamics on a plane and on a sphere. The results were consistent with the theoretical equations which show that Monte Carlo methods can be used to simulate the surface diffusion problems. The anisotropy decay for the rotational diffusion of a molecule on a planar surface is single exponential and the residual anisotropy is zero. However, residual anisotropy is finite for the case of rotational diffusion on a sphere because of the spatial averaging of the anisotropy function. The rotational correlation time in both the cases is (4Drot)(-1) with Drot being the rotational diffusion coefficient. Rotational dynamics of a surface bound dye in a single giant liposome and in sonicated vesicles were studied and the results were explained according to the theoretical equations. A fast component of fluorescence depolarization was also observed for sonicated vesicles which was interpreted as wobbling-in-cylinder dynamics of the surface-bound dye.     

Modeling and computer simulations of tumor growth and tumor response to radiotherapy

A model of tumor growth and tumor response to radiation is introduced in which each tumor cell is taken into account individually. Each cell is assigned a set of radiobiological parameters, and the status of each cell is checked in discrete intervals. Tumor proliferation is governed by the cell cycle times of tumor cells, the growth fraction, the apoptotic capacity of the tumor, and the degree of tumor angiogenesis. The response of tumor cells to radiation is determined by the radiosensitivities and the oxygenation status. Computer simulation is performed on a 3D rigid cubic lattice, starting out from a single tumor cell. Random processes are simulated by Monte Carlo methods. Short cell cycle time, high growth fraction, and tumor angiogenesis all increase tumor proliferation rates. Accelerated time-dose patterns result in lower total doses needed for tumor control, but the extent of dose reduction depends on the kinetics and the radiosensitivities of tumor cells. Tumor angiogenesis alters fully oxygenated and hypoxic fractions within the tumor and subsequently affects the radiation response. It is demonstrated for selected radiobiological parameters that the simulation tools are suitable to quantitatively assess the total doses needed for tumor control. Using the simulation tools, it is feasible to simulate time-dependent effects during fractionated radiotherapy and to compare different time-dose patterns in terms of their tumor control.