Estimation problems associated with stochastic modeling of proliferation and differentiation of O-2A progenitor cells in vitro

Our previous research effort has resulted in a stochastic model that provides an excellent fit to our experimental data on proliferation and differentiation of oligodendrocyte type-2 astrocyte progenitor cells at the clonal level. However, methods for estimation of model parameters and their statistical properties still remain far away from complete exploration. The main technical difficulty is that no explicit analytic expression for the joint distribution of the number of progenitor cells and oligodendrocytes, and consequently for the corresponding likelihood function, is available. In the present paper, we overcome this difficulty by using computer-intensive simulation techniques for estimation of the likelihood function. Since the output of our simulation model is essentially random, stochastic optimization methods are necessary to maximize the estimated likelihood function. We use the Kiefer-Wolfowitz procedure for this purpose. Given sufficient computing resources, the proposed estimation techniques significantly extend the spectrum of problems that become approachable. In particular, these techniques can be applied to more complex branching models of multi-type cell systems with dependent evolutions of different types of cells.     

Cell inactivation by ionizing particles and the shapes of survival curves

Recent experiments concerning the survival of monolayer cells irradiated by different parts of ion Bragg peaks opened a way to a deeper mechanistic understanding of cell inactivation. A new theoretical formula for survival curves has been derived reflecting two basic phases of the given mechanism, i.e. energy transfer to a cell nucleus and subsequent biological effect (depending on the amount of imparted energy). The survival ratio for a given dose has been expressed as a function of inactivation probabilities of individual cells after different numbers of nucleus hits (a given amount of energy being transferred to a cell nucleus in each ion traversal). Having used the experimental data for V79 cells irradiated by protons, deuterons and helium ions in different parts of Bragg peaks preliminary values of these inactivation probabilities for individual cells at different LET values have been established.     

Statistical inference from single channel records: two-state Markov model with limited time resolution

Though stochastic models are widely used to describe single ion channel behaviour, statistical inference based on them has received little consideration. This paper describes techniques of statistical inference, in particular likelihood methods, suitable for Markov models incorporating limited time resolution by means of a discrete detection limit. To simplify the analysis, attention is restricted to two-state models, although the methods have more general applicability. Non-uniqueness of the mean open-time and mean closed-time estimators obtained by moment methods based on single exponential approximations to the apparent open-time and apparent closed-time distributions has been reported. The present study clarifies and extends this previous work by proving that, for such approximations, the likelihood equations as well as the moment equations (usually) have multiple solutions. Such non-uniqueness corresponds to non-identifiability of the statistical model for the apparent quantities. By contrast, higher-order approximations yield theoretically identifiable models. Likelihood-based estimation procedures are developed for both single exponential and bi-exponential approximations. The methods and results are illustrated by numerical examples based on literature and simulated data, with consideration given to empirical distributions and model control, likelihood plots, and point estimation and confidence regions.     

Production and characterization of monoclonal antibodies to rabbit lymphocyte subpopulations. I. Tissue immunofluorescence and flow cytometric analysis

Nine monoclonal antibodies to rabbit T cells and B subpopulations have been generated from three separate fusions of spleen cells from mice immunized with fractionated populations of rabbit lymphocytes. These monoclonal antibodies, as well as a previously described rabbit T cell monoclonal antibody, 9AE10, have been analyzed by immunofluorescence staining on frozen tissue sections of rabbit thymus, spleen, and appendix. This screening method permits rapid identification of the lymphocyte subdomains in each tissue which is not possible by other screening methods. Each monoclonal antibody selected has a unique tissue staining pattern. Flow cytometric analysis of these monoclonal antibodies, using indirect immunofluorescence techniques on thymocytes, splenocytes, and PBL, revealed varying percentages of positive cells and individual mean fluorescence intensities indicating different epitope densities for each antigen. These monoclonal antibodies are now being used to characterize normal lymphocyte function and the role of specific lymphocyte subpopulations in experimental disease models in the rabbit.     

A Semiparametric Joint Model for Longitudinal and Survival Data with Application to Hemodialysis Study

Summary .  In many longitudinal clinical studies, the level and progression rate of repeatedly measured biomarkers on each subject quantify the severity of the disease and that subject's susceptibility to progression of the disease. It is of scientific and clinical interest to relate such quantities to a later time-to-event clinical endpoint such as patient survival. This is usually done with a shared parameter model. In such models, the longitudinal biomarker data and the survival outcome of each subject are assumed to be conditionally independent given subject-level severity or susceptibility (also called frailty in statistical terms). In this article, we study the case where the conditional distribution of longitudinal data is modeled by a linear mixed-effect model, and the conditional distribution of the survival data is given by a Cox proportional hazard model. We allow unknown regression coefficients and time-dependent covariates in both models. The proposed estimators are maximizers of an exact correction to the joint log likelihood with the frailties eliminated as nuisance parameters, an idea that originated from correction of covariate measurement error in measurement error models. The corrected joint log likelihood is shown to be asymptotically concave and leads to consistent and asymptotically normal estimators. Unlike most published methods for joint modeling, the proposed estimation procedure does not rely on distributional assumptions of the frailties. The proposed method was studied in simulations and applied to a data set from the Hemodialysis Study.     

Structure-preserving model reduction of passive and quasi-active neurons

The spatial component of input signals often carries information crucial to a neuron’s function, but models mapping synaptic inputs to the transmembrane potential can be computationally expensive. Existing reduced models of the neuron either merge compartments, thereby sacrificing the spatial specificity of inputs, or apply model reduction techniques that sacrifice the underlying electrophysiology of the model. We use Krylov subspace projection methods to construct reduced models of passive and quasi-active neurons that preserve both the spatial specificity of inputs and the electrophysiological interpretation as an RC and RLC circuit, respectively. Each reduced model accurately computes the potential at the spike initiation zone (SIZ) given a much smaller dimension and simulation time, as we show numerically and theoretically. The structure is preserved through the similarity in the circuit representations, for which we provide circuit diagrams and mathematical expressions for the circuit elements. Furthermore, the transformation from the full to the reduced system is straightforward and depends on intrinsic properties of the dendrite. As each reduced model is accurate and has a clear electrophysiological interpretation, the reduced models can be used not only to simulate morphologically accurate neurons but also to examine computations performed in dendrites.     

Integral equation methods for computing likelihoods and their derivatives in the stochastic integrate-and-fire model

We recently introduced likelihood-based methods for fitting stochastic integrate-and-fire models to spike train data. The key component of this method involves the likelihood that the model will emit a spike at a given time t. Computing this likelihood is equivalent to computing a Markov first passage time density (the probability that the model voltage crosses threshold for the first time at time t). Here we detail an improved method for computing this likelihood, based on solving a certain integral equation. This integral equation method has several advantages over the techniques discussed in our previous work: in particular, the new method has fewer free parameters and is easily differentiable (for gradient computations). The new method is also easily adaptable for the case in which the model conductance, not just the input current, is time-varying. Finally, we describe how to incorporate large deviations approximations to very small likelihoods. Action Editor: Barry J. Richmond     

Bayesian case influence diagnostics for survival models

We propose Bayesian case influence diagnostics for complex survival models. We develop case deletion influence diagnostics for both the joint and marginal posterior distributions based on the Kullback-Leibler divergence (K-L divergence). We present a simplified expression for computing the K-L divergence between the posterior with the full data and the posterior based on single case deletion, as well as investigate its relationships to the conditional predictive ordinate. All the computations for the proposed diagnostic measures can be easily done using Markov chain Monte Carlo samples from the full data posterior distribution. We consider the Cox model with a gamma process prior on the cumulative baseline hazard. We also present a theoretical relationship between our case-deletion diagnostics and diagnostics based on Cox's partial likelihood. A simulated data example and two real data examples are given to demonstrate the methodology.     

Maximum Likelihood Estimates for Binary Data with Random Effects

The purpose of this paper is to present a procedure for obtaining approximate maximum likelihood estimates for compound binary response models. The extra binomial variation is incorporated into the model by adding random effects to the fixed effects on the probit (or logit) scale. Numerical integration techniques are used to arrive at a solution of the likelihood equations. The paper also presents an illustrating numerical example based on a large toxicological data set. The computations are carried out within the GLIM statistical package.     

Transmission of Mechanical Stresses within the Cytoskeleton of Adherent Cells: A Theoretical Analysis Based on a Multi-Component Cell Model

How environmental mechanical forces affect cellular functions is a central problem in cell biology. Theoretical models of cellular biomechanics provide relevant tools for understanding how the contributions of deformable intracellular components and specific adhesion conditions at the cell interface are integrated for determining the overall balance of mechanical forces within the cell. We investigate here the spatial distributions of intracellular stresses when adherent cells are probed by magnetic twisting cytometry. The influence of the cell nucleus stiffness on the simulated nonlinear torque-bead rotation response is analyzed by considering a finite element multi-component cell model in which the cell and its nucleus are considered as different hyperelastic materials. We additionally take into account the mechanical properties of the basal cell cortex, which can be affected by the interaction of the basal cell membrane with the extracellular substrate. In agreement with data obtained on epithelial cells, the simulated behaviour of the cell model relates the hyperelastic response observed at the entire cell scale to the distribution of stresses and strains within the nucleus and the cytoskeleton, up to cell adhesion areas. These results, which indicate how mechanical forces are transmitted at distant points through the cytoskeleton, are compared to recent data imaging the highly localized distribution of intracellular stresses.     

Survival of bowhead whales,Balaena mysticetus,estimated from 1981-1998 photoidentification data

Annual survival probability of bowhead whales, Balaena mysticetus, was estimated using both Bayesian and maximum likelihood implementations of Cormack and Jolly-Seber (JS) models for capture-recapture estimation in open populations and reduced-parameter generalizations of these models. Aerial photographs of naturally marked bowheads collected between 1981 and 1998 provided the data. The marked whales first photographed in a particular year provided the initial 'capture' and 'release' of those marked whales and photographs in subsequent years the 'recaptures'. The Cormack model, often called the Cormack-Jolly-Seber (CJS) model, and the program MARK were used to identify the model with a single survival and time-varying capture probabilities as the most appropriate for these data. When survival was constrained to be one or less, the maximum likelihood estimate computed by MARK was one, invalidating confidence interval computations based on the asymptotic standard error or profile likelihood. A Bayesian Markov chain Monte Carlo (MCMC) implementation of the model was used to produce a posterior distribution for annual survival. The corresponding reduced-parameter JS model was also fit via MCMC because it is the more appropriate of the two models for these photoidentification data. Because the CJS model ignores much of the information on capture probabilities provided by the data, its results are less precise and more sensitive to the prior distributions used than results from the JS model. With priors for annual survival and capture probabilities uniform from 0 to 1, the posterior mean for bowhead survival rate from the JS model is 0.984, and 95% of the posterior probability lies between 0.948 and 1. This high estimated survival rate is consistent with other bowhead life history data.     

Experiments on movement of DNA regions in Escherichia coli evaluated by computer simulation

During the cell cycle of Escherichia coli DNA is replicated and segregated over two prospective daughter cells. Nucleoids as a whole separate gradually in line with cell elongation, but sub-nucleoid DNA regions may behave differently, separating non-gradually. We tested the ability of three models to predict the outcome of a fluorescent in situ hybridisation (FISH) experiment. We did this by comparing computer-simulated data with experimental data. The first model predicts gradual separation in line with cell elongation. The second model predicts that origins stick together for some time after duplication before one copy jumps to the other side of the cell (non-gradual separation). The simulated data of these models are very similar, indicating that FISH is not a suitable method to distinguish between these two models. The third model predicts that origins may be anywhere within the nucleoid(s). We found that simulated data using the third model resemble the experimental data most. However, DNA regions are not randomly localised in the cell, although their localisation is fuzzy. We propose that movement of DNA regions is the result of a combination of factors. Nucleoid segregation (or the forces behind it) dictates the overall direction of movement. Other factors, of which we show that diffusion could be an important one, move DNA in other directions giving rise to non-gradual movement in individual cells and contributing to variation in intracellular position per cell length in a population of cells.     

Estimating parameters and predicting membrane voltages with conductance-based neuron models

Recent results demonstrate techniques for fully quantitative, statistical inference of the dynamics of individual neurons under the Hodgkin–Huxley framework of voltage-gated conductances. Using a variational approximation, this approach has been successfully applied to simulated data from model neurons. Here, we use this method to analyze a population of real neurons recorded in a slice preparation of the zebra finch forebrain nucleus HVC. Our results demonstrate that using only 1,500 ms of voltage recorded while injecting a complex current waveform, we can estimate the values of 12 state variables and 72 parameters in a dynamical model, such that the model accurately predicts the responses of the neuron to novel injected currents. A less complex model produced consistently worse predictions, indicating that the additional currents contribute significantly to the dynamics of these neurons. Preliminary results indicate some differences in the channel complement of the models for different classes of HVC neurons, which accords with expectations from the biology. Whereas the model for each cell is incomplete (representing only the somatic compartment, and likely to be missing classes of channels that the real neurons possess), our approach opens the possibility to investigate in modeling the plausibility of additional classes of channels the cell might possess, thus improving the models over time. These results provide an important foundational basis for building biologically realistic network models, such as the one in HVC that contributes to the process of song production and developmental vocal learning in songbirds.     

A Comparative Survey of Automated Parameter-Search Methods for Compartmental Neural Models

One of the most difficult and time-consuming aspects of building compartmental models of single neurons is assigning values to free parameters to make models match experimental data. Automated parameter-search methods potentially represent a more rapid and less labor-intensive alternative to choosing parameters manually. Here we compare the performance of four different parameter-search methods on several single-neuron models. The methods compared are conjugate-gradient descent, genetic algorithms, simulated annealing, and stochastic search. Each method has been tested on five different neuronal models ranging from simple models with between 3 and 15 parameters to a realistic pyramidal cell model with 23 parameters. The results demonstrate that genetic algorithms and simulated annealing are generally the most effective methods. Simulated annealing was overwhelmingly the most effective method for simple models with small numbers of parameters, but the genetic algorithm method was equally effective for more complex models with larger numbers of parameters. The discussion considers possible explanations for these results and makes several specific recommendations for the use of parameter searches on neuronal models.     

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%.     

Genetic mapping in the presence of genotyping errors

Genetic maps are built using the genotypes of many related individuals. Genotyping errors in these data sets can distort genetic maps, especially by inflating the distances. We have extended the traditional likelihood model used for genetic mapping to include the possibility of genotyping errors. Each individual marker is assigned an error rate, which is inferred from the data, just as the genetic distances are. We have developed a software package, called TMAP, which uses this model to find maximum-likelihood maps for phase-known pedigrees. We have tested our methods using a data set in Vitis and on simulated data and confirmed that our method dramatically reduces the inflationary effect caused by increasing the number of markers and leads to more accurate orders.     

Open capture-recapture models with heterogeneity: I. Cormack-Jolly-Seber model

In open population capture-recapture studies, it is usually assumed that similar animals (e.g., of the same sex and age group) have similar survival rates and capture probabilities. These assumptions are generally perceived to be an oversimplification, and they can lead to incorrect model selection and biased parameter estimates. Allowing for individual variability in survival and capture probabilities among apparently similar animals is now becoming possible, due to advances in closed population models and improved computing power. This article presents a flexible framework of likelihood-based models which allow for individual heterogeneity in survival and capture rates. Heterogeneity is modeled using finite mixtures, which have enough flexibility of distribution shape to accommodate a wide variety of different patterns of individual variation. The models condition on the first capture of each animal, and include as a special case the Cormack-Jolly-Seber model. Model selection is done either using Akaike's information criterion or by likelihood ratio tests, making available checks of different influences on survival rates. Bias in parameter estimates is reduced by including individual heterogeneity. Model selection and bias reduction are important in population studies and for making informed management decisions.