Table of contents

The evolution of multicellular organisms from unicellular counterparts involved a transition in Darwinian individuality from single cells to groups. A particular challenge is to understand the nature of the earliest groups, the causes of their evolution, and the opportunities for emergence of Darwinian properties. Here we outline a conceptual framework based on a logical set of possible pathways for evolution of the simplest self-replicating groups. Central to these pathways is the recognition of a finite number of routes by which genetic information can be transmitted between individual cells and groups. We describe the form and organization of each primordial group state and consider factors affecting persistence and evolution of the nascent multicellular forms. Implications arising from our conceptual framework become apparent when attempting to partition fitness effects at individual and group levels. These are discussed with reference to the evolutionary emergence of individuality and its manifestation in extant multicellular life—including those of marginal Darwinian status.

TRAIL (tumor necrosis factor-related apoptosis-inducing ligand) holds promise as an anti-cancer therapeutic but efficiently induces apoptosis in only a subset of tumor cell lines. Moreover, even in clonal populations of responsive lines, only a fraction of cells dies in response to TRAIL and individual cells exhibit cell-to-cell variability in the timing of cell death. Fractional killing in these cell populations appears to arise not from genetic differences among cells but rather from differences in gene expression states, fluctuations in protein levels and the extent to which TRAIL-induced death or survival pathways become activated. In this study, we ask how cell-to-cell variability manifests in cell types with different sensitivities to TRAIL, as well as how it changes when cells are exposed to combinations of drugs. We show that individual cells that survive treatment with TRAIL can regenerate the sensitivity and death-time distribution of the parental population, demonstrating that fractional killing is a stable property of cell populations. We also show that cell-to-cell variability in the timing and probability of apoptosis in response to treatment can be tuned using combinations of drugs that together increase apoptotic sensitivity compared to treatment with one drug alone. In the case of TRAIL, modulation of cell-to-cell variability by co-drugging appears to involve a reduction in the threshold for mitochondrial outer membrane permeabilization.

Directed cell migration often involves at least two types of cell motility that include multicellular streaming and chain migration. However, what is unclear is how cell contact dynamics and the distinct microenvironments through which cells travel influence the selection of one migratory mode or the other. The embryonic and highly invasive neural crest (NC) are an excellent model system to study this question since NC cells have been observed in vivo to display both of these types of cell motility. Here, we present data from tissue transplantation experiments in chick and in silico modeling that test our hypothesis that cell contact dynamics with each other and the microenvironment promote and sustain either multicellular stream or chain migration. We show that when premigratory cranial NC cells (at the pre-otic level) are transplanted into a more caudal region in the head (at the post-otic level), cells alter their characteristic stream behavior and migrate in chains. Similarly, post-otic NC cells migrate in streams after transplantation into the pre-otic hindbrain, suggesting that local microenvironmental signals dictate the mode of NC cell migration. Simulations of an agent-based model (ABM) that integrates the NC cell behavioral data predict that chain migration critically depends on the interplay of biased cell–cell contact and local microenvironment signals. Together, this integrated modeling and experimental approach suggests new experiments and offers a powerful tool to examine mechanisms that underlie complex cell migration patterns.

Mitochondria not only govern energy production, but are also involved in crucial cellular signalling processes. They are one of the most important organelles determining the Ca²⁺ regulatory pathway in the cell. Several mathematical models explaining these mechanisms were constructed, but only few of them describe interplay between calcium concentrations in endoplasmic reticulum (ER), cytoplasm and mitochondria. Experiments measuring calcium concentrations in mitochondria and ER suggested the existence of cytosolic microdomains with locally elevated calcium concentration in the nearest vicinity of the outer mitochondrial membrane. These intermediate physical connections between ER and mitochondria are called MAM (mitochondria-associated ER membrane) complexes. We propose a model with a direct calcium flow from ER to mitochondria, which may be justified by the existence of MAMs, and perform detailed numerical analysis of the effect of this flow on the type and shape of calcium oscillations. The model is partially based on the Marhl et al model. We have numerically found that the stable oscillations exist for a considerable set of parameter values. However, for some parameter sets the oscillations disappear and the trajectories of the model tend to a steady state with very high calcium level in mitochondria. This can be interpreted as an early step in an apoptotic pathway.

The immune system can detect most invading pathogens. The potential for detection of pathogens is dependent on the somatic diversity of the immune repertoires. While it is known that this somatic diversity is carefully generated, it is unclear how the diversity is distributed in the different genes encoding receptors of immune cells. Utilizing different metrics for richness and diversity at the level of small sequence fragments, we present here an analysis of the entire known human germline repertoire as represented by the sequences from the ImMunoGeneTics database of immune receptors. We have developed a fragment sequence quantification analysis to track variation of repertoires with different degrees of precision. Somatic diversity has previously been functionally characterized mostly by division of the V gene sequences into the more conserved and invariant framework (FR) of the receptor and more varied complementarity determining regions (CDR), that interact with the antigen. We find that CDR and FR can be explicitly identified with our sequence fragment diversity quantification technique. In terms of diversity, CDR and FR are especially distinct in B cell V genes. T cell V genes show less of the CDR/FR periodicity but are more diverse overall. Our analysis further shows that there are other areas of diversity outside the CDR and FR that are found widely dispersed in T cell receptor V genes and more tightly focused in FR1 and FR3 in the B cell receptor V genes. The diversity we observe is not dependent on allelic differences nor is this diversity segregated by individual V gene families. We would thus expect that each individual exhibit a diversity equivalent to that of the entire potential repertoire.

The Raf/MEK/ERK cascade is one of the most studied and important signal transduction pathways. However, existing models largely ignore the existence of isoforms of the constituent kinases and their interactions. Here, we propose a model of the ERK cascade that includes heretofore neglected differences between isoforms of MEK. In particular, MEK1 is subject to a negative feedback from activated ERK, which is further conferred to MEK2 via hetero-dimerization. Specifically, ERK phosphorylates MEK1 at the residue Thr292, hypothetically creating an additional phosphatase binding site, accelerating MEK1 and MEK2 dephosphorylation. We incorporated these recently discovered interactions into a mathematical model of the ERK cascade that reproduces the experimental results of Catalanotti et al (2009 Nature Struct. Mol. Biol.16 294–303) and Kamioka et al (2010 J. Biol. Chem.285 33540–8). Furthermore, the model allows for predictions regarding the differences in the catalytic activity and function of the MEK isoforms. We propose that the MEK1/MEK2 ratio regulates the duration of the response, which increases with the level of MEK2 and decreases with the level of MEK1. In turn, the amplitude of the response is controlled by the total amount of the two isoforms. We confirm the proposed model structure performing a random parameter sampling, which led us to the conclusion that the sampled parameters, selected to properly reproduce wild-type (WT) cell behavior, to allow for qualitative reproduction of differences in behavior WT cells and cell mutants studied experimentally.

Bistable regulatory elements enhance heterogeneity in cell populations and, in multicellular organisms, allow cells to specialize and specify their fate. Our study demonstrates that in a system of bistable genetic switch, the noise characteristics control in which of the two epigenetic attractors the cell population will settle. We focus on two types of noise: the gene switching noise and protein dimerization noise. We found that the change of magnitudes of these noise components for one of the two competing genes introduces a large asymmetry of the protein stationary probability distribution and changes the relative probability of individual gene activation. Interestingly, an increase of noise associated with a given gene can either promote or suppress the activation of the gene, depending on the type of noise. Namely, each gene is repressed by an increase of its gene switching noise and activated by an increase of its protein-product dimerization noise. The observed effect was found robust to the large, up to fivefold deviations of the model parameters. In summary, we demonstrated that noise itself may determine the relative strength of the epigenetic attractors, which may provide a unique mode of control of cell fate decisions.

Cells often have tens of thousands of receptors, even though only a few activated receptors can trigger full cellular responses. Reasons for the overabundance of receptors remain unclear. We suggest that, under certain conditions, the large number of receptors can result in a competition among receptors to be the first to activate the cell. The competition decreases the variability of the time to cellular activation, and hence results in a more synchronous activation of cells. We argue that, in simple models, this variability reduction does not necessarily interfere with the receptor specificity to ligands achieved by the kinetic proofreading mechanism. Thus cells can be activated accurately in time and specifically to certain signals. We predict the minimum number of receptors needed to reduce the coefficient of variation for the time to activation following binding of a specific ligand. Furthermore, we predict the maximum number of receptors so that the kinetic proofreading mechanism still can improve the specificity of the activation. These predictions fall in line with experimentally reported receptor numbers for multiple systems.

The time-ordered product framework of quantum field theory can also be used to understand salient phenomena in stochastic biochemical networks. It is used here to derive Gillespie's stochastic simulation algorithm (SSA) for chemical reaction networks; consequently, the SSA can be interpreted in terms of Feynman diagrams. It is also used here to derive other, more general simulation and parameter-learning algorithms including simulation algorithms for networks of stochastic reaction-like processes operating on parameterized objects, and also hybrid stochastic reaction/differential equation models in which systems of ordinary differential equations evolve the parameters of objects that can also undergo stochastic reactions. Thus, the time-ordered product expansion can be used systematically to derive simulation and parameter-fitting algorithms for stochastic systems.

The immune response to a pathogen has two basic features. The first is the expansion of a few pathogen-specific cells to form a population large enough to control the pathogen. The second is the process of differentiation of cells from an initial naive phenotype to an effector phenotype which controls the pathogen, and subsequently to a memory phenotype that is maintained and responsible for long-term protection. The expansion and the differentiation have been considered largely independently. Changes in cell populations are typically described using ecologically based ordinary differential equation models. In contrast, differentiation of single cells is studied within systems biology and is frequently modeled by considering changes in gene and protein expression in individual cells. Recent advances in experimental systems biology make available for the first time data to allow the coupling of population and high dimensional expression data of immune cells during infections. Here we describe and develop population-expression models which integrate these two processes into systems biology on the multicellular level. When translated into mathematical equations, these models result in non-conservative, non-local advection-diffusion equations. We describe situations where the population-expression approach can make correct inference from data while previous modeling approaches based on common simplifying assumptions would fail. We also explore how model reduction techniques can be used to build population-expression models, minimizing the complexity of the model while keeping the essential features of the system. While we consider problems in immunology in this paper, we expect population-expression models to be more broadly applicable.

In this study, we investigated the rheology of a doublet that is an aggregate of two red blood cells (RBCs). According to previous studies, most aggregates in blood flow consist of RBC doublet-pairs and thus the understanding of doublet dynamics has scientific importance in describing its hemodynamics. The RBC aggregation tendency can be significantly affected by the cell's deformability which can vary under both physiological and pathological conditions. Hence, we conducted a two-dimensional simulation of doublet dynamics under a simple shear flow condition with different deformability between RBCs. To study the dissociation process of the doublet, we employed the aggregation model described by the Morse-type potential function, which is based on the depletion theory. In addition, we developed a new method of updating fluid property to consider viscosity difference between RBC cytoplasm and plasma. Our results showed that deformability difference between the two RBCs could significantly reduce their aggregating tendency under a shear condition of 50 s⁻¹, resulting in disaggregation. Since even under physiological conditions, the cell deformability may be significantly different, consideration of the difference in deformability amongst RBCs in blood flow would be needed for the hemodynamic studies based on a numerical approach.

Particle tracking experiments with high speed digital microscopy yield the positions and trajectories of lipid droplets inside living cells. Angular correlation analysis shows that the lipid droplets have uncorrelated motion at short time scales (τ < 1 ms) followed by anti-persistent motion for lag times in the range of 1 ⩽ τ ⩽ 10 ms. The angular correlation at longer time scales, τ > 10 ms, becomes persistent, indicating directed movement. The motion at all time scales is associated with the lipid droplets being tethered to and driven along the microtubule network. The point at which the angular correlation changes from anti-persistent to persistent motion corresponds to the cross over between sub-diffusive and super diffusive motion, as observed by mean square displacement analysis. Correct analysis of the angular correlations of the detector noise is found to be crucial in modelling the observed phenomena.

Thrombosis significantly contributes to cancer morbidity and mortality. The mechanism behind thrombosis in cancer may be circulating tissue factor (TF), as levels of circulating TF are associated with thrombosis. However, circulating TF antigen level alone has failed to predict thrombosis in patients with cancer. We hypothesize that coagulation factor levels regulate the kinetics of circulating TF-induced thrombosis. Coagulation kinetics were measured as a function of individual coagulation factor levels and TF particle concentration. Clotting times increased when pooled plasma was mixed at or above a ratio of 4:6 with PBS. Clotting times increased when pooled plasma was mixed at or above a ratio of 8:2 with factor VII-depleted plasma, 7:3 with factor IX- or factor X-depleted plasmas, or 2:8 with factor II-, V- or VIII-depleted plasmas. Addition of coagulation factors VII, X, IX, V and II to depleted plasmas shortened clotting and enzyme initiation times, and increased enzyme generation rates in a concentration-dependent manner. Only additions of factors IX and X from low-normal to high-normal levels shortened clotting times and increased enzyme generation rates. Our results demonstrate that coagulation kinetics for TF particles are controlled by factor IX and X levels within the normal physiological range. We hypothesize that individual patient factor IX and X levels may be prognostic for susceptibility to circulating TF-induced thrombosis.

Regulating physical size is an essential problem that biological organisms must solve from the subcellular to the organismal scales, but it is not well understood what physical principles and mechanisms organisms use to sense and regulate their size. Any biophysical size-regulation scheme operates in a noisy environment and must be robust to other cellular dynamics and fluctuations. This work develops theory of filament length regulation inspired by recent experiments on kinesin-8 motor proteins, which move with directional bias on microtubule filaments and alter microtubule dynamics. Purified kinesin-8 motors can depolymerize chemically-stabilized microtubules. In the length-dependent depolymerization model, the rate of depolymerization tends to increase with filament length, because long filaments accumulate more motors at their tips and therefore shorten more quickly. When balanced with a constant filament growth rate, this mechanism can lead to a fixed polymer length. However, the mechanism by which kinesin-8 motors affect the length of dynamic microtubules in cells is less clear. We study the more biologically realistic problem of microtubule dynamic instability modulated by a motor-dependent increase in the filament catastrophe frequency. This leads to a significant decrease in the mean filament length and a narrowing of the filament length distribution. The results improve our understanding of the biophysics of length regulation in cells.

Epigenetic mechanisms of silencing via heritable chromatin modifications play a major role in gene regulation and cell fate specification. We consider a model of epigenetic chromatin silencing in budding yeast and study the bifurcation diagram and characterize the bistable and the monostable regimes. The main focus of this paper is to examine how the perturbations altering the activity of histone modifying enzymes affect the epigenetic states. We analyze the implications of having the total number of silencing proteins, given by the sum of proteins bound to the nucleosomes and the ones available in the ambient, to be constant. This constraint couples different regions of chromatin through the shared reservoir of ambient silencing proteins. We show that the response of the system to perturbations depends dramatically on the titration effect caused by the above constraint. In particular, for a certain range of overall abundance of silencing proteins, the hysteresis loop changes qualitatively with certain jump replaced by continuous merger of different states. In addition, we find a nonmonotonic dependence of gene expression on the rate of histone deacetylation activity of Sir2. We discuss how these qualitative predictions of our model could be compared with experimental studies of the yeast system under anti-silencing drugs.

Forces and stresses generated by the action of myosin minifilaments are analyzed in idealized computer-generated actin bundles, and compared to results for isotropic actin networks. The bundles are generated as random collections of actin filaments in two dimensions with constrained orientations, crosslinked and attached to two fixed walls. Myosin minifilaments are placed on actin filament pairs and allowed to move and deform the network so that it exerts forces on the walls. The vast majority of simulation runs end with contractile minifilament stress, because minifilaments rotate into energetically stable contractile configurations. This process is aided by the bending and stretching of actin filaments, which accomodate minifilament rotation. Stresses for bundles are greater than those for isotropic networks, and antiparallel filaments generate more tension than parallel filaments. The forces transmitted by the actin network to the walls of the simulation cell often exceed the tension in the minifilament itself.

One of the key questions in cell biology is how organelles are passed from parent to daughter cells. To help address this question, I used Brownian dynamics to simulate lipid droplets as model organelles in populations of dividing cells. Lipid droplets are dynamic bodies that can form both de novo and by fission, they can also be depleted. The quantitative interplay among these three events is unknown but would seem crucial for controlling droplet distribution in populations of dividing cells. Surprisingly, of the three main events studied: biogenesis, fission, and depletion, the third played the key role in maintaining droplet organelle number-–and to a lesser extent volume-–in populations of dividing cells where formation events would have seemed paramount. In the case of lipid droplets, this provides computational evidence that they must be sustained, most likely through contacts with the endoplasmic reticulum. The findings also agree with video microscopy experiments over much shorter timescales where droplet depletion in fission yeast cells was not observed. In general, this work shows that organelle maintenance is invaluable and lack thereof cannot necessarily be compensated for by organelle formation. This study provides a time-accurate, physical-based template for long-term cell division studies.

Acute respiratory distress syndrome (ARDS) is acute lung failure secondary to severe systemic inflammation, resulting in a derangement of alveolar mechanics (i.e. the dynamic change in alveolar size and shape during tidal ventilation), leading to alveolar instability that can cause further damage to the pulmonary parenchyma. Mechanical ventilation is a mainstay in the treatment of ARDS, but may induce mechano-physical stresses on unstable alveoli, which can paradoxically propagate the cellular and molecular processes exacerbating ARDS pathology. This phenomenon is called ventilator induced lung injury (VILI), and plays a significant role in morbidity and mortality associated with ARDS. In order to identify optimal ventilation strategies to limit VILI and treat ARDS, it is necessary to understand the complex interplay between biological and physical mechanisms of VILI, first at the alveolar level, and then in aggregate at the whole-lung level. Since there is no current consensus about the underlying dynamics of alveolar mechanics, as an initial step we investigate the ventilatory dynamics of an alveolar sac (AS) with the lung alveolar spatial model (LASM), a 3D spatial biomechanical representation of the AS and its interaction with airflow pressure and the surface tension effects of pulmonary surfactant. We use the LASM to identify the mechanical ramifications of alveolar dynamics associated with ARDS. Using graphical processing unit parallel algorithms, we perform Bayesian inference on the model parameters using experimental data from rat lung under control and Tween-induced ARDS conditions. Our results provide two plausible models that recapitulate two fundamental hypotheses about volume change at the alveolar level: (1) increase in alveolar size through isotropic volume change, or (2) minimal change in AS radius with primary expansion of the mouth of the AS, with the implication that the majority of change in lung volume during the respiratory cycle occurs in the alveolar ducts. These two model solutions correspond to significantly different mechanical properties of the tissue, and we discuss the implications of these different properties and the requirements for new experimental data to discriminate between the hypotheses.

The islets of Langerhans, responsible for controlling blood glucose levels, are dispersed within the pancreas. A universal power law governing the fractal spatial distribution of islets in two-dimensional pancreatic sections has been reported. However, the fractal geometry in the actual three-dimensional pancreas volume, and the developmental process that gives rise to such a self-similar structure, has not been investigated. Here, we examined the three-dimensional spatial distribution of islets in intact mouse pancreata using optical projection tomography and found a power law with a fractal dimension of 2.1. Furthermore, based on two-dimensional pancreatic sections of human autopsies, we found that the distribution of human islets also follows a universal power law with a fractal dimension of 1.5 in adult pancreata, which agrees with the value previously reported in smaller mammalian pancreas sections. Finally, we developed a self-avoiding growth model for the development of the islet distribution and found that the fractal nature of the spatial islet distribution may be associated with the self-avoidance in the branching process of vascularization in the pancreas.

Recently, a large number of studies have been carried out on the early signatures of sudden regime shifts in systems as diverse as ecosystems, financial markets, population biology and complex diseases. The signatures of regime shifts in gene expression dynamics are less systematically investigated. In this paper, we consider sudden regime shifts in the gene expression dynamics described by a fold-bifurcation model involving bistability and hysteresis. We consider two alternative models, models 1 and 2, of competence development in the bacterial population B. subtilis and determine some early signatures of the regime shifts between competence and noncompetence. We use both deterministic and stochastic formalisms for the purpose of our study. The early signatures studied include the critical slowing down as a transition point is approached, rising variance and the lag-1 autocorrelation function, skewness and a ratio of two mean first passage times. Some of the signatures could provide the experimental basis for distinguishing between bistability and excitability as the correct mechanism for the development of competence.