Projecting COVID-19 cases and hospital burden in Ohio.

As the Coronavirus 2019 disease (COVID-19) started to spread rapidly in the state of Ohio, the Ecology, Epidemiology and Population Health (EEPH) program within the Infectious Diseases Institute (IDI) at The Ohio State University (OSU) took the initiative to offer epidemic modeling and decision analytics support to the Ohio Department of Health (ODH). This paper describes the methodology used by the OSU/IDI response modeling team to predict statewide cases of new infections as well as potential hospital burden in the state. The methodology has two components: (1) A Dynamical Survival Analysis (DSA)-based statistical method to perform parameter inference, statewide prediction and uncertainty quantification. (2) A geographic component that down-projects statewide predicted counts to potential hospital burden across the state. We demonstrate the overall methodology with publicly available data. A Python implementation of the methodology is also made publicly available. This manuscript was submitted as part of a theme issue on "Modelling COVID-19 and Preparedness for Future Pandemics".

Necessary and sufficient conditions for exact closures of epidemic equations on configuration model networks.

We prove that it is possible to obtain the exact closure of SIR pairwise epidemic equations on a configuration model network if and only if the degree distribution follows a Poisson, binomial, or negative binomial distribution. The proof relies on establishing the equivalence, for these specific degree distributions, between the closed pairwise model and a dynamical survival analysis (DSA) model that was previously shown to be exact. Specifically, we demonstrate that the DSA model is equivalent to the well-known edge-based Volz model. Using this result, we also provide reductions of the closed pairwise and Volz models to a single equation that involves only susceptibles. This equation has a useful statistical interpretation in terms of times to infection. We provide some numerical examples to illustrate our results.

Incorporating age and delay into models for biophysical systems.

In many biological systems, chemical reactions or changes in a physical state are assumed to occur instantaneously. For describing the dynamics of those systems, Markov models that require exponentially distributed inter-event times have been used widely. However, some biophysical processes such as gene transcription and translation are known to have a significant gap between the initiation and the completion of the processes, which renders the usual assumption of exponential distribution untenable. In this paper, we consider relaxing this assumption by incorporating age-dependent random time delays (distributed according to a given probability distribution) into the system dynamics. We do so by constructing a measure-valued Markov process on a more abstract state space, which allows us to keep track of the 'ages' of molecules participating in a chemical reaction. We study the large-volume limit of such age-structured systems. We show that, when appropriately scaled, the stochastic system can be approximated by a system of partial differential equations (PDEs) in the large-volume limit, as opposed to ordinary differential equations (ODEs) in the classical theory. We show how the limiting PDE system can be used for the purpose of further model reductions and for devising efficient simulation algorithms. In order to describe the ideas, we use a simple transcription process as a running example. We, however, note that the methods developed in this paper apply to a wide class of biophysical systems.

High dimensional model representation of log-likelihood ratio: binary classification with expression data.

BACKGROUND: Binary classification rules based on a small-sample of high-dimensional data (for instance, gene expression data) are ubiquitous in modern bioinformatics. Constructing such classifiers is challenging due to (a) the complex nature of underlying biological traits, such as gene interactions, and (b) the need for highly interpretable glass-box models. We use the theory of high dimensional model representation (HDMR) to build interpretable low dimensional approximations of the log-likelihood ratio accounting for the effects of each individual gene as well as gene-gene interactions. We propose two algorithms approximating the second order HDMR expansion, and a hypothesis test based on the HDMR formulation to identify significantly dysregulated pairwise interactions. The theory is seen as flexible and requiring only a mild set of assumptions. RESULTS: We apply our approach to gene expression data from both synthetic and real (breast and lung cancer) datasets comparing it also against several popular state-of-the-art methods. The analyses suggest the proposed algorithms can be used to obtain interpretable prediction rules with high prediction accuracies and to successfully extract significantly dysregulated gene-gene interactions from the data. They also compare favorably against their competitors across multiple synthetic data scenarios. CONCLUSION: The proposed HDMR-based approach appears to produce a reliable classifier that additionally allows one to describe how individual genes or gene-gene interactions affect classification decisions. Both real and synthetic data analyses suggest that our methods can be used to identify gene networks with dysregulated pairwise interactions, and are therefore appropriate for differential networks analysis.

Ligand-Free Estrogen Receptor α (ESR1) as Master Regulator for the Expression of CYP3A4 and Other Cytochrome P450 Enzymes in the Human Liver.

Cytochrome P450 3A4 isoform (CYP3A4) transcription is controlled by hepatic transcription factors (TFs), but how TFs dynamically interact remains uncertain. We hypothesize that several TFs form a regulatory network with nonlinear, dynamic, and hierarchical interactions. To resolve complex interactions, we have applied a computational approach for estimating Sobol's sensitivity indices (SSI) under generalized linear models to existing liver RNA expression microarray data (GSE9588) and RNA-seq data from genotype-tissue expression (GTEx), generating robust importance ranking of TF effects and interactions. The SSI-based analysis identified TFs and interacting TF pairs, triplets, and quadruplets involved in CYP3A4 expression. In addition to known CYP3A4 TFs, estrogen receptor α (ESR1) emerges as key TF with the strongest main effect and as the most frequently included TF interacting partner. Model predictions were validated using small interfering RNA (siRNA)/short hairpin RNA (shRNA) gene knockdown and clustered regularly interspaced short palindromic repeats (CRISPR)-mediated transcriptional activation of ESR1 in biliary epithelial Huh7 cells and human hepatocytes in the absence of estrogen. Moreover, ESR1 and known CYP3A4 TFs mutually regulate each other. Detectable in both male and female hepatocytes without added estrogen, the results demonstrate a role for unliganded ESR1 in CYP3A4 expression consistent with unliganded ESR1 signaling reported in other cell types. Added estrogen further enhances ESR1 effects. We propose a hierarchical regulatory network for CYP3A4 expression directed by ESR1 through self-regulation, cross regulation, and TF-TF interactions. We also demonstrate that ESR1 regulates the expression of other P450 enzymes, suggesting broad influence of ESR1 on xenobiotics metabolism in human liver. Further studies are required to understand the mechanisms underlying role of ESR1 in P450 regulation. SIGNIFICANCE STATEMENT: This study focuses on identifying key transcription factors and regulatory networks for CYP3A4, the main drug metabolizing enzymes in liver. We applied a new computational approach (Sobol's sensitivity analysis) to existing hepatic gene expression data to determine the role of transcription factors in regulating CYP3A4 expression, and used molecular genetics methods (siRNA/shRNA gene knockdown and CRISPR-mediated transcriptional activation) to test these interactions in life cells. This approach reveals a robust network of TFs, including their putative interactions and the relative impact of each interaction. We find that ESR1 serves as a key transcription factor function in regulating CYP3A4, and it appears to be acting at least in part in a ligand-free fashion.

Thymus-derived regulatory T cells contribute to tolerance to commensal microbiota.

Peripheral mechanisms preventing autoimmunity and maintaining tolerance to commensal microbiota involve CD4(+) Foxp3(+) regulatory T (Treg) cells generated in the thymus or extrathymically by induction of naive CD4(+) Foxp3(-) T cells. Previous studies suggested that the T-cell receptor repertoires of thymic Treg cells and induced Treg cells are biased towards self and non-self antigens, respectively, but their relative contribution in controlling immunopathology, such as colitis and other untoward inflammatory responses triggered by different types of antigens, remains unresolved. The intestine, and especially the colon, is a particularly suitable organ to study this question, given the variety of self-, microbiota- and food-derived antigens to which Treg cells and other T-cell populations are exposed. Intestinal environments can enhance conversion to a regulatory lineage and favour tolerogenic presentation of antigens to naive CD4(+) T cells, suggesting that intestinal homeostasis depends on microbiota-specific induced Treg cells. Here, to identify the origin and antigen-specificity of intestinal Treg cells, we performed single-cell and high-throughput sequencing of the T-cell receptor repertoires of CD4(+) Foxp3(+) and CD4(+) Foxp3(-) T cells, and analysed their reactivity against specific commensal species. We show that thymus-derived Treg cells constitute most Treg cells in all lymphoid and intestinal organs, including the colon, where their repertoire is heavily influenced by the composition of the microbiota. Our results suggest that thymic Treg cells, and not induced Treg cells, dominantly mediate tolerance to antigens produced by intestinal commensals.

Quasi-Steady-State Approximations Derived from the Stochastic Model of Enzyme Kinetics.

The paper outlines a general approach to deriving quasi-steady-state approximations (QSSAs) of the stochastic reaction networks describing the Michaelis-Menten enzyme kinetics. In particular, it explains how different sets of assumptions about chemical species abundance and reaction rates lead to the standard QSSA, the total QSSA, and the reverse QSSA. These three QSSAs have been widely studied in the literature in deterministic ordinary differential equation settings, and several sets of conditions for their validity have been proposed. With the help of the multiscaling techniques introduced in Ball et al. (Ann Appl Probab 16(4):1925-1961, 2006), Kang and Kurtz (Ann Appl Probab 23(2):529-583, 2013), it is seen that the conditions for deterministic QSSAs largely agree (with some exceptions) with the ones for stochastic QSSAs in the large-volume limits. The paper also illustrates how the stochastic QSSA approach may be extended to more complex stochastic kinetic networks like, for instance, the enzyme-substrate-inhibitor system.

The sum of standardized residuals: Goodness-of-fit test for binary response models.

We propose a new goodness-of-fit statistic for evaluating generalized linear models with binary responses on the basis of the sum of standardized residuals. We derive the asymptotic distribution of the sum of standardized residuals statistic and argue that, despite its relative simplicity, it typically outperforms many of the more sophisticated currently used goodness-of-fit statistics.

Non-linear interactions between candidate genes of myocardial infarction revealed in mRNA expression profiles.

BACKGROUND: Alterations in gene expression are key events in disease etiology and risk. Poor reproducibility in detecting differentially expressed genes across studies suggests individual genes may not be sufficiently informative for complex diseases, such as myocardial infarction (MI). Rather, dysregulation of the 'molecular network' may be critical for pathogenic processes. Such a dynamic network can be built from pairwise non-linear interactions. RESULTS: We investigate non-linear interactions represented in mRNA expression profiles that integrate genetic background and environmental factors. Using logistic regression, we test the association of individual GWAS-based candidate genes and non-linear interaction terms (between these mRNA expression levels) with MI. Based on microarray data in CATHGEN (CATHeterization in GENetics) and FHS (Framingham Heart Study), we find individual genes and pairs of mRNAs, encoded by 41 MI candidate genes, with significant interaction terms in the logistic regression model. Two pairs replicate between CATHGEN and FHS (CNNM2|GUCY1A3 and CNNM2|ZEB2). Analysis of RNAseq data from GTEx (Genotype-Tissue Expression) shows that 20 % of these disease-associated RNA pairs are co-expressed, further prioritizing significant interactions. Because edges in sparse co-expression networks formed solely by the 41 candidate genes are unlikely to represent direct physical interactions, we identify additional RNAs as links between network pairs of candidate genes. This approach reveals additional mRNAs and interaction terms significant in the context of MI, for example, the path CNNM2|ACSL5|SCARF1|GUCY1A3, characterized by the common themes of magnesium and lipid processing. CONCLUSIONS: The results of this study support a role for non-linear interactions between genes in MI and provide a basis for further study of MI systems biology. mRNA expression profiles encoded by a limited number of candidate genes yield sparse networks of MI-relevant interactions that can be expanded to include additional candidates by co-expression analysis. The non-linear interactions observed here inform our understanding of the clinical relevance of gene-gene interactions in the pathophysiology of MI, while providing a new strategy in developing clinical biomarker panels.

Gene expression profiling of brain samples from patients with Lewy body dementia.

Dementia with Lewy Bodies (DLB) is the second most common neurodegenerative disorder in the elderly. The development and progression of DLB remain unclear. In this study we used next generation sequencing to assess RNA expression profiles and cellular processes associated with DLB in the anterior cingulate cortex, a brain region affected by DLB pathology. The expression measurements were made in autopsy brain tissues from 8 DLB subjects and 10 age-matched controls using AmpliSeq technology with ion torrent sequencing. The analysis of RNA expression profiles revealed 490 differentially expressed genes, among which 367 genes were down-regulated and 123 were up-regulated. Functional enrichment analysis of genes differentially expressed in DLB indicated downregulation of genes associated with myelination, neurogenesis, and regulation of nervous system development. miRNA binding sites enriched in these mRNAs yielded a list of candidate miRNAs participating in DLB pathophysiology. Our study provides a comprehensive picture of gene expression landscape in DLB, identifying key cellular processes associated with DLB pathology.

Double asymptotics for the chi-square statistic.

Consider distributional limit of the Pearson chi-square statistic when the number of classes mn increases with the sample size n and [Formula: see text]. Under mild moment conditions, the limit is Gaussian for λ = ∞, Poisson for finite λ > 0, and degenerate for λ = 0.

Neurodegeneration-associated instability of ribosomal DNA.

Homologous recombination (HR)-mediated instability of the repetitively organized ribosomal DNA (rDNA) has been proposed as a mediator of cell senescence in yeast triggering the DNA damage response. High individual variability in the content of human rDNA suggests that this genomic region remained relatively unstable throughout evolution. Therefore, quantitative real-time polymerase chain reaction was used to determine the genomic content of rDNA in post mortem samples of parietal cortex from 14 young and 9 elderly individuals with no diagnosis of a chronic neurodegenerative/neurological disease. In addition, rDNA content in that brain region was compared between 10 age-matched control individuals and 10 patients with dementia with Lewy bodies (DLB) which involves neurodegeneration of the cerebral cortex. Probing rRNA-coding regions of rDNA revealed no effects of aging on the rDNA content. Elevated rDNA content was observed in DLB. Conversely, in the DLB pathology-free cerebellum, lower genomic content of rDNA was present in the DLB group. In the parietal cortex, such a DLB-associated instability of rDNA was not accompanied by any major changes of cytosine-phosphate-guanine methylation of the rDNA promoter. As increased cerebro-cortical rDNA content was previously reported in Alzheimer's disease, neurodegeneration appears to be associated with instability of rDNA. The hypothetical origins and consequences of this phenomenon are discussed including possibilities that the DNA damage-induced recombination destabilizes rDNA and that differential content of rDNA affects heterochromatin formation, gene expression and/or DNA damage response. This article is part of a Special Issue entitled: Role of the Nucleolus in Human Disease.

Analyzing allele specific RNA expression using mixture models.

BACKGROUND: Measuring allele-specific RNA expression provides valuable insights into cis-acting genetic and epigenetic regulation of gene expression. Widespread adoption of high-throughput sequencing technologies for studying RNA expression (RNA-Seq) permits measurement of allelic RNA expression imbalance (AEI) at heterozygous single nucleotide polymorphisms (SNPs) across the entire transcriptome, and this approach has become especially popular with the emergence of large databases, such as GTEx. However, the existing binomial-type methods used to model allelic expression from RNA-seq assume a strong negative correlation between reference and variant allele reads, which may not be reasonable biologically. RESULTS: Here we propose a new strategy for AEI analysis using RNA-seq data. Under the null hypothesis of no AEI, a group of SNPs (possibly across multiple genes) is considered comparable if their respective total sums of the allelic reads are of similar magnitude. Within each group of "comparable" SNPs, we identify SNPs with AEI signal by fitting a mixture of folded Skellam distributions to the absolute values of read differences. By applying this methodology to RNA-Seq data from human autopsy brain tissues, we identified numerous instances of moderate to strong imbalanced allelic RNA expression at heterozygous SNPs. Findings with SLC1A3 mRNA exhibiting known expression differences are discussed as examples. CONCLUSION: The folded Skellam mixture model searches for SNPs with significant difference between reference and variant allele reads (adjusted for different library sizes), using information from a group of "comparable" SNPs across multiple genes. This model is particularly suitable for performing AEI analysis on genes with few heterozygous SNPs available from RNA-seq, and it can fit over-dispersed read counts without specifying the direction of the correlation between reference and variant alleles.

RNA sequencing of transcriptomes in human brain regions: protein-coding and non-coding RNAs, isoforms and alleles.

BACKGROUND: We used RNA sequencing to analyze transcript profiles of ten autopsy brain regions from ten subjects. RNA sequencing techniques were designed to detect both coding and non-coding RNA, splice isoform composition, and allelic expression. Brain regions were selected from five subjects with a documented history of smoking and five non-smokers. Paired-end RNA sequencing was performed on SOLiD instruments to a depth of >40 million reads, using linearly amplified, ribosomally depleted RNA. Sequencing libraries were prepared with both poly-dT and random hexamer primers to detect all RNA classes, including long non-coding (lncRNA), intronic and intergenic transcripts, and transcripts lacking poly-A tails, providing additional data not previously available. The study was designed to generate a database of the complete transcriptomes in brain region for gene network analyses and discovery of regulatory variants. RESULTS: Of 20,318 protein coding and 18,080 lncRNA genes annotated from GENCODE and lncipedia, 12 thousand protein coding and 2 thousand lncRNA transcripts were detectable at a conservative threshold. Of the aligned reads, 52 % were exonic, 34 % intronic and 14 % intergenic. A majority of protein coding genes (65 %) was expressed in all regions, whereas ncRNAs displayed a more restricted distribution. Profiles of RNA isoforms varied across brain regions and subjects at multiple gene loci, with neurexin 3 (NRXN3) a prominent example. Allelic RNA ratios deviating from unity were identified in > 400 genes, detectable in both protein-coding and non-coding genes, indicating the presence of cis-acting regulatory variants. Mathematical modeling was used to identify RNAs stably expressed in all brain regions (serving as potential markers for normalizing expression levels), linked to basic cellular functions. An initial analysis of differential expression analysis between smokers and nonsmokers implicated a number of genes, several previously associated with nicotine exposure. CONCLUSIONS: RNA sequencing identifies distinct and consistent differences in gene expression between brain regions, with non-coding RNA displaying greater diversity between brain regions than mRNAs. Numerous RNAs exhibit robust allele selective expression, proving a means for discovery of cis-acting regulatory factors with potential clinical relevance.

Association Between Prediagnostic Serum 25-Hydroxyvitamin D Concentration and Glioma.

There are no previous studies of the association between prediagnostic serum vitamin D concentration and glioma. Vitamin D has immunosuppressive properties; as does glioma. It was, therefore, our hypothesis that elevated vitamin D concentration would increase glioma risk. We conducted a nested case-control study using specimens from the Janus Serum Bank cohort in Norway. Blood donors who were subsequently diagnosed with glioma (n = 592), between 1974 and 2007, were matched to donors without glioma (n = 1112) on date and age at blood collection and sex. We measured 25-hydroxyvitamin D [25(OH)D], an indicator of vitamin D availability, using liquid chromatography coupled with mass spectrometry. Seasonally adjusted odds ratios (ORs) and 95% confidence intervals (95% CIs) were estimated for each control quintile of 25(OH)D using conditional logistic regression. Among men diagnosed with high grade glioma >56, we found a negative trend (P = .04). Men diagnosed ≤ 56 showed a borderline positive trend (P = .08). High levels (>66 nmol/L) of 25(OH)D in men >56 were inversely related to high grade glioma from ≥2 yr before diagnosis (OR = 0.59; 95% CI = 0.38, 0.91) to ≥15 yr before diagnosis (OR = 0.61; 95% CI = 0.38,0.96). Our findings are consistent long before glioma diagnosis and are therefore unlikely to reflect preclinical disease.

A mathematical model of pre-diagnostic glioma growth.

Due to their location, the malignant gliomas of the brain in humans are very difficult to treat in advanced stages. Blood-based biomarkers for glioma are needed for more accurate evaluation of treatment response as well as early diagnosis. However, biomarker research in primary brain tumors is challenging given their relative rarity and genetic diversity. It is further complicated by variations in the permeability of the blood brain barrier that affects the amount of marker released into the bloodstream. Inspired by recent temporal data indicating a possible decrease in serum glucose levels in patients with gliomas yet to be diagnosed, we present an ordinary differential equation model to capture early stage glioma growth. The model contains glioma-glucose-immune interactions and poses a potential mechanism by which this glucose drop can be explained. We present numerical simulations, parameter sensitivity analysis, linear stability analysis and a numerical experiment whereby we show how a dormant glioma can become malignant.

Missing heritability of common diseases and treatments outside the protein-coding exome.

Genetic factors strongly influence risk of common human diseases and treatment outcomes but the causative variants remain largely unknown; this gap has been called the 'missing heritability'. We propose several hypotheses that in combination have the potential to narrow the gap. First, given a multi-stage path from wellness to disease, we propose that common variants under positive evolutionary selection represent normal variation and gate the transition between wellness and an 'off-well' state, revealing adaptations to changing environmental conditions. In contrast, genome-wide association studies (GWAS) focus on deleterious variants conveying disease risk, accelerating the path from off-well to illness and finally specific diseases, while common 'normal' variants remain hidden in the noise. Second, epistasis (dynamic gene-gene interactions) likely assumes a central role in adaptations and evolution; yet, GWAS analyses currently are poorly designed to reveal epistasis. As gene regulation is germane to adaptation, we propose that epistasis among common normal regulatory variants, or between common variants and less frequent deleterious variants, can have strong protective or deleterious phenotypic effects. These gene-gene interactions can be highly sensitive to environmental stimuli and could account for large differences in drug response between individuals. Residing largely outside the protein-coding exome, common regulatory variants affect either transcription of coding and non-coding RNAs (regulatory SNPs, or rSNPs) or RNA functions and processing (structural RNA SNPs, or srSNPs). Third, with the vast majority of causative variants yet to be discovered, GWAS rely on surrogate markers, a confounding factor aggravated by the presence of more than one causative variant per gene and by epistasis. We propose that the confluence of these factors may be responsible to large extent for the observed heritability gap.

Dynamic SARS-CoV-2 surveillance model combining seroprevalence and wastewater concentrations for post-vaccine disease burden estimates.

BACKGROUND: Despite wide scale assessments, it remains unclear how large-scale severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) vaccination affected the wastewater concentration of the virus or the overall disease burden as measured by hospitalization rates. METHODS: We used weekly SARS-CoV-2 wastewater concentration with a stratified random sampling of seroprevalence, and linked vaccination and hospitalization data, from April 2021-August 2021 in Jefferson County, Kentucky (USA). Our susceptible ( S ), vaccinated ( V ), variant-specific infected ( I 1 and I 2 ), recovered ( R ), and seropositive ( T ) model ( S V I 2 R T ) tracked prevalence longitudinally. This was related to wastewater concentration. RESULTS: Here we show the 64% county vaccination rate translate into about a 61% decrease in SARS-CoV-2 incidence. The estimated effect of SARS-CoV-2 Delta variant emergence is a 24-fold increase of infection counts, which correspond to an over 9-fold increase in wastewater concentration. Hospitalization burden and wastewater concentration have the strongest correlation (r = 0.95) at 1 week lag. CONCLUSIONS: Our study underscores the importance of continuing environmental surveillance post-vaccine and provides a proof-of-concept for environmental epidemiology monitoring of infectious disease for future pandemic preparedness.

It is unclear how large-scale COVID-19 vaccination impacts wastewater concentration or overall disease burden. Here, we developed a mathematical surveillance model that allows estimation of overall vaccine impact based on the amount of SARS-CoV-2 in wastewater, seroprevalence and the number of cases admitted to hospitals between April 2021­August 2021 in Jefferson County, Kentucky USA. We found that a 64% vaccination coverage correlated to a 61% decrease in COVID-19 cases. The emergence of the SARS-CoV-2 Delta variant during the time of the surveillance directly correlated with a sharp increase in infection incidence as well as viral counts in wastewater. The hospitalization burden was closely reflected by the viral count found in the wastewater, indicating that post-vaccine environmental surveillance can be an effective method of estimating changing disease prevalence in future pandemics.

Inference for discretely observed stochastic kinetic networks with applications to epidemic modeling.

We present a new method for Bayesian Markov Chain Monte Carlo-based inference in certain types of stochastic models, suitable for modeling noisy epidemic data. We apply the so-called uniformization representation of a Markov process, in order to efficiently generate appropriate conditional distributions in the Gibbs sampler algorithm. The approach is shown to work well in various data-poor settings, that is, when only partial information about the epidemic process is available, as illustrated on the synthetic data from SIR-type epidemics and the Center for Disease Control and Prevention data from the onset of the H1N1 pandemic in the United States.

Bayesian multivariate Poisson abundance models for T-cell receptor data.

A major feature of an adaptive immune system is its ability to generate B- and T-cell clones capable of recognizing and neutralizing specific antigens. These clones recognize antigens with the help of the surface molecules, called antigen receptors, acquired individually during the clonal development process. In order to ensure a response to a broad range of antigens, the number of different receptor molecules is extremely large, resulting in a huge clonal diversity of both B- and T-cell receptor populations and making their experimental comparisons statistically challenging. To facilitate such comparisons, we propose a flexible parametric model of multivariate count data and illustrate its use in a simultaneous analysis of multiple antigen receptor populations derived from mammalian T-cells. The model relies on a representation of the observed receptor counts as a multivariate Poisson abundance mixture (m PAM). A Bayesian parameter fitting procedure is proposed, based on the complete posterior likelihood, rather than the conditional one used typically in similar settings. The new procedure is shown to be considerably more efficient than its conditional counterpart (as measured by the Fisher information) in the regions of m PAM parameter space relevant to model T-cell data.