A stochastic approximation algorithm with Markov chain Monte-Carlo method for incomplete data estimation problems

We propose a general procedure for solving incomplete data estimation problems. The procedure can be used to find the maximum likelihood estimate or to solve estimating equations in difficult cases such as estimation with the censored or truncated regression model, the nonlinear structural measurement error model, and the random effects model. The procedure is based on the general principle of stochastic approximation and the Markov chain Monte-Carlo method. Applying the theory on adaptive algorithms, we derive conditions under which the proposed procedure converges. Simulation studies also indicate that the proposed procedure consistently converges to the maximum likelihood estimate for the structural measurement error logistic regression model.

Equilibrium distributions of microsatellite repeat length resulting from a balance between slippage events and point mutations

We describe and test a Markov chain model of microsatellite evolution that can explain the different distributions of microsatellite lengths across different organisms and repeat motifs. Two key features of this model are the dependence of mutation rates on microsatellite length and a mutation process that includes both strand slippage and point mutation events. We compute the stationary distribution of allele lengths under this model and use it to fit DNA data for di-, tri-, and tetranucleotide repeats in humans, mice, fruit flies, and yeast. The best fit results lead to slippage rate estimates that are highest in mice, followed by humans, then yeast, and then fruit flies. Within each organism, the estimates are highest in di-, then tri-, and then tetranucleotide repeats. Our estimates are consistent with experimentally determined mutation rates from other studies. The results suggest that the different length distributions among organisms and repeat motifs can be explained by a simple difference in slippage rates and that selective constraints on length need not be imposed.

Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach

A maximum likelihood estimator based on the coalescent for unequal migration rates and different subpopulation sizes is developed. The method uses a Markov chain Monte Carlo approach to investigate possible genealogies with branch lengths and with migration events. Properties of the new method are shown by using simulated data from a four-population n-island model and a source–sink population model. Our estimation method as coded in migrate is tested against genetree; both programs deliver a very similar likelihood surface. The algorithm converges to the estimates fairly quickly, even when the Markov chain is started from unfavorable parameters. The method was used to estimate gene flow in the Nile valley by using mtDNA data from three human populations.

Cortical activity flips among quasi-stationary states.

Parallel recordings of spike trains of several single cortical neurons in behaving monkeys were analyzed as a hidden Markov process. The parallel spike trains were considered as a multivariate Poisson process whose vector firing rates change with time. As a consequence of this approach, the complete recording can be segmented into a sequence of a few statistically discriminated hidden states, whose dynamics are modeled as a first-order Markov chain. The biological validity and benefits of this approach were examined in several independent ways: (i) the statistical consistency of the segmentation and its correspondence to the behavior of the animals; (ii) direct measurement of the collective flips of activity, obtained by the model; and (iii) the relation between the segmentation and the pair-wise short-term cross-correlations between the recorded spike trains. Comparison with surrogate data was also carried out for each of the above examinations to assure their significance. Our results indicated the existence of well-separated states of activity, within which the firing rates were approximately stationary. With our present data we could reliably discriminate six to eight such states. The transitions between states were fast and were associated with concomitant changes of firing rates of several neurons. Different behavioral modes and stimuli were consistently reflected by different states of neural activity. Moreover, the pair-wise correlations between neurons varied considerably between the different states, supporting the hypothesis that these distinct states were brought about by the cooperative action of many neurons.

Sequence analysis of the cis-regulatory regions of the bithorax complex of Drosophila.

The bithorax complex (BX-C) of Drosophila, one of two complexes that act as master regulators of the body plan of the fly, has now been entirely sequenced and comprises approximately 315,000 bp, only 1.4% of which codes for protein. Analysis of this sequence reveals significantly overrepresented DNA motifs of unknown, as well as known, functions in the non-protein-coding portion of the sequence. The following types of motifs in that portion are analyzed: (i) concatamers of mono-, di-, and trinucleotides; (ii) tightly clustered hexanucleotides (spaced < or = 5 bases apart); (iii) direct and reverse repeats longer than 20 bp; and (iv) a number of motifs known from biochemical studies to play a role in the regulation of the BX-C. The hexanucleotide AGATAC is remarkably overrepresented and is surmised to play a role in chromosome pairing. The positions of sites of highly overrepresented motifs are plotted for those that occur at more than five sites in the sequence, when < 0.5 case is expected. Expected values are based on a third-order Markov chain, which is the optimal order for representing the BXCALL sequence.

Tertiary structure of an amyloid immunoglobulin light chain protein: a proposed model for amyloid fibril formation.

An immunoglobulin light chain protein was isolated from the urine of an individual (BRE) with systemic amyloidosis. Complete amino acid sequence of the variable region of the light chain (VL) protein established it as a kappa I, which when compared with other kappa I amyloid associated proteins had unique residues, including Ile-34, Leu-40, and Tyr-71. To study the tertiary structure, BRE VL was expressed in Escherichia coli by using a PCR product amplified from the patient BRE's bone marrow DNA. The PCR product was ligated into pCZ11, a thermal-inducible replication vector. Recombinant BRE VL was isolated, purified to homogeneity, and crystallized by using ammonium sulfate as the precipitant. Two crystal forms were obtained. In crystal form I the BRE VL kappa domain crystallizes as a dimer with unit cell constants isomorphous to previously published kappa protein structures. Comparison with a nonamyloid VL kappa domain from patient REI, identified significant differences in position of residues in the hypervariable segments plus variations in framework region (FR) segments 40-46 (FR2) and 66-67 (FR3). In addition, positional differences can be seen along the two types of local diads, corresponding to the monomer-monomer and dimer-dimer interfaces. From the packing diagram, a model for the amyloid light chain (AL) fibril is proposed based on a pseudohexagonal spiral structure with a rise of approximately the width of two dimers per 360 degree turn. This spiral structure could be consistent with the dimensions of amyloid fibrils as determined by electron microscopy.

A general model of invariant chain association with class II major histocompatibility complex proteins.

The binding of invariant chain to major histocompatibility complex (MHC) proteins is an important step in processing of MHC class II proteins and in antigen presentation. The question of how invariant chain can bind to all MHC class II proteins is central to understanding these processes. We have employed molecular modeling to predict the structure of class II-associated invariant chain peptide (CLIP)-MHC protein complexes and to ask whether the predicted mode of association could be general across all MHC class II proteins. CLIP fits identically into the MHC class II alleles HLA-DR3, I-Ak, I-Au, and I-Ad, with a consistent pattern of hydrogen bonds, contacts, and hydrophobic burial and without bad contacts. Our model predicts the burial of CLIP residues Met-91 and Met-99 in the deep P1 and P9 anchor pockets and other detailed interactions, which we have compared with available data. The predicted pattern of I-A allele-specific effects on CLIP binding is very similar to that observed experimentally by alanine-scanning mutations of CLIP. Together, these results indicate that CLIP may bind in a single, general way across products of MHC class II alleles.

A mouse model for X-linked adrenoleukodystrophy

X-linked adrenoleukodystrophy (X-ALD) is a peroxisomal disorder with impaired β-oxidation of very long chain fatty acids (VLCFAs) and reduced function of peroxisomal very long chain fatty acyl-CoA synthetase (VLCS) that leads to severe and progressive neurological disability. The X-ALD gene, identified by positional cloning, encodes a peroxisomal membrane protein (adrenoleukodystrophy protein; ALDP) that belongs to the ATP binding cassette transporter protein superfamily. Mutational analyses and functional studies of the X-ALD gene confirm that it and not VLCS is the gene responsible for X-ALD. Its role in the β-oxidation of VLCFAs and its effect on the function of VLCS are unclear. The complex pathology of X-ALD and the extreme variability of its clinical phenotypes are also unexplained. To facilitate understanding of X-ALD pathophysiology, we developed an X-ALD mouse model by gene targeting. The X-ALD mouse exhibits reduced β-oxidation of VLCFAs, resulting in significantly elevated levels of saturated VLCFAs in total lipids from all tissues measured and in cholesterol esters from adrenal glands. Lipid cleft inclusions were observed in adrenocortical cells of X-ALD mice under the electron microscope. No neurological involvement has been detected in X-ALD mice up to 6 months. We conclude that X-ALD mice exhibit biochemical defects equivalent to those found in human X-ALD and thus provide an experimental system for testing therapeutic intervention.

From residue matching patterns to protein folding topographies: General model and bovine pancreatic trypsin inhibitor

A coarse-grained model for protein-folding dynamics is introduced based on a discretized representation of torsional modes. The model, based on the Ramachandran map of the local torsional potential surface and the class (hydrophobic/polar/neutral) of each residue, recognizes patterns of both torsional conformations and hydrophobic-polar contacts, with tolerance for imperfect patterns. It incorporates empirical rates for formation of secondary and tertiary structure. The method yields a topological representation of the evolving local torsional configuration of the folding protein, modulo the basins of the Ramachandran map. The folding process is modeled as a sequence of transitions from one contact pattern to another, as the torsional patterns evolve. We test the model by applying it to the folding process of bovine pancreatic trypsin inhibitor, obtaining a kinetic description of the transitions between the contact patterns visited by the protein along the dominant folding pathway. The kinetics and detailed balance make it possible to invert the result to obtain a coarse topographic description of the potential energy surface along the dominant folding pathway, in effect to go backward or forward between a topological representation of the chain conformation and a topographical description of the potential energy surface governing the folding process. As a result, the strong structure-seeking character of bovine pancreatic trypsin inhibitor and the principal features of its folding pathway are reproduced in a reasonably quantitative way.

Monitoring the assembly of Ig light-chain amyloid fibrils by atomic force microscopy

Aggregation of Ig light chains to form amyloid fibrils is a characteristic feature of light-chain amyloidosis, a light-chain deposition disease. A recombinant variable domain of the light chain SMA was used to form amyloid fibrils in vitro. Fibril formation was monitored by atomic force microscopy imaging. Single filaments 2.4 nm in diameter were predominant at early times; protofibrils 4.0 nm in diameter were predominant at intermediate times; type I and type II fibrils 8.0 nm and 6.0 nm in diameter, respectively, were predominant at the endpoints. The increase in number of fibrils correlated with increased binding of the fluorescent dye thioflavin T. The fibrils and protofibrils showed a braided structure, suggesting that their formation involves the winding of protofibrils and filaments, respectively. These observations support a model in which two filaments combine to form a protofibril, two protofibrils intertwine to form a type I fibril, and three filaments form a type II fibril.

Efficiency of DNA replication in the polymerase chain reaction

A detailed quantitative kinetic model for the polymerase chain reaction (PCR) is developed, which allows us to predict the probability of replication of a DNA molecule in terms of the physical parameters involved in the system. The important issue of the determination of the number of PCR cycles during which this probability can be considered to be a constant is solved within the framework of the model. New phenomena of multimodality and scaling behavior in the distribution of the number of molecules after a given number of PCR cycles are presented. The relevance of the model for quantitative PCR is discussed, and a novel quantitative PCR technique is proposed.

A model for spontaneous B-lineage lymphomas in IgHμ-HOX11 transgenic mice

HOX11, a divergent homeodomain-containing transcription factor, was isolated from the breakpoint of the nonrandom t(10;14)(q24;q11) chromosome translocation found in human T cell acute lymphoblastic leukemias. The translocation places the HOX11 coding sequence under the transcriptional control of TCR α/δ regulatory elements, resulting in ectopic expression of a normal HOX11 protein in thymocytes. To investigate the oncogenic potential of HOX11, we targeted its expression in lymphocytes of transgenic mice by placing the human cellular DNA under the transcriptional control of Ig heavy chain or LCK regulatory sequences. Only IgHμ-HOX11 mice expressing low levels of HOX11 were viable. During their second year of life, all HOX11 transgenic mice became terminally ill with more than 75% developing large cell lymphomas in the spleen, which frequently disseminated to thymus, lymph nodes, and other nonhematopoietic tissues. Lymphoma cells were predominantly clonal IgM+IgD+ mature B cells. Repopulation of severe combined immunodeficient mice with cells from hyperplastic spleens indicated that the HOX11 tumor phenotype was transplantable. Before tumor development, expression of the transgene did not result in perturbations in lymphopoiesis; however, lymphoid hyperplasia involving the splenic marginal zones was present in 20% of spleens. Our studies provide direct evidence that expression of HOX11 in lymphocytes leads to malignant transformation. These mice are a useful model system to study mechanisms involved in transformation from B-lineage hyperplasia to malignant lymphoma and for testing novel approaches to therapy. They represent a novel animal model for non-Hodgkin’s lymphoma of peripheral mature B cell origin.

A Polymer Model for the Structural Organization of Chromatin Loops and Minibands in Interphase Chromosomes

A quantitative model of interphase chromosome higher-order structure is presented based on the isochore model of the genome and results obtained in the field of copolymer research. G1 chromosomes are approximated in the model as multiblock copolymers of the 30-nm chromatin fiber, which alternately contain two types of 0.5- to 1-Mbp blocks (R and G minibands) differing in GC content and DNA-bound proteins. A G1 chromosome forms a single-chain string of loop clusters (micelles), with each loop ∼1–2 Mbp in size. The number of ∼20 loops per micelle was estimated from the dependence of geometrical versus genomic distances between two points on a G1 chromosome. The greater degree of chromatin extension in R versus G minibands and a difference in the replication time for these minibands (early S phase for R versus late S phase for G) are explained in this model as a result of the location of R minibands at micelle cores and G minibands at loop apices. The estimated number of micelles per nucleus is close to the observed number of replication clusters at the onset of S phase. A relationship between chromosomal and nuclear sizes for several types of higher eukaryotic cells (insects, plants, and mammals) is well described through the micelle structure of interphase chromosomes. For yeast cells, this relationship is described by a linear coil configuration of chromosomes.

Regulation of product chain length by isoprenyl diphosphate synthases

An analysis of the x-ray structure of homodimeric avian farnesyl diphosphate synthase (geranyltransferase, EC 2.5.1.10) coupled with information about conserved amino acids obtained from a sequence alignment of 35 isoprenyl diphosphate synthases that synthesize farnesyl (C15), geranylgeranyl (C20), and higher chain length isoprenoid diphosphates suggested that the side chains of residues corresponding to F112 and F113 in the avian enzyme were important for determining the ultimate length of the hydrocarbon chains. This hypothesis was supported by site-directed mutagenesis to transform wild-type avian farnesyl diphosphate synthase (FPS) into synthases capable of producing geranylgeranyl diphosphate (F112A), geranylfarnesyl (C25) diphosphate (F113S), and longer chain prenyl diphosphates (F112A/F113S). An x-ray analysis of the structure of the F112A/F113S mutant in the apo state and with allylic substrates bound produced the strongest evidence that these mutations caused the observed change in product specificity by directly altering the size of the binding pocket for the growing isoprenoid chain in the active site of the enzyme. The proposed binding pocket in the apo mutant structure was increased in depth by 5.8 Å as compared with that for the wild-type enzyme. Allylic diphosphates were observed in the holo structures, bound through magnesium ions to the aspartates of the first of two conserved aspartate-rich sequences (D117–D121), with the hydrocarbon tails of all the ligands growing down the hydrophobic pocket toward the mutation site. A model was constructed to show how the growth of a long chain prenyl product may proceed by creation of a hydrophobic passageway from the FPS active site to the outside surface of the enzyme.

Transition-state structure as a unifying basis in protein-folding mechanisms: Contact order, chain topology, stability, and the extended nucleus mechanism

I attempt to reconcile apparently conflicting factors and mechanisms that have been proposed to determine the rate constant for two-state folding of small proteins, on the basis of general features of the structures of transition states. Φ-Value analysis implies a transition state for folding that resembles an expanded and distorted native structure, which is built around an extended nucleus. The nucleus is composed predominantly of elements of partly or well-formed native secondary structure that are stabilized by local and long-range tertiary interactions. These long-range interactions give rise to connecting loops, frequently containing the native loops that are poorly structured. I derive an equation that relates differences in the contact order of a protein to changes in the length of linking loops, which, in turn, is directly related to the unfavorable free energy of the loops in the transition state. Kinetic data on loop extension mutants of CI2 and α-spectrin SH3 domain fit the equation qualitatively. The rate of folding depends primarily on the interactions that directly stabilize the nucleus, especially those in native-like secondary structure and those resulting from the entropy loss from the connecting loops, which vary with contact order. This partitioning of energy accounts for the success of some algorithms that predict folding rates, because they use these principles either explicitly or implicitly. The extended nucleus model thus unifies the observations of rate depending on both stability and topology.