Stochastic journeys of cell progenies through compartments and the role of self-renewal, symmetric and asymmetric division.

Division and differentiation events by which cell populations with specific functions are generated often take place as part of a developmental programme, which can be represented by a sequence of compartments. A compartment is the set of cells with common characteristics; sharing, for instance, a spatial location or a phenotype. Differentiation events are transitions from one compartment to the next. Cells may also die or divide. We consider three different types of division events: (i) where both daughter cells inherit the mother's phenotype (self-renewal), (ii) where only one of the daughters changes phenotype (asymmetric division), and (iii) where both daughters change phenotype (symmetric division). The self-renewal probability in each compartment determines whether the progeny of a single cell, moving through the sequence of compartments, is finite or grows without bound. We analyse the progeny stochastic dynamics with probability generating functions. In the case of self-renewal, by following one of the daughters after any division event, we may construct lifelines containing only one cell at any time. We analyse the number of divisions along such lines, and the compartment where lines terminate with a death event. Analysis and numerical simulations are applied to a five-compartment model of the gradual differentiation of hematopoietic stem cells and to a model of thymocyte development: from pre-double positive to single positive (SP) cells with a bifurcation to either SP4 or SP8 in the last compartment of the sequence.

A story of viral co-infection, co-transmission and co-feeding in ticks: how to compute an invasion reproduction number.

With a single circulating vector-borne virus, the basic reproduction number incorporates contributions from tick-to-tick (co-feeding), tick-to-host and host-to-tick transmission routes. With two different circulating vector-borne viral strains, resident and invasive, and under the assumption that co-feeding is the only transmission route in a tick population, the invasion reproduction number depends on whether the model system of ordinary differential equations possesses the property of neutrality. We show that a simple model, with two populations of ticks infected with one strain, resident or invasive, and one population of co-infected ticks, does not have Alizon's neutrality property. We present model alternatives that are capable of representing the invasion potential of a novel strain by including populations of ticks dually infected with the same strain. The invasion reproduction number is analysed with the next-generation method and via numerical simulations.

Quantifying in vitro B. anthracis growth and PA production and decay: a mathematical modelling approach.

Protective antigen (PA) is a protein produced by Bacillus anthracis. It forms part of the anthrax toxin and is a key immunogen in US and UK anthrax vaccines. In this study, we have conducted experiments to quantify PA in the supernatants of cultures of B. anthracis Sterne strain, which is the strain used in the manufacture of the UK anthrax vaccine. Then, for the first time, we quantify PA production and degradation via mathematical modelling and Bayesian statistical techniques, making use of this new experimental data as well as two other independent published data sets. We propose a single mathematical model, in terms of delay differential equations (DDEs), which can explain the in vitro dynamics of all three data sets. Since we did not heat activate the B. anthracis spores prior to inoculation, germination occurred much slower in our experiments, allowing us to calibrate two additional parameters with respect to the other data sets. Our model is able to distinguish between natural PA decay and that triggered by bacteria via proteases. There is promising consistency between the different independent data sets for most of the parameter estimates. The quantitative characterisation of B. anthracis PA production and degradation obtained here will contribute towards the ambition to include a realistic description of toxin dynamics, the host immune response, and anti-toxin treatments in future mechanistic models of anthrax infection.

The reproduction number and its probability distribution for stochastic viral dynamics.

We consider stochastic models of individual infected cells. The reproduction number, R, is understood as a random variable representing the number of new cells infected by one initial infected cell in an otherwise susceptible (target cell) population. Variability in R results partly from heterogeneity in the viral burst size (the number of viral progeny generated from an infected cell during its lifetime), which depends on the distribution of cellular lifetimes and on the mechanism of virion release. We analyse viral dynamics models with an eclipse phase: the period of time after a cell is infected but before it is capable of releasing virions. The duration of the eclipse, or the subsequent infectious, phase is non-exponential, but composed of stages. We derive the probability distribution of the reproduction number for these viral dynamics models, and show it is a negative binomial distribution in the case of constant viral release from infectious cells, and under the assumption of an excess of target cells. In a deterministic model, the ultimate in-host establishment or extinction of the viral infection depends entirely on whether the mean reproduction number is greater than, or less than, one, respectively. Here, the probability of extinction is determined by the probability distribution of R, not simply its mean value. In particular, we show that in some cases the probability of infection is not an increasing function of the mean reproduction number.

Multi-variate model of T cell clonotype competition and homeostasis.

Diversity of the naive T cell repertoire is maintained by competition for stimuli provided by self-peptides bound to major histocompatibility complexes (self-pMHCs). We extend an existing bi-variate competition model to a multi-variate model of the dynamics of multiple T cell clonotypes which share stimuli. In order to understand the late-time behaviour of the system, we analyse: (i) the dynamics until the extinction of the first clonotype, (ii) the time to the first extinction event, (iii) the probability of extinction of each clonotype, and (iv) the size of the surviving clonotypes when the first extinction event takes place. We also find the probability distribution of the number of cell divisions per clonotype before its extinction. The mean size of a new clonotype at quasi-steady state is an increasing function of the stimulus available to it, and a decreasing function of the fraction of stimuli it shares with other clonotypes. Thus, the probability of, and time to, extinction of a new clonotype entering the pool of T cell clonotypes is determined by the extent of competition for stimuli it experiences and by its initial number of cells.

Why are cell populations maintained via multiple compartments?

We consider the maintenance of 'product' cell populations from 'progenitor' cells via a sequence of one or more cell types, or compartments, where each cell's fate is chosen stochastically. If there is only one compartment then large amplification, that is, a large ratio of product cells to progenitors comes with disadvantages. The product cell population is dominated by large families (cells descended from the same progenitor) and many generations separate, on average, product cells from progenitors. These disadvantages are avoided using suitably constructed sequences of compartments: the amplification factor of a sequence is the product of the amplification factors of each compartment, while the average number of generations is a sum over contributions from each compartment. Passing through multiple compartments is, in fact, an efficient way to maintain a product cell population from a small flux of progenitors, avoiding excessive clonality and minimizing the number of rounds of division en route. We use division, exit and death rates, estimated from measurements of single-positive thymocytes, to choose illustrative parameter values in the single-compartment case. We also consider a five-compartment model of thymocyte differentiation, from double-negative precursors to single-positive product cells.

Counting generations in birth and death processes with competing Erlang and exponential waiting times.

Lymphocyte populations, stimulated in vitro or in vivo, grow as cells divide. Stochastic models are appropriate because some cells undergo multiple rounds of division, some die, and others of the same type in the same conditions do not divide at all. If individual cells behave independently, then each cell can be imagined as sampling from a probability density of times to division and death. The exponential density is the most mathematically and computationally convenient choice. It has the advantage of satisfying the memoryless property, consistent with a Markov process, but it overestimates the probability of short division times. With the aim of preserving the advantages of a Markovian framework while improving the representation of experimentally-observed division times, we consider a multi-stage model of cellular division and death. We use Erlang-distributed (or, more generally, phase-type distributed) times to division, and exponentially distributed times to death. We classify cells into generations, using the rule that the daughters of cells in generation n are in generation [Formula: see text]. In some circumstances, our representation is equivalent to established models of lymphocyte dynamics. We find the growth rate of the cell population by calculating the proportions of cells by stage and generation. The exponent describing the late-time cell population growth, and the criterion for extinction of the population, differs from what would be expected if N steps with rate [Formula: see text] were equivalent to a single step of rate [Formula: see text]. We link with a published experimental dataset, where cell counts were reported after T cells were transferred to lymphopenic mice, using Approximate Bayesian Computation. In the comparison, the death rate is assumed to be proportional to the generation and the Erlang time to division for generation 0 is allowed to differ from that of subsequent generations. The multi-stage representation is preferred to a simple exponential in posterior distributions, and the mean time to first division is estimated to be longer than the mean time to subsequent divisions.

Diffusion in a disk with inclusion: Evaluating Green's functions.

We give exact Green's functions in two space dimensions. We work in a scaled domain that is a circle of unit radius with a smaller circular "inclusion", of radius a, removed, without restriction on the size or position of the inclusion. We consider the two cases where one of the two boundaries is absorbing and the other is reflecting. Given a particle with diffusivity D, in a circle with radius R, the mean time to reach the absorbing boundary is a function of the initial condition, given by the integral of Green's function over the domain. We scale to a circle of unit radius, then transform to bipolar coordinates. We show the equivalence of two different series expansions, and obtain closed expressions that are not series expansions.

Step-by-step comparison of ordinary differential equation and agent-based approaches to pharmacokinetic-pharmacodynamic models.

Mathematical models in oncology aid in the design of drugs and understanding of their mechanisms of action by simulation of drug biodistribution, drug effects, and interaction between tumor and healthy cells. The traditional approach in pharmacometrics is to develop and validate ordinary differential equation models to quantify trends at the population level. In this approach, time-course of biological measurements is modeled continuously, assuming a homogenous population. Another approach, agent-based models, focuses on the behavior and fate of biological entities at the individual level, which subsequently could be summarized to reflect the population level. Heterogeneous cell populations and discrete events are simulated, and spatial distribution can be incorporated. In this tutorial, an agent-based model is presented and compared to an ordinary differential equation model for a tumor efficacy model inhibiting the pERK pathway. We highlight strengths, weaknesses, and opportunities of each approach.

Quantification of Type I Interferon Inhibition by Viral Proteins: Ebola Virus as a Case Study.

Type I interferons (IFNs) are cytokines with both antiviral properties and protective roles in innate immune responses to viral infection. They induce an antiviral cellular state and link innate and adaptive immune responses. Yet, viruses have evolved different strategies to inhibit such host responses. One of them is the existence of viral proteins which subvert type I IFN responses to allow quick and successful viral replication, thus, sustaining the infection within a host. We propose mathematical models to characterise the intra-cellular mechanisms involved in viral protein antagonism of type I IFN responses, and compare three different molecular inhibition strategies. We study the Ebola viral protein, VP35, with this mathematical approach. Approximate Bayesian computation sequential Monte Carlo, together with experimental data and the mathematical models proposed, are used to perform model calibration, as well as model selection of the different hypotheses considered. Finally, we assess if model parameters are identifiable and discuss how such identifiability can be improved with new experimental data.

Quantifying T Cell Cross-Reactivity: Influenza and Coronaviruses.

If viral strains are sufficiently similar in their immunodominant epitopes, then populations of cross-reactive T cells may be boosted by exposure to one strain and provide protection against infection by another at a later date. This type of pre-existing immunity may be important in the adaptive immune response to influenza and to coronaviruses. Patterns of recognition of epitopes by T cell clonotypes (a set of cells sharing the same T cell receptor) are represented as edges on a bipartite network. We describe different methods of constructing bipartite networks that exhibit cross-reactivity, and the dynamics of the T cell repertoire in conditions of homeostasis, infection and re-infection. Cross-reactivity may arise simply by chance, or because immunodominant epitopes of different strains are structurally similar. We introduce a circular space of epitopes, so that T cell cross-reactivity is a quantitative measure of the overlap between clonotypes that recognize similar (that is, close in epitope space) epitopes.

A Stochastic Intracellular Model of Anthrax Infection With Spore Germination Heterogeneity.

We present a stochastic mathematical model of the intracellular infection dynamics of Bacillus anthracis in macrophages. Following inhalation of B. anthracis spores, these are ingested by alveolar phagocytes. Ingested spores then begin to germinate and divide intracellularly. This can lead to the eventual death of the host cell and the extracellular release of bacterial progeny. Some macrophages successfully eliminate the intracellular bacteria and will recover. Here, a stochastic birth-and-death process with catastrophe is proposed, which includes the mechanism of spore germination and maturation of B. anthracis. The resulting model is used to explore the potential for heterogeneity in the spore germination rate, with the consideration of two extreme cases for the rate distribution: continuous Gaussian and discrete Bernoulli. We make use of approximate Bayesian computation to calibrate our model using experimental measurements from in vitro infection of murine peritoneal macrophages with spores of the Sterne 34F2 strain of B. anthracis. The calibrated stochastic model allows us to compute the probability of rupture, mean time to rupture, and rupture size distribution, of a macrophage that has been infected with one spore. We also obtain the mean spore and bacterial loads over time for a population of cells, each assumed to be initially infected with a single spore. Our results support the existence of significant heterogeneity in the germination rate, with a subset of spores expected to germinate much later than the majority. Furthermore, in agreement with experimental evidence, our results suggest that most of the spores taken up by macrophages are likely to be eliminated by the host cell, but a few germinated spores may survive phagocytosis and lead to the death of the infected cell. Finally, we discuss how this stochastic modelling approach, together with dose-response data, allows us to quantify and predict individual infection risk following exposure.

Competitive binding of STATs to receptor phospho-Tyr motifs accounts for altered cytokine responses.

Cytokines elicit pleiotropic and non-redundant activities despite strong overlap in their usage of receptors, JAKs and STATs molecules. We use IL-6 and IL-27 to ask how two cytokines activating the same signaling pathway have different biological roles. We found that IL-27 induces more sustained STAT1 phosphorylation than IL-6, with the two cytokines inducing comparable levels of STAT3 phosphorylation. Mathematical and statistical modeling of IL-6 and IL-27 signaling identified STAT3 binding to GP130, and STAT1 binding to IL-27Rα, as the main dynamical processes contributing to sustained pSTAT1 levels by IL-27. Mutation of Tyr613 on IL-27Rα decreased IL-27-induced STAT1 phosphorylation by 80% but had limited effect on STAT3 phosphorgylation. Strong receptor/STAT coupling by IL-27 initiated a unique gene expression program, which required sustained STAT1 phosphorylation and IRF1 expression and was enriched in classical Interferon Stimulated Genes. Interestingly, the STAT/receptor coupling exhibited by IL-6/IL-27 was altered in patients with systemic lupus erythematosus (SLE). IL-6/IL-27 induced a more potent STAT1 activation in SLE patients than in healthy controls, which correlated with higher STAT1 expression in these patients. Partial inhibition of JAK activation by sub-saturating doses of Tofacitinib specifically lowered the levels of STAT1 activation by IL-6. Our data show that receptor and STATs concentrations critically contribute to shape cytokine responses and generate functional pleiotropy in health and disease.

Fusion and fission events regulate endosome maturation and viral escape.

Endosomes are intracellular vesicles that mediate the communication of the cell with its extracellular environment. They are an essential part of the cell's machinery regulating intracellular trafficking via the endocytic pathway. Many viruses, which in order to replicate require a host cell, attach themselves to the cellular membrane; an event which usually initiates uptake of a viral particle through the endocytic pathway. In this way viruses hijack endosomes for their journey towards intracellular sites of replication and avoid degradation without host detection by escaping the endosomal compartment. Recent experimental techniques have defined the role of endosomal maturation in the ability of enveloped viruses to release their genetic material into the cytoplasm. Endosome maturation depends on a family of small hydrolase enzymes (or GTPases) called Rab proteins, arranged on the cytoplasmic surface of its membrane. Here, we model endosomes as intracellular compartments described by two variables (its levels of active Rab5 and Rab7 proteins) and which can undergo coagulation (or fusion) and fragmentation (or fission). The key element in our approach is the "per-cell endosomal distribution" and its dynamical (Boltzmann) equation. The Boltzmann equation allows us to derive the dynamics of the total number of endosomes in a cell, as well as the mean and the standard deviation of its active Rab5 and Rab7 levels. We compare our mathematical results with experiments of Dengue viral escape from endosomes. The relationship between endosomal active Rab levels and pH suggests a mechanism that can account for the observed variability in viral escape times, which in turn regulate the viability of a viral intracellular infection.

Quantification of Ebola virus replication kinetics in vitro.

Mathematical modelling has successfully been used to provide quantitative descriptions of many viral infections, but for the Ebola virus, which requires biosafety level 4 facilities for experimentation, modelling can play a crucial role. Ebola virus modelling efforts have primarily focused on in vivo virus kinetics, e.g., in animal models, to aid the development of antivirals and vaccines. But, thus far, these studies have not yielded a detailed specification of the infection cycle, which could provide a foundational description of the virus kinetics and thus a deeper understanding of their clinical manifestation. Here, we obtain a diverse experimental data set of the Ebola virus infection in vitro, and then make use of Bayesian inference methods to fully identify parameters in a mathematical model of the infection. Our results provide insights into the distribution of time an infected cell spends in the eclipse phase (the period between infection and the start of virus production), as well as the rate at which infectious virions lose infectivity. We suggest how these results can be used in future models to describe co-infection with defective interfering particles, which are an emerging alternative therapeutic.

Stochastic dynamics of Francisella tularensis infection and replication.

We study the pathogenesis of Francisella tularensis infection with an experimental mouse model, agent-based computation and mathematical analysis. Following inhalational exposure to Francisella tularensis SCHU S4, a small initial number of bacteria enter lung host cells and proliferate inside them, eventually destroying the host cell and releasing numerous copies that infect other cells. Our analysis of disease progression is based on a stochastic model of a population of infectious agents inside one host cell, extending the birth-and-death process by the occurrence of catastrophes: cell rupture events that affect all bacteria in a cell simultaneously. Closed expressions are obtained for the survival function of an infected cell, the number of bacteria released as a function of time after infection, and the total bacterial load. We compare our mathematical analysis with the results of agent-based computation and, making use of approximate Bayesian statistical inference, with experimental measurements carried out after murine aerosol infection with the virulent SCHU S4 strain of the bacterium Francisella tularensis, that infects alveolar macrophages. The posterior distribution of the rate of replication of intracellular bacteria is consistent with the estimate that the time between rounds of bacterial division is less than 6 hours in vivo.

Fate of a Naive T Cell: A Stochastic Journey.

The homeostasis of T cell populations depends on migration, division and death of individual cells (1). T cells migrate between spatial compartments (spleen, lymph nodes, lung, liver, etc.), where they may divide or differentiate, and eventually die (2). The kinetics of recirculation influences the speed at which local infections are detected and controlled (3). New experimental techniques have been developed to measure the lifespan of cells, and their migration dynamics; for example, fluorescence-activated cell sorting (4), in vitro time-lapse microscopy (5), or in vivo stable isotope labeling (e.g., deuterium) (6). When combined with mathematical and computational models, they allow estimation of rates of migration, division, differentiation and death (6, 7). In this work, we develop a stochastic model of a single cell migrating between spatial compartments, dividing and eventually dying. We calculate the number of division events during a T cell's journey, its lifespan, the probability of dying in each compartment and the number of progeny cells. A fast-migration approximation allows us to compute these quantities when migration rates are larger than division and death rates. Making use of published rates: (i) we analyse how perturbations in a given spatial compartment impact the dynamics of a T cell, (ii) we study the accuracy of the fast-migration approximation, and (iii) we quantify the role played by direct migration (not via the blood) between some compartments.

First passage events in biological systems with non-exponential inter-event times.

It is often possible to model the dynamics of biological systems as a series of discrete transitions between a finite set of observable states (or compartments). When the residence times in each state, or inter-event times more generally, are exponentially distributed, then one can write a set of ordinary differential equations, which accurately describe the evolution of mean quantities. Non-exponential inter-event times can also be experimentally observed, but are more difficult to analyse mathematically. In this paper, we focus on the computation of first passage events and their probabilities in biological systems with non-exponential inter-event times. We show, with three case studies from Molecular Immunology, Virology and Epidemiology, that significant errors are introduced when drawing conclusions based on the assumption that inter-event times are exponentially distributed. Our approach allows these errors to be avoided with the use of phase-type distributions that approximate arbitrarily distributed inter-event times.

Some deterministic and stochastic mathematical models of naïve T-cell homeostasis.

Humans live for decades, whereas mice live for months. Over these long timescales, naïve T cells die or divide infrequently enough that it makes sense to approximate death and division as instantaneous events. The population of T cells in the body is naturally divided into clonotypes; a clonotype is the set of cells that have identical T-cell receptors. While total numbers of cells, such as naïve CD4+ T cells, are large enough that ordinary differential equations are an appropriate starting point for mathematical models, the numbers of cells per clonotype are not. Here, we review a number of basic mathematical models of the maintenance of clonal diversity. As well as deterministic models, we discuss stochastic models that explicitly track the integer number of naïve T cells in many competing clonotypes over the lifetime of a mouse or human, including the effect of waning thymic production. Experimental evaluation of clonal diversity by bulk high-throughput sequencing has many difficulties, but the use of single-cell sequencing is restricted to numbers of cells many orders of magnitude smaller than the total number of T cells in the body. Mathematical questions associated with extrapolating from small samples are therefore key to advances in understanding the diversity of the repertoire of T cells. We conclude with some mathematical models on how to advance in this area.

A Novel Stochastic Multi-Scale Model of Francisella tularensis Infection to Predict Risk of Infection in a Laboratory.

We present a multi-scale model of the within-phagocyte, within-host and population-level infection dynamics of Francisella tularensis, which extends the mechanistic one proposed by Wood et al. (2014). Our multi-scale model incorporates key aspects of the interaction between host phagocytes and extracellular bacteria, accounts for inter-phagocyte variability in the number of bacteria released upon phagocyte rupture, and allows one to compute the probability of response, and mean time until response, of an infected individual as a function of the initial infection dose. A Bayesian approach is applied to parameterize both the within-phagocyte and within-host models using infection data. Finally, we show how dose response probabilities at the individual level can be used to estimate the airborne propagation of Francisella tularensis in indoor settings (such as a microbiology laboratory) at the population level, by means of a deterministic zonal ventilation model.