François Bertaux

Mᵐᵉ Alessandra

M Mark

M Paul

M Emmanuel

M Jérémie

M Dirk

M Gregory

Alessandra Carbone

Mark Chaplain

Paul Macklin

Nick Jagiella

Géraldine Cellière

Paul Van Liedekerke

Noémie Boisset

François I Fages

Jim Pioche

single cells to spatially organized cell populations

Cell-intrinsic, non-environmental sources of cell-to-cell variability, such as stochastic gene expression, are increasingly recognized to play an important role in the dynamics of tissues, tumors, microbial communities… However, they are usually ignored or oversimplified in theoretical models of cell populations.

In this thesis, we propose a cell-based (each cell is represented individually), multi-scale (cellular decisions are controlled by biochemical reaction pathways simulated in each cell) approach to model the dynamics of cell populations. The main novelty compared to traditional approaches is that the fluctuations of protein levels driven by stochastic gene expression are systematically accounted for (i.e., for every protein in the modeled pathways). This enables to investigate the joint effect of cell-intrinsic and environmental sources of cell-tocell variability on cell population dynamics. Central to our approach is a parsimonious and principled parameterization strategy for stochastic gene expression models.

The approach is applied on two case studies. First, it is used to investigate the resistance of HeLa cells to the anti-cancer agent TRAIL, which can induce apoptosis specifically in cancer cells. A single-cell model of TRAIL-induced apoptosis is constructed and compared to existing quantitative, single-cell experimental data. The model explains fractional killing and correctly predicts transient cell fate inheritance and reversible resistance, two observed properties that are out of reach of previous models of TRAIL-induced apoptosis, which do not capture the dynamics of cell-to-cell variability. In a second step, we integrate this model into multi-cellular simulations to study TRAIL resistance in virtual scenarios constructed to help bridging the gap between standard in-vitro assays and the response of in-vivo tumors.

More precisely, we consider the long-term response of multi-cellular spheroids to repeated TRAIL treatments. Analysis of model simulations points to a novel, mechanistic explanation for transient resistance acquisition, which involves the targeted degradation of activated proteins and a differential turnover between pro-and anti-apoptotic proteins.

Second, we apply our approach to a synthetic spatial pa erning system in yeast cells developed by collaborators. Focusing first on a sensing circuit responding to a messenger molecule, we construct a single-cell model that accurately capture the response kinetics of the circuit as observed in flow cytometry data. We then integrate this model into multicellular simulations and show that the response of spatially-organized micro-colonies sub-mi ed to gradients of messenger molecules is correctly predicted. Finally, we incorporate a model of a killing circuit and compare the predicted pa erns of dead or alive cells with experimental data, yielding insights into how the circuit parameters translate into multicellular organization phenotypes.

Our modeling approach has the potential to accelerate the obtention of more quantitative and predictive models of cell populations that encompass the molecular scale.

Résumé

Les sources intrinsèques d'héterogénéité cellulaire, comme l'expression stochastique des gènes, sont de plus en plus reconnues comme jouant un rôle important dans la dynamique des tissus, tumeurs, communautés microbiennes… Cependant, elles sont souvent ignorées ou représentées de manière simpliste dans les modèles théoriques de populations de cellules.

Dans ce e thèse, nous proposons une approche cellule-centrée (chaque cellule est représentée de manière individuelle), multi-échelle (les décisions cellulaires sont placées sous le contrôle de voies de signalisation biochimiques simulées dans chaque cellule) pour modéliser la dynamique de populations de cellules. La nouveauté principale de ce e approche réside dans la prise en compte systématique (pour toutes les protéines modélisées) des fluctuations du niveau des protéines résultant de l'expression stochastique des gènes. Cela permet d'étudier l'effet combiné des causes intrinsèques et environnementales d'héterogénéité cellulaire sur la dynamique de la population de cellules. Un élément central de notre approche est une stratégie parsimonieuse pour a ribuer les paramètres de modèles d'expression stochastique des gènes. Nous appliquons ce e approche à deux cas d'étude. Nous considérons en premier la resistance à l'agent anti-cancer TRAIL, qui peut induire l'apoptose sélectivement dans les cellules cancéreuses. Nous construisons d'abord un modèle 'cellule unique' de l'apoptose induite par TRAIL et le comparons à des données existantes quantitatives. Le modèle explique la mort fractionnelle (le fait que seul une fraction des cellules meurent à la suite d'un traitement) et prédit correctement l'héritabilité transiente du destin cellulaire ainsi que l'acquisition transiente de résistance, deux propriétés observées mais hors de portée des modèles pré-existants, qui ne capturent pas la dynamique de l'héterogénéité cellulaire. Dans une seconde étape, nous intégrons ce modèle dans des simulations multi-cellulaires pour étudier la résistance à TRAIL dans des scénarios virtuels intermédiaires entre les études classiques in-vitro et la réponse de tumeurs in-vivo. Plus précisément, nous considérons la réponse en temps long de sphéroides multi-cellulaires à des traitements répétés de TRAIL. L'analyse de nos simulations permet de proposer une explication originale et méchanistique de l'acquisition transiente de résistance, impliquant la dégradation ciblée des protéines activées et un différentiel dans le renouvellement des protéines pro-et anti-apoptotiques. Nous appliquons aussi notre approche à un système synthétique de création de motifs développé dans des levures par des collaborateurs. Nous nous concentrons d'abord sur un circuit senseur d'une molécule messager pour lequel nous construisons un modèle cellule unique qui capture de manière fine la dynamique de réponse du circuit telle qu'observée par cytométrie en flux. Nous intégrons ensuite ce modèle dans des des simulations multi-cellulaires et montrons que la réponse de micro-colonies organisées spatialement et soumises à des gradients de molécule messager est correctement prédite. Finalement, nous incorporons un modèle d'un circuit de mort et comparons les motifs prédits de cellules mortes/vivantes avec des données expérimentales, nous perme ons de mieux comprendre comment les paramètres du circuit se traduisent en phénotypes d'organisation multi-cellulaire.

Notre approche peut contribuer à l'obtention de modèles de populations de cellules de plus en plus quantitatifs, prédictifs et qui englobent l'échelle moléculaire.

Chapter

Introduction

.

Motivations

Understanding the behavior of multi-cellular assemblies such as tissues, tumors or microbial communities is a central goal of biology, and improvements of this understanding often translate into medical or biotechnological advances. Mathematical modeling is an invaluable tool to help realizing this goal (Thorne et al.

; ). By formalizing into equations our current understanding on the mechanisms at play, it allows to evaluate the quality of this understanding, to identify knowledge gaps or inconsistencies with observations. It can also point out to new, highly informative experiments to perform, and permits to test hypotheses on yet unclear mechanisms.

Cell-based multi-scale modeling, a natural formalism to describe the dy-

namics of tissues

Although many models of tissues that abstract cellular identity have brought important biological insights (consider for example the seminal work on the chemical basis of morphogenesis ) provides a natural and appealing way to formalize our understanding of the dynamics of tissues. Indeed, this paradigm adopts a cell-centric perspective, reflecting the way modern biologists think about life and how biological knowledge is organized. Turing himself considered such cell-based modeling in his paper on morphogenesis:

In this section a mathematical model of the growing embryo will be described. […] the cell theory is recognized but the cells are idealized into geometrical points. […] one proceeds as with a physical theory and defines an entity called 'the state of the system'.

One then describes how that state is to be determined from the state at a moment very shortly before. […] the description of the state consists of two parts, the mechanical and the chemical. The mechanical part of the state describes the positions, masses, velocities and elastic properties of the cells, and the forces between them. […] The chemical part of the state is given […] as the chemical composition of each separate cell; the diffusibility of each substance between each two adjacent cells must also be given. […] This account of the problem omits many features, e.g. electrical properties and the internal structure of the cell. But even so it is a problem of formidable mathematical complexity. One cannot at present hope to make any progress with the understanding of such systems except in very simplified cases.

The cell-centric perspective provides a natural decomposition of the tissue modeling problem into sub-problems. A first sub-problem concerns the modeling of the internal dynamics of cells and how it drives cellular decisions. Then, one should describe how the composition of the cell environment impacts the cell internal dynamics. By essence, molecular systems biology a empts to explicit those two models in terms of molecular interactions (Aebersold ).

To enable the description of tissue dynamics, it is then needed to also model ) how cells modify their environment and ) how contacts between neighbor cells impacts their internal dynamics. Finally, one should describe the own dynamics of the extra-cellular environment, in which, for example, processes such as molecular diffusion, transport, or mechanical relaxation could take place. Modeling the la er sub-problems (impact of cells on environment, cell-cell interactions and the own dynamics of the environment) and coupling them to the former (internal dynamics of cells and how they are influenced by the environment) is the central goal of multi-cellular systems biology.

When all the sub-models are coherently integrated, the resulting model is by definition a cell-based multi-scale model. Such model provides a natural and intuitive formalization of how a tissue is believed to change with time.

A bo om-up strategy for the construction of predictive models of tissues

To this decomposition of the modeling problem also corresponds a natural, bo om-up strategy to build, extend and refine models of tissue dynamics (Galle et al. ).

The first step is to build predictive models of how individual cells evolve as a function of their environment. For this purpose, models of cell internal dynamics and their response to cell environment can be constructed based on existing knowledge. Their predictive power can then be assessed by precisely probing cellular processes in well defined environmental conditions and at the single-cell level (Spiller et al. ). Iteration between model re-finements and experiments can be performed until a deemed sufficient predictive power is obtained.

In a second step, sub-models for how cells influence their environment and how cells in contact interact can be constructed and integrated with single-cell models. The resulting cell-based multi-scale models can also be iteratively refined by considering experimental systems of increasingly complex cell-environment and cell-cell couplings.

Recent progress renders this strategy technically within our reach

While this strategy for formalizing, constructing and improving our knowledge of tissue dynamics is natural and appealing, its practical applicability is hindered by important technical limitations. Among those limitations, the most critical are of experimental nature.

Indeed, while our ability to understand and predict the evolution of individual cells is a cornerstone of cell-based multi-scale modeling, it was for a long time not possible to experimentally probe the behavior of single cells with precision and scope (Elowi et al.

).

Thus, current knowledge on cell biology was in large part obtained from population measurements, in which a large number of seemingly similar cells are submi ed to a treatment of interest, and then are lyzed such that their molecular constituents can be pooled together for analysis. In those conditions, interpretation of the results is difficult: biologists naturally and intuitively reason with a cell-level perspective, but available measurements only reflect a population average of the behavior of many cells (Batchelor et al.

).

In the past decade, tremendous progress in experimental techniques was realized and enabled experimentalists to probe the behavior of single cells with molecular detail, and in cer- The utilization of those techniques enabled the investigation of cellular processes at the single-cell level and led to major insights into the underlying mechanisms. An illustrative example is the study of extrinsic apoptosis (i.e. programmed cell death induced by external factors) in mammalian cells. Using fluorescent reporters to measure the activity of caspases (proteins that control the initiation and execution of apoptosis via their protease activity), it was found that the initiation phase was slow and highly variable between cells while the second phase, leading to irreversible commitment to death, happened very suddenly and showed li le cell-to-cell variability (Albeck, Burke, Aldridge, et al. ; Hellwig et al. ).

Such behavior, which is thought to be fundamental to allow a clear, unambiguous cellular decision between life and death, is masked in population-level measurements of caspases activity kinetics.

Aside of experimental limitations, our ability to manipulate cell-based multi-scale models is restricted by the computational cost of numerical simulations. Indeed, because those models are in general very difficult to study analytically, computational costs directly limit the complexity and the size of the models we can investigate. However, while this prevent for example the simulation of many cells with full molecular details (the first whole-cell model (Karr et al. ) only simulated an isolated M. genitalium cell, although hundreds of times), it is already possible to simulate very large models and the available computational power is still growing at a high rate.

In summary, the main technical limitation to the implementation of the bo om-up, cellcentric strategy to understand the dynamics of tissues, that is the difficulty to probe the behavior of single cells in controlled environments, was markedly overcome in the past decade.

However, as will be now discussed, the concommitently emerging picture of how single cells behave also seriously challenges our ability to model this behavior. In other words, out of sight behind important technical barriers, were standing at least equally important methodological barriers.

The methodological challenge posed by cell-to-cell variability

The relevance of the cell-based multi-scale modeling paradigm depends on our ability to predict the evolution of single cells in well controlled environmental conditions. In cases of failure, the modeling of tissue dynamics could still be a empted and achieved with approaches that abstract cellular identity. However, it would likely prevent the much needed gap bridging between the molecular and tissue scales (Wolkenhauer et al.

).

Paradoxically, as we are more and more able to observe and measure single cell behaviors, they seem less and less predictable. Indeed, there is a growing awareness that cell-to-cell variability is widespread and constitutes a fundamental property of cells. It refers to the fact that a priori similar cells cultivated in the same environmental conditions could develop important phenotypic heterogeneity and take very different cellular decisions in response to signals (Balázsi et al. ).

This phenomenon challenges our ability to build predictive models of single cell behavior.

To facilitate the description of those challenges, we will consider an illustrative example. A population of cancer cells is submi ed to a cytotoxic drug. We assume that all cells always sense the same environment, i.e. that cell-extrinsic sources of cell-to-cell variability can be excluded. Although at the time of treatment there are no observable differences between cells, % die rapidly after treatment while others survive.

A first a itude a modeler could adopt in this situation is to consider that because 'identical' cells (from the viewpoint of the observer) respond differently to the treatment, the model should include some intrinsic noise (or stochasticity) to account for cell fate variability. Concretely, if the model for cell dynamics and its response to drug is phenomenologic, one will simply a ribute a probability of responding to the drug and apply it to each cell. If on the contrary the model is mechanistic (i.e. a empts to explicit chemical reactions at play), one can incorporate stochasticity into an originally deterministic model so as to account for molecular noise (Wilkinson ).

This approach was in many cases successful to explain cell-to-cell variability observa- ) after exposure to very low levels of signal ligands.

Intrinsic noise is not the only possible explanation to cell-to-cell variability. Because live monitoring of the complete molecular state of single cells is way beyond our scope, it is possible that the observed apparition of cell-to-cell variability (in our example, variability in the decision between life and death) simply reflects pre-existing, hidden differences between cells (Ricicova et al. ).

What could indicate the presence of hidden heterogeneity? First, if a mechanistic model of how the drug acts to induce cell death is available, it could fail to reproduce sufficient cell fate variability under realistic assumptions for its parameters. In this case, one will search for plausible cellular quantities that could vary between cells and are implicated in the drug response. Those quantities might be already in the model but assumed not to vary between cells or they could lie at the boundaries of model scope, hence requiring a model extension.

Indication of hidden heterogeneity could also exist in the experimental data and can be detected without mechanistic knowledge. For example, if sibling cells are more likely to respond the same way to the drug than non-related cells, it shows that some level of hidden heterogeneity contributes to cell fate variability and that the related cellular quantities are, to some extent, conserved at division.

To model the impact of pre-existing heterogeneity (often referred to as extrinsic noise in such context) on the observed variability, the standard approach (Gaudet et al.

; Kallenberger et al. ; Shokhirev et al. ) is to sample the value of the initially variable cellular quantities from adequate distributions and to use it as initial condition for a model of cell response (this model can be deterministic or include some stochasticity representing intrinsic noise).

Situations in which pre-existing heterogeneity was demonstrated to be a major cause of cell-to-cell variability in response to a stimulus are numerous. For example, cancer cells treated with the apoptosis inducing ligand TRAIL commit to death after a highly variable delay that could be linked to pre-existing differences in the level of extrinsic apoptosis signaling proteins (Spencer et al.

). The predominant role of pre-existing heterogeneity was demonstrated by the observation that sister cells commit to death almost at the same time. In addition, when a deterministic model of TRAIL signaling was combined with randomly distributed initial protein levels, predicted variability in death commitment delays were similar to observations. Such non-genetic variability in response to TRAIL could have implications for cancer therapy, because many TRAIL-like drugs (i.e., that target the same receptors) were or are still being developed and tested in clinical trials (Dimberg et al.

).

Another therapeutically important example is the identification of bacterial persistence (or growth bistability, i.e. the pre-existence of a small fraction of non-growing cells) as a mechanism of resistance to antibiotics (Balaban et al.

).

In summary, modeling cell-to-cell variability in response to a stimulus is challenging because both intrinsic noise during the response and pre-existing, potentially un-observed heterogeneity could contribute. To address this, models of cell response (that can include stochasticity, see for example (Tay et al.

)) were combined with distributions of cell state to represent the pre-existing cell-to-cell differences. In many cases, this modeling approach enabled to successfully explain observations of cell-to-cell variability. However, as we will see, it suffers from a major shortcoming that prevent to use such models as single-cell models to incorporate into cell-based, multi-scale models of tissue dynamics.

Towards a dynamic and integrated view on cell-to-cell variability

The main limitation of the previous approach is that pre-existing heterogeneity is seen as a static parameter of the cell population. The origin of this heterogeneity is ignored, as well as the possibility that this heterogeneity might emerge from dynamic fluctuations in single cells.

However, with the notable exception of genetic information, many cellular quantities exhibit fluctuation timescales that can be of the order of the cell cycle or lower. For example, in mammalian cells, the mixing time (i.e. the typical time for single cells to reach high or low level of the population distribution) of protein levels was observed for a dozen of genes to range between one and three days (Sigal et al.

). In general, sources of diversification are inevitably at play within cells, because many constituents are in low copy number and thus subjected to molecular noise. While regulation mechanisms can be implemented to buffer the effect of molecular noise, completely supressing noise would have a tremendous energetic cost and would not be evolutionnary favored (Fraser et al. ).

When a cell population is submi ed to a stimulation, such as a TRAIL treatment, cells will respond differently depending on their initial state. If the response has an impact on the proliferative behavior of the cell or on its potential to undergo cell death, selection will operate and the distribution of cell states in the population will be affected. In the case of TRAIL-induced apoptosis, if the amount of a given protein has a strong and positive impact on the execution of apoptosis, its level in the population is expected to decrease after treatment (independently of the modifications that could affect this protein during signaling).

Because it does not account for diversification of cell state in absence of signaling, the approach that models pre-existing heterogeneity via distributions will grossly over-estimate selection effects: if a second treatment is applied on the same population, it will be predicted to have the same, diminished effect independently on the time interval between treatments, while in reality one would expect the recovery of the initial population distribution (and hence of the treatment impact) as time goes on. For example, transient resistance acquisition after treatment with anti-cancer drugs was observed in several occasions (for TRAIL (Flusberg et al.

), for inhibitors of EGFR signaling (Sun et al.

) or chemotherapeutic agents (Sharma et al.

)), and in each case pre-existing differences were shown to contribute to cell fate variability. Models that use the distribution approach to account for those differences would not be able to reproduce transient resistance acquisition.

This approach is therefore not suited for long-term simulation of cell populations submi ed to arbitrary stimulations, while this capability is required by the cell-based multi-scale paradigm for describing tissue dynamics. Instead, single-cell models should account for the diversification of the factors that in the presence of drug promote cell fate variability. In short, variability that was modeled as being extrinsic should be recasted as intrinsic noise, allowing to fully represent the dynamic nature of cell-to-cell variability in cell populations.

For example, if there is a set of proteins for which their levels are known to impact cell fate decision in response to a cytotoxic drug, and if the natural cell-to-cell variability in those levels is such that cell fate variability is expected; the single-cell fluctuations of those levels from which population variability originates should be explicitely modeled.

This would enable to study the interplay between diversification, selection and druginduced cellular responses. Such model could be used to predict the long-term outcome of repeated drug treatments on cancer cell populations cultivated in homogeneous and controlled conditions. Those predictions can permit the identification of efficient strategies for drug choice, combination, dosing and scheduling that can be tested and validated experimentally.

When such models are successfully obtained, they can be integrated in cell-based multiscale models in order to describe increasingly complex population contexts (i.e. with increased cell-environment and cell-cell couplings). For example, the response to repeated treatments of in-vivo tumors is likely to differ from the response of the same population cultivated in-vitro in homegeneous conditions, because of heterogeneous drug penetration or nutrients supply, secretion of factors by the surrounding tissues, and potentially many other factors.

By comparing model predictions with observations from adequate experimental systems (such as in-vitro grown multi-cellular spheroids, which provides an interesting intermediate between monolayer cultures and in-vivo systems (Hirschhaeuser et al. )), an increasingly comprehensive, multi-scale picture of the cell population dynamics could be achieved.

Integrating top-down and bo om-up description of tissues: expected insights from synthetic biology

From now on, we considered the problem of how to improve and formalize our understanding of the dynamics of existing, natural tissues and discussed the applicability of a natural, cell-centric, bo om-up strategy to do so. However, synthetic biology, a recent but fast growing field of biological research, aims to craft new, artificial living systems with desired properties and behavior (Benner & Sismour ).

For this purpose, synthetic biologists rely on the genetic modification of existing cells.

Therefore, the first successes of synthetic biology concerned the implementation of cell-level desired behavior, such as the production of high value molecules or the sensory detection of certain compounds. But nowadays synthetic biologists aim to extend their scope by implementing tissue-level behaviors (Purnick & Weiss ).

This goal is highly challenging because one needs to find the molecular-level, genetic modifications to perform in single cells that would enable a tissue-level desired behavior.

In other words, multi-cellular synthetic biology has to solve the "inverse" problem that multicellular systems biology aim to address. Thus, two distinct tasks should be performed: first, one needs to identify candidate genetic implementations that are potentially able to generate the desired tissue-level behavior; second, those implementations should be constructed and the emerging behavior experimentally tested. Because the second task is tedious, costly and time-consuming, success greatly depends on the quality with which the first task (the design problem) is achieved. The set of genetic implementations to test should be minimized while their likelihood to be valid should be maximized.

Without a rational, model-based approach, meeting those constraints might be impos-sible in practice as soon as the desired behavior has a certain degree of complexity. Cellbased multi-scale modeling is well suited for such model-based design because it allows to decompose the design problem into two sub-problems. First, using abstract, phenomenological models of cell behavior, cell-cell and cell-environment couplings, one can identify high-level cell behaviors that would lead to the desired tissue-level behavior. Then, one would search for genetic implementations that can program those high-level cell behaviors.

This second step corresponds to the simpler (but still very challenging), historical design problem addressed by synthetic biologists.

Because this top-down design strategy rely on abstract, high-level phenomenological models, deviations between the real and desired cell and tissue behavior are inevitably expected. Thus, it should be complemented by a bo om-up modeling approach to characterize the real behavior of a given genetic implementation of the system. This would allow to debug and tune the system until the desired behavior is achieved. In addition, by comparing and linking the top-down and bo om-up model descriptions of the system, one could gain insights on the functioning of natural tissues. For example, links between highlevel tissue behavior and molecular networks into cells could shed light on the evolution of multi-cellular organisms from uni-cellular organisms.

.

Contributions

The principal contribution of this thesis is a cell-based modeling approach to simulate the dynamics of multi-cellular populations with molecular details. The approach can be conceptually divided in two main steps. First, in each cell, cellular decisions are driven by quantitative models of biochemical reaction pathways. Importantly, stochastic fluctuations in protein levels arising from transcriptional noise are modeled for all the proteins involved in those pathways. Second, those single-cell models are integrated into multi-cellular simulations in order to investigate the emergent population behavior when cell-cell and cellenvironment couplings are added.

Cell-based, multi-scale modeling approach encompassing molecular details have already ). However, a main novelty of our approach is the fact that protein fluctuations are accounted for all modeled proteins. This is important because those fluctuations intrinsically generate cell-to-cell variability, even when cells sense homogenous and constant environments. If the environment is changing, this variability will be dynamically reshaped by the interplay between environment sensing and protein fluctuations. Moreover, if the environment is also heterogenous among cells, both cell-intrinsic and environmental sources of cell-to-cell variability will impact the behavior of the whole population. Systematically modeling those fluctuations is challenging, because the parameters governing the stochas-tic expression of proteins are gene-specific, largely unknown and difficult to infer (Sherman

& Cohen

). Building upon analytical results on a popular stochastic gene expression model, we have developped a parsimonious and principled method to enable this systematic parameterization.

We applied our cell-based multi-scale modeling approach to two biological problems.

The first problem concerns the resistance of cancer cell populations to the therapeutical agent TRAIL, which induces apoptosis selectively in cancer cells. The second problem is the modeling of a synthetic yeast pa erning system currently developped in the Weiss Laboratory at MIT.

In both cases, accounting for protein fluctuations appeared highly needed. For TRAILinduced apoptosis, evidence for an important role of cell-to-cell variability in protein levels was accumulating (Spencer et al.

; Rehm et al.

); yet, despite the existence of a detailed kinetic model of TRAIL apoptosis signaling (Albeck, Burke, Spencer, et al. ), no model was able to explain the observed fractional killing, or observations on the transient correlation between sister cell fate (Spencer et al.

). For the synthetic yeast pa erning system, high-throughput flow cytometry data revealed highly heterogenous response kinetics for a sensing circuit central to system functioning. Such heterogeneity is likely to impact the overall system behavior.

In both cases, integrating single-cell models into multi-cellular simulations also appeared critical. For TRAIL-induced apoptosis, it is required to investigate the long-term response of cancer cell populations repeatedly treated by TRAIL, and to study the contribution of cell-intrinsic and environmental effects in TRAIL resistance in realistic scenarios (i.e., that ressemble the treatment of in-vivo tumors). For the synthetic yeast pa erning system, the need to simulate spatially organized cell population is obvious as the long-term objective of the project is to obtain controlled, dynamic spatial pa erns within and between microcolonies.

In both cases, our method to account for protein fluctuations was successful. For TRAILinduced apoptosis, comparison with available single-cell data (Spencer et al.

) shows that it permits to greatly increase the prediction scope of a previous model of TRAIL apoptosis signaling (Albeck, Burke, Spencer, et al. ). Moreover, the obtained model also led to the finding that fractional killing and transient resistance acquisition after TRAIL treatment does not necessarily involve the parallel activation of others, pro-survival pathways as was predominantly proposed (Falschlehner et al.

) but could simply emerge from an interplay between protein fluctuations, selection and targeted degradation of pro-apoptotic forms. Together, those contributions resulted in a publication into a reference journal in the computational biology field (Bertaux et al.

). For the synthetic yeast pa erning system, it enabled the construction of a single-cell model of the sensing circuit that accurately capture the heterogenous response kinetics observed by flow cytometry.

In both cases, integrating the corresponding single-cell models into multi-cellular simulations also yielded important insights. For TRAIL-induced apoptosis, it predicted that resistance acquisition by non-mutational mechanisms can be very important, in agreement with observations (Flusberg et al. ), and shed light on the molecular mechanisms involved, such as the targeted degradation of activated pro-apoptotic form and an imbalance between the turnover of pro-and anti-apoptotic proteins. In addition, simulations of the treatment of multi-cellular spheroids suggested that limited TRAIL penetration is unlikely to be a driving cause of resistance, but that it can exacerbate the impact of cell-intrinsic resistance acquisition. For the synthetic yeast pa erning system, the multi-cellular simulations correctly predicted the spatio-temporal response of micro-colonies submi ed to gradients of messenger molecule, demonstrating that biological models can be quantitative from single cells to spatially organized populations.

Another contribution of this thesis is the developpment of software tools (for example, Chapter presents the results obtained by applying our cell-based multi-scale modeling approach to guide the development of an artificial yeast pa erning system. We start by describing the context of the collaboration, the genetic implementation of the system and the experimental data that will be confronted to our model. Focusing on a sensing circuit, which is a central component of the full system, we construct a single-cell model that quantitatively capture the behavior of the system as quantified by flow cytometry data. Then, we integrate into multi-cellular simulations and ask wether we can correctly predict the behavior of spatially organized micro-colonies responding to a gradient of messenger molecules.

a
Finally, we extend the model to account for the combined behavior of the sensing circuit with a cell death circuit and compare model predictions with experimental data.

In Chapter , we discuss methodological problems related to the description and simulation of cell-based multi-scale models such as the ones used in this thesis. We present algorithms for the simulation combined signal transduction and stochastic gene expression and discuss their computational efficiency. We then discuss a few important points regarding the simulation of cell populations. Finally, we consider the problem of minimizing the development time and code errors asociated to the simulation of cell-based multi-scale models and present a tool allowing the automated generation of simulation code from a programmatic model description.

Chapter concludes this thesis. We summarize the results and try to relate them with the current efforts of the quantitative biology and system biology communities. We then propose future directions of research. ). This selectivity, together with the independence of extrinsic apoptosis with respect to the mutational status of p (p mutations are a common cause of resistance to chemotherapy), suggests that using TRAIL-induced apoptosis could be an efficient strategy in cancer therapy (Johnstone et al. ).

Chapter

This led to dozens clinical trials involving drugs targeting TRAIL death receptors (such as TRAIL itself or antibodies) either in monotherapy or in combination with other anti-cancer agents. Most trials were disappointing: although the drugs were well tolerated, antitumor efficacy was usually limited (Johnstone et al. ; Dimberg et al. ). Because understanding the mechanisms behind TRAIL resistance could enable to design be er therapeutical strategies, a tremendous research effort was and is still undertaken to characterize TRAIL-induced apoptosis.

How does TRAIL work: molecular events from receptor binding to death commitment Most of the molecular events leading from TRAIL exposure to cell death are known (Johnstone et al. ). Figure . provides a stylized, simplified representation of those events. First, the binding of TRAIL to the death receptors DR or DR promotes their crosslinking. On the intra-cellular domain, this leads to the formation of the death-inducing signaling complex (DISC), which recruit and activate initiator caspases like caspase-(C ) (Falschlehner et al.

).

Once activated, initiator caspases cleave and activate effector caspases such as caspase-(C ). Effector caspases cleave essential structural proteins, inhibitors of DNase, and DNA repair proteins (PARP). Thus, as soon as the activity of effector caspases becomes important, the cell is irreversibly commi ed to apoptosis and the associated morphological changes usually occur within to minutes (Albeck, Burke, Aldridge, et al. ). The cellular effect of effector caspase activation is regulated by factors such as XIAP, which blocks the proteolytic activity of caspase-by binding tightly to its active site and promotes its degradation via ubiquitination (Deveraux et al.

). In addition to the direct activation of effector caspases, initiator caspases also activate Bid that then activates Bax (Luo et al. ).

If not kept in check by inhibitors, most notably Bcl , activated Bax directly contributes to the formation of pores in the mitochondria outer membrane, leading to mitochondrial outer membrane permeabilization (MOMP) (Kim et al. ). In most cells, MOMP is required to reach lethal activity of effector caspases. Following MOMP, critical apoptosis regulators, such as Smac and cytochrome c (CyC), translocate into the cytoplasm. Smac binds to and inactivates XIAP, thus relieving the inhibition of effector caspases by XIAP (Du et al.

).

Cytochrome c combines with Apaf-to form the apoptosome that in turn activates the initiator caspase-(C ) that activates effector caspases.

In summary, death signals are transmitted to effector caspases by initiator caspases either directly or indirectly via a mitochondrial pathway. The displayed interactions are abtracted from an established kinetic model (Albeck, Burke, Spencer, et al. 2008).

Non-genetic cell fate variability in TRAIL-induced apoptosis

Resistance of cancer cell populations to TRAIL-induced apoptosis is widely observed in both in-vitro and in-vivo se ings (Lemke, Karstedt, Zinngrebe, et al. ). Obviously, partial killing after TRAIL exposure has important therapeutical implications. Characterizing the mechanisms by which cells die in response to TRAIL is not sufficient: it is even more critical to understand why some cells can escape TRAIL-induced apoptosis.

The fate of a cell in response to TRAIL is not fully determined by its genotype. Indeed, even monoclonal populations cultivated in-vitro as monolayers usually exhibit partial responses to TRAIL treatment: only a fraction of the cells undergo apoptosis, even at saturating TRAIL doses (Spencer et al.

). Selection of subclones harboring new mutations is unlikely because the amount of such fractional killing seems to be a stable property of a given cell line: when two treatments are separated by at least a week, the same fraction of cells is killed in both cases (Flusberg & Sorger ).

However, if TRAIL treatment is repeated within a few days after initial exposure, less ).

A second type of defense against death signaling is the active degradation of activated proapoptotic proteins. Such degradation is achieved by the ubiquitin-proteasome pathway and mediated by proteins that bind the pro-apoptotic target and prime them for degradation (Vucic et al.

). For example, XIAP promotes the degradation of caspase-(Suzuki et al.

) and TRAF promotes the degradation of caspase- (Gonzalvez et al. ). Note that because this mechanism requires binding to the pro-apoptotic protein, it could have a dual action against death signaling: competitive binding and degradation (as for XIAP). Because of the variety and complexity of the molecular mechanisms that seem involved in determining a cell response to TRAIL exposure, there is no hope to gain a comprehensive, mechanistic picture of the cell decision between life and death without adopting a quantitative, system-level approach. In addition, because non-genetic cell-to-cell variability seems to impact this decision, both experimental and theoretical investigations probably cannot rely solely on population level measurements and should acknowledge the single-cell level.

Co-treatment with a protein synthesis inhibitor: a valuable tool to investigate core TRAIL apoptosis signaling

A common practice in the experimental investigation of extrinsic apoptosis is to block protein synthesis when cells are stimulated with death signals. For example, the first singlecell studies of TRAIL-induced caspases activation, in which fluorescent reporters (FRET) allowed to monitor caspase activity dynamics, used cycloheximide (CHX, a very potent inhibitor of protein synthesis) in co-treatment with TRAIL (Albeck, Burke, Aldridge, et al.

; Hellwig et al. ).

The principal motivation for this is to remove the influence of TRAIL-induced protein synthesis by pro-survival signaling. Therefore, changes in protein levels and activity are necessary mediated by protein-protein interactions and protein degradation. This facilitates the interpretation of experimental results because the core apoptosis signaling is functioning in isolation with respect to transcription-or translation-mediated regulations. In addition, because gene expression is generally noisy (gene and mRNAs are in low copy number in cells), blocking protein synthesis removes a source of noise in cell response to TRAIL.

Those single-cell experiments led to important insights into the timing of initiator and activator caspases activitation in response to TRAIL and CHX (Albeck, Burke, Aldridge, et al. ; Hellwig et al. ) in HeLa cells. First, the pre-MOMP delay is highly variable between cells. During that period, initiator caspases (caspase-) activity steadily rises while effector caspases (caspase-) are not active. After MOMP, a brutal, snap-action activation of effector caspases occur in all cells and with highly similar kinetics. Further experiments confirmed the role of XIAP and targeted degradation of caspase-in preventing caspaseactivity during the pre-MOMP period.

Those precise single-cell observations were confronted to an ODE, mass-action kinetic model of TRAIL-induced apoptosis (Albeck, Burke, Spencer, et al. ). The model is exhaustive in the sense that it represents all events from TRAIL binding to receptors until effector caspases activation, including MOMP-related events. It describes reactions between protein species, including native proteins. Despite this large size, it is not completely mechanistic, as some events such as receptor oligomerization upon TRAIL binding and DISC assembly are simplified into a few reactions. The model has important capabilities: it can reproduce the snap-action control over effector caspases activation and its regulation by Bcl and smac. Because it was constructed to reproduce experiments involving TRAIL and CHX co-treatments, it does not include protein synthesis.

In summary, single-cell experiments on TRAIL-induced apoptosis conducted in absence of protein synthesis, when complemented with kinetic modeling, enabled to gain a compre-hensive and relatively detailed picture of the molecular events that connect receptor binding to cell death commitment. However, because in those conditions all cells eventually die, direct insights on the mechanisms behind fractional killing were out of reach.

The role of protein level variability in death time variability

Fractional killing is the manifestation of cell-to-cell variability in the decision between life and death. As we already mentionned, cell-to-cell variability is also observed among cells that die: they commit to death after a highly variable delay from one another, including conditions in which protein synthesis is blocked (Albeck, Burke, Spencer, et al. ; Spencer et al. ; Rehm et al. ). Therefore, such variability cannot be explained by differences in translation-dependent TRAIL-induced survival responses.

Rather, it was proposed to originate mostly from pre-existing differences in the levels of proteins composing the apoptosis pathway. This explanation was supported by modeling: when taken as initial conditions of the previously described model from Albeck and colleagues (which is deterministic), differences in protein levels are sufficient to explain observed variability in death times (Spencer et al.

Adopting a dynamic view on protein level variability by modeling fluctuations rather than distributions

To quantitatively assess the loss of death time synchrony between sister cells as their age at treatment time increases, the previously mentioned modeling approach, where only the consequences of protein fluctuations (cell-to-cell variability at a fixed time) are accounted for, is inadequate. Instead, protein fluctuations themselves should be modeled.

The need to account for protein fluctuations is even more stringent when considering observations after treatment with TRAIL alone. In that case, protein synthesis is not blocked and thus can impact the decision between life and death. Even if TRAIL does not change protein synthesis (via induction of survival pathways), significant differences with TRAIL and CHX treatments are expected, as the constitutively noisy protein synthesis will interact with signaling reactions.

Thus, before complexifying the model to account for eventual regulations via survival pathways induction, it is critical to assess how much can be explained when protein synthesis is not altered by TRAIL signaling. In this chapter, we investigate this question by enabling a kinetic model of TRAIL-induced apoptosis with stochastic protein turnover for all proteins, following a generic and principled approach. It enriches the model with a fundamental property as the dynamics of cell-to-cell variability is represented, allowing disentangling the effects of constitutive protein fluctuations, signaling protein-protein reactions and potentially induced changes in protein synthesis.

. Model construction using a generic and principled ap- ).

One approach to account for those differences is to incorporate protein level variability as random initial conditions of an ODE model describing the signaling reactions ("extrinsic noise approach") (Gaudet et al. ). However, variability is imposed at time zero and then behavior is deterministic: thus, as we argued in this thesis introduction, it is not appropriate to study transduction on long time scales, during which protein levels dynamically fluctuate (Sigal et al. ).

A more natural manner to account for protein level variability is to represent their stochastic synthesis and degradation ("intrinsic noise approach"). Although several studies did account for cell-to-cell differences in protein levels in an extrinsic, static manner via random initial conditions (Spencer et al.

; Gaudet et

Modeling the fluctuations of individual proteins: a popular model of stochastic gene expression

We model protein fluctuations arising from gene expression noise with a stochastic process describing mRNA level fluctuations (promoter activity switches, mRNA production and degradation are stochastic events), and deterministic processes for protein translation and degradation (Figure . ). The stochastic reactions are triggered at a constant rate (or propensity), therefore, the time between two stochastic events is exponentially distributed. For example, the time needed for an inactive gene to become active is exponentially distributed with average T off = 1 k on and the time needed to produce one mRNA when the gene is active is also exponentially distributed with average equal to

Integration of individual fluctuation models within kinetic models of signal transduction pathways

Starting from an ODE, kinetic model describing the protein-protein reactions involved in signal transduction, each species corresponding to a protein in its native form is equipped with its own instance of the stochastic protein turnover model.

Those stochastic models are then integrated within the ODE model of signal transduction (Figure . ). The resulting model is therefore hybrid and features a gene expression layer where reactions are stochastic and a signal transduction layer where reactions are deterministic.

Inferring the rates of a stochastic protein turnover model is difficult

This approach introduces new parameters ( rate parameters per protein) that need to be inferred. While the rates of such stochastic protein turnover models are rarely directly mea-

E[G] = k on k on + k off = k on γ g (γ g = k on + k off ) ( . ) E[m] = E[G] k sm γ m ( .
)

E[P ] = E[m] k sp γ p ( . )
Coefficients of variation (standard deviation divided by mean, abbreviated CV from now on) of the gene status, mRNA level and protein level

CV [G] 2 = 1 -E[G] E[G] ( . ) CV [m] 2 = γ m γ m + γ g CV [G] 2 + 1 E[m]
( . )

CV [P ] 2 = γ m γ p (γ g + γ m + γ p ) (γ g + γ m )(γ g + γ p )(γ m + γ p ) CV [G] 2 + γ p γ p + γ m 1 E[m] ( . )
Usually, only the steady-state distribution of protein levels is experimentally available, for example using specific antibodies and flow cytometry. It has been shown that the moments of this distribution, and even its shape, is usually not sufficient to identify a single set of parameters of the stochastic protein turnover model that generates this distribution at steady-state (Sherman & Cohen ). This unidentifiability cannot be fully resolved by the additional knowledge of the protein degradation rate γ p .

To infer the gene activity bursting kinetics with precision, another approach requiring genetic engineering and live-cell microscopy has been proposed (Suter et al.

). It relies on modifying the sequence of the gene to obtain both a short-lived mRNA and protein, the la er having luciferase activity enabling its quantification in live cells. Using maximum likelihood methods on obtained single-cell trajectories, estimates for k on and k off can be recovered. The amount of experimental work needed by this method, together with the fact that additional estimations of the endogeneous mRNA and protein stabilities are required to describe fluctuations of the endogeneous protein, renders it in practice intractable for the systematic inference of protein fluctuation models to integrate in large size signal transduction pathways.

Two statistics to characterize protein fluctuations

The difficulty to infer all the rates of a stochastic protein turnover model instance could mean that their full knowledge is not required to describe the fluctuations of protein level with a satisfying accuracy. To investigate this idea, it is needed to define a small number of measures or statistics that provide a characterization of the fluctuations generated by a stochastic protein turnover model instance.

A first natural choice is the relative wideness (CV) of the protein level state distribution.

Its expression was given in ( .

Approximating the fluctuations of long-lived proteins by standard models

Intuitively, fluctuations of short-lived proteins are expected to be more sensitive to the precise kinetics of bursting than long-lived proteins, because the protein level at a given time reflects previous protein synthesis along a time window related to the protein half-life. Therefore, the level of long-lived proteins only informs about a cumulated amount of protein synthesis events, but not about when exactly those events occurred. On the other hand, the level of short-lived proteins is more closely related to the instantaneous rate of protein synthesis. It is this argument that led (Suter et al.

) to use short-lived mRNAs and proteins to infer the transcriptional bursting kinetics associated to a promoter sequence.

To investigate the validity of this hypothesis, we explored the relationship between the two statistics CV and τ and the bursting parameters T on = 1 k off and T off = 1 k on in two different se ings of protein and mRNA stability (Figure . ). In the first se ing, the protein is stable (its half-life is hours and is set by dilution), and the mRNA half-life is hours. It corresponds to typical stabilities found in mammalian cells (Schwanhäusser et al.

). In the second se ing, both the protein and mRNA are short-lived (half-lives of and hours respectively). Imposing k on , k off , γ m and γ p is not sufficient to determine the fluctuations of a stochastic protein turnover model: k sm remains unconstrained and can impact CV and τ . It was chosen such that the mean mRNA level E[m] equals , corresponding to the median value for the mRNAs quantified in a genome-wide study in mammalian cells (Schwanhäusser et al.

).

The range of both CV and τ are narrower for the first se ing compared to the second, pro- viding support to the hypothesis that protein fluctuations are less sensitive to the bursting rates when the protein is long-lived. Otherwise stated, for typical protein and mRNA halflives, a large set of promoter rate combinations leads to fluctuations of similar amplitude and speed as characterized by CV and τ .

Interestingly, in that case the obtained mixing time is around hours, in the middle of the range of experimentally estimated values for twenty endogenous proteins in human cells (Sigal et al. ). This non-trivial consistency result adds to the relevance of using the stochastic protein turnover model to describe protein fluctuations in human cells.

Therefore, using a standard parameterization for the stochastic protein turnover model can provide a good approximation of protein fluctuations for most proteins. By standard parameterization we mean a parameterization that respect typical constraints on protein and mRNA stabilities, mean mRNA level and promoter switching times, as will be illustrated concretely on the example of TRAIL-induced apoptosis. Fluctuations and steady-state distribution of such a model are presented in Figure . . This result is a cornerstone of our approach: it enables a parsimonious parameterization strategy in which only short-lived proteins necessitate particular a ention, while others are equipped with standard stochastic protein turnover models.

Application to TRAIL-induced apoptosis

We applied our approach to TRAIL-induced apoptosis. We used the EARM kinetic model mentionned earlier (Spencer et al. ) promoter switching rates, mRNA half-life and protein half-life (Figure . ).

Importantly, all parameters have been constrained based on experimental data and analytical results, with the exception of four parameters ("ON" and "OFF" promoter switching rates, mRNA and protein half-lives for FLIP and Mcl-). Note that because of their similar protein and mRNA half-life, we used the same couple of promoter switching rates for both proteins in our exploration of parameter space. Thus, our parsimonious parameterization strategy for the stochastic protein turnover models drastically limited the number of introduced degrees of freedom. This enabled us to systematically explore realistic ranges for remaining parameters.

Aside of the parameterization of the introduced protein fluctuation models, we made several modifications to the parameters of the model used in (Spencer et al.

). First, the half-life of the death ligand TRAIL, that was originally set to ~ hours, was set to hours, following an available experimental measurement (Youn et al.

).

Another discrepancy between the original parameters and biological knowledge was the degradation rates of the non-native forms. Almost all of them were assumed to also have a half-life of ~ hours, while as we already mentionned, most activated apoptotic proteins are known to be targeted for proteasomal degradation and hence exhibit much shorter half- (Spencer et al., 2009) Figure 2.7: Building stochastic protein turnover models for TRAIL-induced apoptosis. Routine followed to choose rates of all 17 native proteins in the EARM kinetic model of TRAIL-induced apoptosis. Typical values from multi-genes studies in mammalian cells are used to constrain rate values. Specific attention is given to FLIP and Mcl-1 because they are known to be short-lived, and thus more prone to exhibit large variations.

Finally, we did not change the kinetic rates of the protein-protein reactions, except the rate of pC cleavage by C , set to zero to represent the assumption that the caspase feedback loop (the ability of caspase to cleave pro-caspase , which in turn can cleave pro-caspase ) is absent (see authors comment to (Albeck, Burke, Spencer, et al. ), on editor's website).

A complete list of parameter values that we used is given in Appendix .

The impact of those changes will be discussed in the discussion closing this Chapter.

Importantly, most results presented in this Chapter are robust with respect to those changes.

However, as we will see in the next Chapter, the long-term response of cell populations treated with TRAIL are strongly impacted by the rates at which activated apoptotic proteins are degraded.

To study the influence of stochastic protein turnover on fractional killing and reversible resistance, we sought to confront our model with existing quantitative data about TRAILinduced apoptosis in HeLa cells. Those experiments, described in detail later, can be classified into two groups based on the type of information they contain: ) quantification of the variability in cell fate, ) characterization of the transient memory in cell state. While previous approaches using ODE models with distributions for initial protein levels (capturing a static description of cell-to-cell variability) (Spencer et al.

; Gaudet et al. ) are potentially able to reproduce the first type of data, a dynamic view on cell-to-cell variability as proposed in our model is needed to account for both types of data. We adopted the following strategy: first, search for models able to reproduce observations on cell fate variability; and second test whether valid models can robustly predict observed behaviors where transient memory ma ers.

. We first asked if the observed cell fate variability could be reproduced. In the model, it is only determined by protein levels at treatment time (behavior is deterministic as synthesis is assumed to be fully blocked by CHX and noise in signaling reactions is neglected), and differences between sister cells are only caused by protein synthesis noise occurring between division and treatment (in agreement with the fact that recently divided sisters died almost synchronously, we assumed an equal repartition of protein content at division).

We found that excellent agreement with observed MOMP time variability can be ob- Also, cell fate inheritance between sister cells was markedly changed: only young sister cells that underwent MOMP rapidly were importantly correlated (Figure . -B, grey curve).

We asked whether the observed cell fate variability, including fractional killing, could be reproduced in-silico. Within our modeling assumptions, absence of co-treatment with CHX makes a fundamental difference: as synthesis continues, the effect of gene expression noise during TRAIL-induced apoptosis could be investigated, and comparison with Mcl-1 promoter switching rates are standard (Ton=0.1 hours and Toff=2.6 hours) but the short half-life of their mRNA and protein is accounted for (2 hours and 0.5 hours respectively). We reproduce here for this model all the results presented in Figures 2.8, 2.9, 2.16 and 2.17 ). This property is thus a model prediction. For the top frame, promoter switching rates are standard (as in Fig. 2). Because mRNA and protein half-lives are short, protein level fluctuates more rapidly and the steady-state distribution is changed (it is wider and the mode is in 0) compared to the standard stochastic protein turnover model (Fig. 2). On the bottom frame, the steadystate distribution becomes bimodal because the promoter switching rates are low compared to mRNA and protein degradation. In both cases, fluctuations and distribution are shown for Mcl-1; they are similar for FLIP as only the protein synthesis rate changes to account for a different mean protein level. exploration, but further analysis showed that assuming low promoter switching rates for

Mcl-alone was sufficient to obtain quantitative agreement for MOMP distributions surviving fractions, and that sister cells MOMP time correlations were still correctly predicted (Figure . ). Thus, comparison with transient cell fate inheritance data supports that large, rare fluctuations of Mcl-could be responsible for the observed cell fate variability. ). The switching rates of the FLIP promoter were assumed to be standard (Ton=0.1 hours and Toff=2.6 hours). All the results presented in the main text for the "fitted" model are reproduced here. Note that because fewer cells were simulated compared to main text figures (5.10 4 instead of 10 5 for sister cell experiments), sister correlation curves appears slightly noisier.

.

Model predicts reversible resistance

Recently, reversible resistance was observed among various cell lines (Flusberg et al. ).

Cell populations were submi ed to two consecutive TRAIL treatments. The duration between treatments was varied from day to week (Figure . -A). One-day survivors were significantly more resistant than the initial population, but such resistance was significantly decreased or even lost in one-week survivors. Thus, cells surviving a first TRAIL treatment are transiently resistant.

Remarkably, in-silico reproduction of those (Figure . -B) showed that our model predicts the presence of reversible resistance (Figure . -C): one-day survivors exhibit a dosedependent increase of resistance to a second TRAIL treatment, which disappears after to days. This is surprising since our model does not include induced regulation mediated by survival pathways. Moreover, the presence of reversible resistance is a robust property of the model as it is also obtained when assuming standard promoter switching rates for

Mcl-and FLIP ("non-fi ed" model, Figure . ).

However, agreement with related experimental data (Flusberg et .

Molecular determinants of fractional killing and reversible resistance

What are the mechanisms behind cell escape to TRAIL-induced apoptosis, either on the short-term (fractional killing) or the long-term (reversible resistance)? Using the fact that in-silico, all protein, mRNA levels and gene activity states can be monitored in single cells, we investigated those questions at the molecular level.

To study the influence of pre-existing differences on cell fate, we compared at the time of stimulation the sub-population of 'future survivors' with the whole population ( To investigate the determinants of reversible resistance, we tracked the temporal evolution of protein levels in surviving cells (Figure . -A,C). The protein level composition of one day survivors contrasts with the protein content observed in future survivors: almost all protein levels differ importantly from the naïve population composition, while that was the case only for Mcl-and FLIP in future survivors. This is expected as all proteins are partly activated during signal transduction, leading to a higher degradation (active forms have a shorter half-life).

Therefore, the distinction between the causes of cell survival and the consequences of partial apoptosis induction cannot be easily resolved by the sole observation of protein levels in survivors. When signaling stops, recovery of protein levels is expected to follow exponential kinetics governed by the turnover rate ( synchronously when treated with TRAIL and a protein synthesis inhibitor, while in contrast, unrelated cells died after highly variable durations. This demonstrated that TRAIL signaling is mostly deterministic when protein synthesis is blocked and that the timing of death is determined by the cell internal state at the time of treatment.

Second, they observed that such synchrony in sister cells death is gradually lost as the time between division and treatment increases. This showed that the cell 'TRAIL sensitivity state' (the part of cell internal state involved in death timing determination) naturally fluctuates over a dozen of hours. In addition, the modeling results in (Spencer et al.

)

highly suggested that such state is mainly composed by the various levels of the proteins acting in the extrinsic apoptosis pathway.

In parallel, important progress on the characterization of the stochasticity in gene expression has been made: the two-state transcriptional bursting model was shown to permit high accuracy and several approaches to infer its parameters were proposed, enabling the quan-

Modeling protein fluctuations in TRAIL-induced apoptosis

In this Chapter, we presented how we merged those two approaches by integrating such stochastic models of gene expression within an existing kinetic model of TRAIL-induced apoptosis (Albeck, Burke, Spencer, et al. ) in a systematic and principled manner. Doing so provides advantages compared to previous approaches to account for cell-to-cell variability in protein levels (Spencer et al.

; Gaudet et al. ; Kallenberger et al. ).

First, variability is not considered as an "input" parameter but arises naturally from stochastic fluctuations. The dynamics of this variability is thus intrinsically represented within the system, allowing investigating the effects of transient memory in protein levels. Second, the influence of protein synthesis noise during TRAIL-induced apoptosis could also be investigated. Importantly, we followed a parsimonious parameterization strategy, motivated by the fact that fluctuations of long-lived proteins are rather insensitive to the precise kinetics of transcriptional bursting, enabling us to equip most proteins (long-lived proteins) with reasonably accurate fluctuation models even in absence of gene expression data for each and every promoter.

The sister cells experiment for which cells were treated with TRAIL and CHX provided ideal data to validate our modeling approach: in that case, behavior is mostly deterministic as soon as treatment starts and only fluctuations occurring before treatment are responsible for death time variability and de-correlation between sister cells. Moreover, gene regulation via survival pathways induction is ineffective as protein synthesis is blocked.

Because our model was able to quantitatively reproduce the MOMP time distribution and then accurately predicted sister cells correlation, our modeling approach appears as a promising tool to investigate the effect of protein fluctuations on signal transduction, despite the limitations inherent to its simplicity (for example, stochastic gene expression events were assumed to be independent between proteins, neglecting the fact that levels of different proteins can be partially correlated (Gaudet et al. ; Sigal et al. ), possibly because of common transcription factors or coordinated chromatin-state transitions). Of note, good agreement was readily obtained when assuming for FLIP and Mcl-standard pro-moter switching rates (but short protein and mRNA half-lives, in agreement with available knowledge -"non-fi ed" model, Figure . ).

Finally, the transposition of fluctuation timescales from individual proteins into 'TRAIL sensitivity states' is not trivial: while in our model, stable proteins levels are mixed in about hours, cells were switching between 'fast dying' and 'slow dying' phenotypes more rapidly (about -hours). As combinatorial and non-linear effects are at play, mechanistic models of protein-proteins reactions are needed to link protein-level timescales with more high-level phenotypic transitions (Gupta et al. ).

Questioning the role of survival pathways

Several studies reported that TRAIL can induce survival pathways (Chaudhary et al. ;

Son et al.

; Sun et al.

). How such induced changes affect signal transduction and eventually stop apoptotic signaling remains unclear. On the other hand, the contribution of constitutive protein synthesis noise, which is responsible for pre-existing differences between cells, has not been evaluated.

Although it does not exclude the existence of other mechanisms, an important result of our study is that fractional killing can be obtained without assuming any TRAIL-induced regulation. Alternatively, we find that because of its fast turnover, constitutive expression ). Improving how DISC assembly is modeled might thus be needed to elucidate the precise role of FLIP in fractional killing and reversible resistance, especially for cell lines that express higher FLIP amounts than HeLa.

Origins of reversible resistance: joint effect of selection and stochastic protein turnover

A second significant result reported here is that our model predicts the phenomenon of reversible resistance, showing that constitutively noisy protein synthesis, protein-protein interactions and protein degradation are by themselves sufficient to explain a dose-dependent, significant increase of resistance in recent survivors and its gradual loss within -days.

This result is consistent with the observation that NF-κB blockade does not change resistance acquisition after TRAIL treatment (Flusberg et al. ) (in MCF A cells; HeLa cells have not been tested). In-silico analysis at the molecular level revealed that reversible resistance as predicted by the model was shaped by a complex interplay between ) selection based on protein levels and transcriptional activity, ) protein turnover and ) residual signaling activity.

As opposed to the death process, which involves a sharp and complete activation of effector caspases, our results suggest that recovery in cells that did not commit to death is a slow and complex process. While one should not conclude from our results that parallel activation of survival pathways by TRAIL plays no role in reversible resistance, our results show that the sole contribution of protein level fluctuations occurring within the extrinsic apoptosis pathway can partly lead to reversible resistance. Thus, protein fluctuations should be accounted for to gain quantitative insights into reversible resistance.

Impact of caspase feedback loop and non-native forms degradation on model behavior

As previously justified in the section presenting model construction, we assumed in all our simulations ) absence of the caspase feedback loop (i.e., the ability for caspase to cleave pro-caspase , which can then activate pro-caspase ) and ) non-native forms for which a degradation rate could not be specifically a ributed (all non-native forms except TRAIL, the species representing mitochondrial pores and the complexes involving FLIP and Mcl-) are degraded with a hours half-life.

We investigated the impact of those choices on our main results by repeating the corresponding in-silico experiments for different values for the non-native forms half-life and status of the caspase feedback loop. The results of the short-term, single-cell experiments are almost un-affected when the non-native forms are longer-lived ( hours) and similar when they are shorter-lived ( hours) (Figure . ). Here, the presence or absence of the caspase feedback loop had no effects.

Regarding model predictions for transient resistance acquisition, the quantified resistance gain in one-day survivors is not strongly affected. Still, the presence of the caspase feedback loop systematically increases it, while it decreases when non-native forms half-life is below hours (Figure . ). Note that here resistance acquisition was quantified by comparing cell numbers at treatment time and hours after (to follow Flusberg and colleagues) and hence does not necessarily reflect the total killing efficacy of a TRAIL treatment. In fact, we found the long-term evolution of population size for cells treated with TRAIL to strongly depend on non-native forms degradation. The source of this dependence and its potential consequences for the outcome of repeated TRAIL treatments are investigated in details in Chapter .

Accounting for gene expression noise appears necessary to investigate signal transduction

While here we focused on TRAIL-induced apoptosis, our modeling approach is generic and can be applied to other signal transduction pathways. Our results showed that even in absence of induced gene regulation, gene expression noise interacts with signaling dynamics on a non-trivial manner. Thus, even in contexts where the influence of induced generegulation is indisputable, its sound quantification probably requires to investigate first the role of constitutive gene expression noise.

Only then models could be enriched parsimoniously with well-characterized regulatory links until all observations are successfully explained. Significant advances to allow such detailed characterization of gene regulation occurred recently (Tay et al.

; Molina et al.

; Neuert et al.

). Following such approaches could significantly extend the reach of models of signal transduction towards accurate, single-cell level description of populations submi ed to varying signaling contexts. On the methodological side, we illustrate in Chapter how to extend our modeling approach to account for regulated stochastic gene expression and demonstrate that it allows to capture induction dynamics of an artificial yeast signal transduction pathway with single-cell resolution. Gaining a quantitative understanding of how those different factors shape TRAIL resistance in-vivo would permit to rationally design treatment strategies to maximize their effi- . The role of time: could resistance acquisition protect tumor monolayers from repeated TRAIL treatments?

Modeling populations of single cells submi ed to TRAIL treatments

In this section, we investigate the long-term consequences of fractional killing and transient resistance acquisition as predicted by our model. To allow the simulation of the longterm evolution of cell populations submi ed to arbitrary TRAIL stimulations, the single-cell TRAIL model is coupled with a simple model of cell proliferation (Figure . ). This model accounts for cell-to-cell variability in the cell cycle distribution in order to reproduce the loss of synchrony in cell divisions that is observed in cancer cell micro-colonies cultured in-vitro.

Cells of the population are assumed to experience the same concentration of TRAIL in their environment. This unique environmental TRAIL concentration can be time-varying and is assumed not to be influenced by cells (Figure . ). This situation is meant to represent the long-term culture of cancer cell populations as monolayers.

Single-cell model of TRAIL-induced apoptosis Cell proliferation model Population of individual cells

TRAIL apoptosis signaling

A reference model for HeLa cells

As a basis for our analysis, we consider one of the most parsimonious parameterization of ). The non-native forms are assumed to be degraded with a half-life of hours when no specific values can be set. The caspase feedback loop is assumed to be absent. Finally, we brought a minor modification to the previously tested parameterizations: cPARP (which value provides a readout for the integrated amount of caspase-activity and is used for determining irreversible commitment to apoptosis) is assumed to be stable (hence 'degraded' only via dilution).

In the previous Chapter, model-data agreement in terms of cell killing was tested for a single TRAIL dose (in absence of CHX co-treatment) and for short-term data (up to hours after treatment). Because we are now interested in the long-term cytotoxic effect of arbitrary TRAIL dose temporal profiles, we first tested the ability of this model parameterization to reproduce existing data that spans a wide range of TRAIL doses and in which killing efficiency is measured after hours of treatment (Lemke, Karstedt, Abd El Hay, et al. ).

We obtained a reasonable agreement for both wild-type and Mcl-shRNA expressing HeLa cells dose-response curves (Figure . ). Thus, the predicted killing efficiencies for HeLa

In-silico quantification of initial resistance: reproducing a standard doseresponse curve

The standard way to evaluate the potential of drugs as anti-cancer cytotoxic agents is to perform dose-response experiments on cancer cell lines cultured in-vitro. We performed such experiment in-silico by proceeding as follows. First, we started with small populations of single cells that were let grown without TRAIL until they reached a size of cells. Then, we applied different doses of TRAIL to those populations and simulated their evolution for days.

Figure . shows the resulting curves for the temporal evolution of population size. For each dose, two independent replicates were simulated, to provide indication on the variability inherent to the finite size of our population. Here, this variability exists but is rather limited (for the initial population size we used). As expected, low doses only moderately affect cell growth, while higher doses are able to decrease population size. Both the extent and duration of cell killing increases with TRAIL. From those curves, one can extract a standard dose-response curve, which is the common output of standard cytotoxicity assays (Fallahi-

Sichani et al.

). A typical end-point duration for anti-cancer drug cytotoxicity assays is days. The simulated dose-response curve corresponding to that end-point is also shown Cell number TRAIL=0.05ng/ml TRAIL=0.1ng/ml TRAIL=0.5ng/ml TRAIL=1.0ng/ml TRAIL=5.0ng/ml TRAIL=10.0ng/ml TRAIL=50.0ng/ml TRAIL=100.0ng/ml TRAIL=500.0ng/ml TRAIL=1000.0ng/ml TRAIL=5000.0ng/ml in Figure . . What informations on the ability of TRAIL to eradicate or at least control a population of growing HeLa cells is provided by such dose-response curves? Here, it suggests that there exists a dose (between and ng/mL) such that treating every hours would be just sufficient to maintain the population size below its initial size (as could also be seen from the corresponding cell number curves). Moreover, higher doses would apparently be sufficient to rapidly eradicate the population, as a -fold reduction of cell number could be achieved in days.

Dose-response cell number curves Population of individual cells

Simulate for 72 hours

Apply different TRAIL concentrations

Long-term resistance to repeated TRAIL treatments cannot be predicted from the response to single treatments

The previous considerations ignored the possibility that the population of cells surviving a first TRAIL treatment can be transiently more resistant than the initial population. But Resistance acquisition between the first treatment and the subsequent treatments is clearly visible. For example, while a dose of ng/mL resulted in more than a -fold reduction in cell number days after the first treatment, all the following treatments killed cells but not enough to balance growth, causing a net growth in population size on the long-term. Resistance does not accumulate at each treatment: rather, it seems to stabilize after the first one or two treatments. Thus, it is possible to define a "long-term" dose response curve, by computing for each dose the net growth rate between the second and last pulse. This net growth rate can be translated into a ratio of cell number difference with respect to a control for a days time period, allowing comparison with the "standard" dose-response curve (Figure . ).

This comparison highlights the potential therapeutical consequences of non-genetic resistance acquisition as predicted by our model: the dose that is needed to "control" the population size is times higher that what is predicted from the standard cytotoxicity dose-response curve. The maximum killing efficiency (E inf with the notations of (Fallahi-Sichani et al.

)) is also markedly changed. Therefore, our results suggest that when choosing the dosage of an anti-cancer therapeutic strategy, relying solely on standard, single-treatment Cell number TRAIL=0.05ng/ml TRAIL=0.1ng/ml TRAIL=0.5ng/ml TRAIL=1.0ng/ml TRAIL=5.0ng/ml TRAIL=10.0ng/ml TRAIL=50.0ng/ml TRAIL=100.0ng/ml TRAIL=500.0ng/ml TRAIL=1000.0ng/ml TRAIL=5000.0ng/ml !"#$#%&'(")&(*+$#$,-') Extraction of a "long-term" dose-response curve from the long-term repeated TRAIL treatments simulations (green). The "standard" dose-response is shown in red for comparison. The long-term dose response curve does not extend to low TRAIL values because in those cases population sizes grew too fast to allow an accurate computation of the long-term net growth rate with reasonable computational times. The TRAIL dose that would allow to balance proliferation on the long-term is highlighted and compared to its equivalent as predicted by the standard dose response curve.

Influence of scheduling on the outcome of long-term repeated TRAIL treatments: despite transient resistance acquisition, waiting between treatments is not favorable

In the previous analysis, we considered a fixed period of days between treatments. It is a priori unclear whether increasing or decreasing this period could improve the long-term outcome of the repeated treatments. On one hand, if the period is too long, after a treatment, killing might stop to be important before the next treatment (hence cell proliferation would dominate in the interval). On the other hand, as suggested recently (Flusberg et al. ), because of transient resistance acquisition, increasing this period could enable to recover more sensitivity and thus more killing efficiency for the next treatment.

To study the influence of scheduling, we took the treatment of ng/mL every days as a reference (yellow curves in Figure . ) and tested different periods between treatments while keeping total TRAIL dose delivered constant (Figure . ). First, while for unfrequent treatments, the net growth rate of the population displays large variations in between treatments (strong killing followed by dominant proliferation), the amplitude of those variations decreases with the treatment period. Second, the efficiency of the repeated TRAIL treatments to limit population growth increases with treatment frequency until reaching a plateau: the hours or hours period treatment strategies were equally good.

Hence, our model predicts that the trade-off between waiting for the initial sensitivity to recover and the cell proliferation that occurs during this waiting time is such that waiting is not favored but that maintaining sustained (hence lower) killing is preferable. This important model prediction can be tested experimentally.

Distinguishing population-level or cell-level resistance acquisition

We found that non-genetic resistance acquisition as predicted by our model could have a The second type of mechanism involved in resistance acquisition operates at the population level. Even if the individual cells "TRAIL resistance" is not changed in surviving cells, cells with low surviving probability are by definition eliminated in higher proportions than cells with high surviving probabilities. Hence, the amount of cells with high surviving probabilities is expected to increase by a pure selection effect. In that case, we speak of population-level resistance acquisition as it does not involve changes of resistance of individual cells.

In our model, transcriptional noise is assumed not to be affected by TRAIL. Therefore, changes in the average mRNA levels in a population necessarily imply a selection effect. In Influence of scheduling on the long-term outcome of repeated TRAIL treatments. Different doses are applied at various frequencies such that in average, a dose of 100 ng/mL is delivered every 3 days. In other words, the same total amount of TRAIL is delivered either at strong doses un-frequently or at low doses frequently. Top: temporal evolution of alive cell number. Bottom: temporal evolution of TRAIL concentration. In each case, two independent replicates are shown. The inset in the top plot shows the endpoint cell number as a function of treatment period (marks indicate simulations that were stopped earlier because cell number was too high). Note that for treatment periods below 48 hours, the long-term killing efficiency is approximately constant.

are shown) being transiently enriched after TRAIL treatments while transcripts coding for pro-apoptotic proteins (pC and Bid are shown) are transiently depleted after treatments.

However, the quantitative impact on average mRNA levels is small ( % at maximum) and is in general almost completely erased hours after each treatment. Still, it cannot yet be excluded that the combined effect of those small differences can have an impact on the population resistance, as they all go in the same direction (enrichment/depletion of proapoptotic/anti-apoptotic mRNAs).

As opposed to mRNAs, the levels of native proteins can change because of cell-level effects caused by TRAIL signaling. For example, pro-apoptotic proteins can be cleaved and then degraded, causing a decrease in the level of the native, uncleaved protein. Thus, the temporal evolution of native protein level population averages reflects both selection and cell-level effects. Those evolutions are shown in Figure . . Almost all protein levels decrease after treatment before recovering, even for anti-apoptotic proteins such as XIAP. The example of XIAP is interesting because mRNA levels are enriched (Figure . , bo om-right)

while protein levels are depleted, demonstrating that cell-level decrease via degradation is dominant over population-level selection, at least regarding XIAP evolution.

Overall, those results suggest that protein degradation in surviving cells is a main driver of TRAIL resistance acquisition as predicted by our model. A heuristic explanation based on the model structure for our observations is that in order to survive, cells need to limit the transformation of death signals into effector caspase activation; and this is achieved by binding of anti-apoptotic proteins to activated pro-apoptotic proteins followed by targeted degradation, resulting in decreased amounts of both pro-and anti-apoptotic proteins. However, because some anti-apoptotic proteins exhibit fast turnover, the balance between provs anti-apoptotic proteins is transiently in favor of anti-apoptotic proteins, resulting in an increased resistance to subsequent treatments.

To rigorously demonstrate the dominant role of TRAIL signaling related degradation of proteins in resistance acquisition, we ran a long-term, repeated TRAIL treatments simulation in which its influence was artificially removed (Figure . , green curve). More precisely, before applying the second (or third, etc…) TRAIL treatment, we reset the state of each (alive) cell to what it could have been if no TRAIL had been added at the first (or second, etc…) TRAIL treatment. By state we mean the complete state of the intra-cellular model, i.e. the levels of protein species (native or non-native) and the level of mRNAs and promoter activity status for the native proteins. Hence, only selection effects are kept.

As expected, although it is still visible, long-term resistance acquisition is importantly decreased compared to the normal simulation (red curve). 3.4 was performed for the dose 100 ng/mL (red curve), except that states of cells alive at the i+1 th treatment were reset to what they could have been if no TRAIL had been added at the i th treatment (green curve). Horizontal thin black lines indicate times of treatments. The dashed black line highlights the evolution of cell number at treatment times that is expected in absence of resistance acquisition.

Key role of targeted protein degradation in determining both initial and long-term resistance

The previous results were obtained with a single model parameterization. To investigate further the role of targeted degradation in TRAIL resistance, we varied the aggregate parameter in our model that used for the half-life of most non-native forms (all except Flip:R* & Mcl-:tBid, for which we use the degradation of the non-native form; the abstract species representing mitochondrial pores; and cPARP, which is our readout for effector caspases activity).

In each case, we simulated the long-term dose-response as done previously for our reference model parameterization. The resulting cell number curves for the highest dose ( ng/mL) are shown in Figure . (upper plot). The rate of non-native forms degradation has a strong impact on both the initial and long-term resistance to TRAIL: compare for examples the curves corresponding to a half-life of . and . hours. The resistance to TRAIL increases with the degradation of non-native forms.

From those curves, it seems that this sensitization effect applies on both the initial and long-term resistance. To analyze further the relationship between the initial and long-term resistance, we computed in each case the initial killing efficiency (i.e., the effective killing rate obtained between the first and the second treatment) and the long-term killing efficiency (i.e., the effective killing rate obtained between the second and last treatment). Suprisingly, we found that the ratio between those two quantities was independent of the non-native forms degradation rate and hence of the initial resistance (Figure . , bo om plots). In summary, our model predicts that altering the rate at which non-native forms are degraded can strongly decrease TRAIL resistance, and that initial and long-term resistance are equally affected.

The most direct way to alter those rates is to treat cells with proteasome inhibitors such as bortezomib (Velcade as commercial name, this drug being already approved as an anticancer drug). Indeed, it is now recognized that TRAIL resistant cell lines can be strongly In all plots, different colors correspond to different values for the aggregate parameter describing the degradation of non-native forms. Top: cell number curves for 1000 ng/mL TRAIL treatments repeated every 3 days. In each case, two independent replicates are shown. From those curves, the initial and longterm killing efficiencies can be computed. Bottom left: Long-term versus initial killing efficiency. Bottom right: ratio between the initial and long-term killing efficiencies ("resistance factor") versus the initial killing efficiency. The grey dashed lines highlight the apparent linear relationship between initial and long-term killing efficiencies.

A -parameters phenomenological model captures the observed dynamics of resistance under TRAIL stimulation

The previous results show that the predicted resistance acquisition when TRAIL treatments are repeated is not driven by targeted degradation, as altering it equally impacts initial and long-term resistance. But then, what parameters are implicated in resistance acquisition?

We asked whether altering the balance between anti-and pro-apoptotic protein expression levels could impact resistance acquisition. To answer this question, we repeated the previous analysis except that we either ) changed the mean expression levels of all anti-apoptotic proteins while keeping the mean expression levels of pro-apoptotic proteins unchanged or ) the opposite. Again, in both cases, those alterations were able to impact initial and longterm resistance, but not their ratio (not shown).

Next, we tested the impact of the 'dilution' rate, i.e. the rates at which the stable native proteins (all but Flip and Mcl-) are degraded. Interestingly, we observed the opposite behavior compared to previous variations of non-native forms degradation or protein expression levels: resistance acquisition is strongly affected (it increases as dilution rate decreases), while initial resistance is almost un-affected (Figure . ). 3.10 except that the turnover rate of stable proteins is varied (i.e., dilution rate). The fold change corresponding to the reference model is indicated in the colorbar (1 corresponds to a half-life of 27 hours, 10 to a half-life of 2.7 hours, etc.). In each case, the corresponding protein synthesis rates are scaled by the same factor such that the mean expression level is unchanged. Grey dashed lines indicate the linear relationship between initial and long-term killing efficiency that was valid for variation in non-native forms degradation (see Figure 3.10). It does not hold here, instead, resistance acquisition is strongly affected while initial resistance does not change much.

The fact that some types of parameter modifications can change the overall TRAIL resistance without changing resistance acquisition itself while others can change resistance acquisition without affecting much initial resistance suggests that the model behavior can be described with a "factorized", low-parametric phenomenological model in which the separate contributions of initial (or 'baseline') resistance and resistance acquisition could be easily identified. To investigate this hypothesis, we formalized our heuristic explanation of the model behavior into a simple, -parameters model (depicted in Figure . ). This model features a single, population-level variable G describing the instantaneous ability of the population to transform death signals into actual killing (transmission capacity).

For simplicity, we use a binary representation of the death signal (DS): it is either present (DS = 1) or absent (DS = 0). The death rate of the population is then given by λ(t) =

G(t) DS(t): G(t

) is simply the 'potential' death rate that is 'realized' if death signals are present at t. The equation describing the dynamics of G is the following:

dG dt = σ -γ G(t) -k DS(t) G(t)
(G = death rate obtained in presence of death signals)

In this equation, σ and γ represent the normal synthesis and degradation rates of the transmission capacity G. The additional term k DS(t) G(t) accounts for the fact that transforming death signals into killing has a cost in terms of transmission capacity (for example because cells with low transmission capacity will be positively selected or because degradation will decrease the balance of pro-apoptotic proteins over anti-apoptotic proteins). In summary, the killing dynamics in cell populations facing TRAIL stimulation is described with only three parameters: σ, γ and k.

This equation is easy to solve analytically for any temporal pa ern of presence/absence of death signals (i.e. the treatment strategy). In absence of death signals, G converges exponentially to G off = σ γ , and hence the population would have the maximum response (i.e., death rate λ off = G off ) when death signals are added. On the contrary, in presence of death signals the transmission capacity converges exponentially towards G on = σ γ+k ≤ G off . Thus, the long-term death rate λ off = G off in presence of death signals is lower than the initial death rate λ off obtained when a naïve population is treated.

A prediction of this phenomenological model is that the "resistance factor" (i.e. the relative change of killing efficiency of subsequent treatments compared to the initial treatment, as shown in Figure . , bo om right) would converge exponentially to as the period between treatments is increased. Interestingly, this behavior is also observed for the original, mechanistic model (Figure . ). For the mechanistic model, repeated treatments with 1000 ng/mL TRAIL were performed for different treatment periods. In each case, the resistance ratio (killing efficiency of late treatments divided by killing efficiency of the first treatment, where killing efficiency is the ratio between cell number at the time of next treatment and cell number at treatment time) was computed. Parameters of the phenomenological model were adjusted to reproduce the same dependency, assuming that each treatment leads to 3 hours of death signal presence. Those parameters are σ = 1, k = 2 and γ = ln 2 55 hrs -1 . For the mechanistic model, increase of the resistance factor at low TRAIL periods is expected, as delays between TRAIL exposure and cell death result in death caused by the first treatment being attributed to the second treatment, etc.

Another prediction of the phenomenological model is that on the long-term, sustaining death signals is the best strategy to limit cell proliferation (this can be demonstrated analytically). Again, this is consistent with the original model predictions (see Figure . ). Hence, our simple model can capture essential features of the original model behavior regarding resistance to long-term repeated treatments.

In principle, it should be possible to map each parameter in the original, mechanistic model into its contribution to the phenomenological parameters σ, γ and k. It seems intuitive that altering protein synthesis rates would impact σ only, while changing the turnover rate of native proteins should impact γ mostly, and those intuitions are consistent with our observations. Interpreting the impact of non-native forms degradation in terms of the phenomenological parameters σ, γ and k seems less immediate and would require further investigation.

. The role of space: could geometrical effects protect tumor spheroids from repeated TRAIL treatments?

In the previous section, we investigated the role of time in the resistance of cancer cell populations treated with TRAIL by performing long-term single-cell resolved population dynamics simulations. In this section, we extend the model to investigate the role of cell-extrinsic heterogeneity that arises from the -dimensional organization of cells growing as spheroids.

We focus on the effect of limited TRAIL penetration within the spheroid that is expected if its diffusion is slow and its degradation is high.

Modeling cell growth and apoptosis for multi-cellular spheroids

For this purpose, it is needed to explicit the spatial localization of cells as well as the spatial distribution of TRAIL. This leads to a cell-based, D multi-scale model composed of several components highlighted in Figure . . To model adhesion/repulsion between cells in contact, we use the modified Her -model as in (Galle et al. ). In this model, the adhesion/repulsion force for a given contact The diffusion and degradation of TRAIL is also solved. Diffusion is assumed to be small when cell density is high, degradation is assumed to be high when cell density is high. In each cell, an instance of the singlecell model of TRAIL-induced apoptosis is simulated and governs cell decision to enter apoptotis. Apoptotic cells have their radius decreasing at a constant rate.

between two cells i and j derives from the following potential energy:

V i,j = 2 5 Ẽij ( R i + R j -d ij ) 5/2 R i R j R i + R j -σA ij ( . )
where

Ẽij =   3 4   1 -ν 2 i E i + 1 -ν 2 j E j     -1 and A ij = π R i R j R i + R j ( R i + R j -d ij )
In this equation, the first term corresponds to elastic repulsion. R i denotes the radius of cell i, d ij the distance between the two cell centers, E i the Young modulus of cell i and ν i its Poisson ratio. Ẽij is an effective Young modulus for this particular contact. The Young modulus has the dimension of a pressure and measures cell elasticity, while the Poisson ratio is a dimensionless parameter smaller than 0.5 measuring volume conservation upon deformation (i.e. compressibility). The second term corresponds to adhesion between the two cells. σ is the adhesivity (dimension of an energy per unit surface) and A ij the contact area between the two spheres (not the sphere overlap surface but an approximation of the true contact surface valid for small deformations (Landau et al.

)). The corresponding force is directed along the axis joining the two cell centers and has the following magnitude:

F adh/rep ij = Ẽij ( R i + R j -d ij ) 3/2 R i R j R i + R j -πσ R i R j R i + R j ( . )
We model cell micro-motility by a force which direction and magnitude stochastically switches at a rate 1 τ M M such that the mean duration between micro-motility force changes for a given cell is τ M M (the persistence time). The direction of this force chosen randomly (isotropic) and its magnitude is sampled from a gaussian distribution. Thus, parameters describe the micro-motility model: the persistence time, the average micro-motility force, and its standard deviation representing cell-to-cell variability.

Finally, each cell motion is slowed down by a friction force that follows the expression:

-→ F f riction i = -4πR 2 i µ -→ v i ( . )
Here for simplicity we have aggregated cell-cell/cell-ECM, normal/parallel friction coef-ficients into a single friction coefficient c i and used the functional form c = µ S c (Galle et al. ) where µ is the friction constant and S c the surface of contact.

It is generally accepted that inertia terms are negligible in comparison to friction terms (Galle et al. ). Indeed, given that cell mass m ∼ 10 -12 kg, the time constant with which velocity adapts to changes in forces other than friction is m 4πR 2 i µ ∼ 0.001 µ -1 s (for a cell radius of 10 µm), that is less than 10 -11 s for the value of µ used by (Galle et al.

). Therefore, the equation of motion for cell i can be wri en as:

4πR 2 i µ -→ v i = ∑ j,i in contact -→ F adh/rep j→i + -→ F motility i ( . )
For a given configuration of cells, this equation is straight-forward to solve as soon as cell-cell contacts have been detected. It provides a velocity vector for each cell that is used to compute its displacement for a given time step. In our simulations, the time step is adaptively chosen based on computed velocities such that no cell moves more than 0.1 µm, ensuring a certain accuracy and that no contacts are missed. All the parameters needed to simulate the motion of cells are given in Table . To model cell growth and division, as in the previous section each cell is a ributed a cell cycle duration sampled from a normal distribution. During the cell cycle, cell radius is assumed to increase linearly with a rate ensuring that volume doubles. At division, the mother cell is replaced by two half-volume daughter cells that are just in contact. The axis of division is assumed to be randomly isotropic. Parameters needed to simulate cell growth and division are summarized in Table . To model TRAIL penetration, we discretize a cuboïdal domain surrounding the spheroid with a regular grid, in which TRAIL diffusion and degradation is computed using an explicit first order scheme. To account for the fact that TRAIL diffusion (and degradation) is expected to be slower (higher) in regions of high cell density, we a ribute at each step a different diffusion coefficient ( and TRAIL half-life) dependending on the presence (D cells and HL cells T RAIL ) or absence (D medium and HL medium T RAIL ) of cells in that voxel. Because our aim is to study the potential impact of insufficient TRAIL penetration on the outcome of repeated TRAIL treatments, we made conservative assumptions for D cells and

HL cells

T RAIL values to ensure the existence of gradients within compact spheroids of more than cells. For comparison, the diffusion coefficient of water in water is 2000 µm 2 .s -1 and a typical diffusion coefficient for an antibody in a tissue is 10 µm 2 .s -1 (Thurber et al.

).

TRAIL treatments are modeled by using time-dependent Dirichlet boundary conditions that follows the same evolution as for 'monolayer' simulations (exponential decrease with a hours half-life in-between administrations). The parameters needed to simulate TRAIL penetration are summarized in Table .

Boundaries dimension

L medium µm

Finally, our single-cell model of TRAIL-induced apoptosis is simulated in each cell as previously except that the TRAIL concentration depends on the cell position. When a cell commits to apoptosis (i.e., its internal amount of cleaved PARP reaches a threshold), we assume that its radius is reduced at a given rate representing the speed of the cell lysis process (when a dying cell radius goes below a threshold, it is simply removed from the simulation). This parameter has potentially an impact on the spatial re-organization of the spheroid when cell death is high, which in turn could affect TRAIL penetration and hence subsequent killing.

Simulation of repeated treatments on spheroids

In the previous section, we found that because of resistance acquisition as predicted by our reference model of TRAIL-induced apoptosis in HeLa cells the dose required to control the size of a cell population (when repeated every days) growing as a monolayer was ng/mL, while ng/mL would have appeared sufficient from standard killing assay.

We asked whether insufficient TRAIL penetration could also impact the long-term outcome of repeated treatments. For this purpose, we ran several simulations of the multi-scale model described above as follows. We start by seeding a small population ( cells) within a small volume until it reaches a size of cells. Of note, this allows to obtain a realistic correlation between internal cell state (variable because of protein fluctuations) and spatial position: neighbors cells are more likely to be close siblings and hence to have correlated cell state.

We then apply TRAIL every days for more than days. A example visualization of the resulting spheroid evolution is shown in Different doses between and ng/mL were tested. In addition, two values for the radius reduction rate for apoptotic cells were tested ( and µm per hour). First, for ng/mL the long-term evolution of population size of the spheroids was undistinguishable from the corresponding 'monolayer' simulation (Figure . , top). Thus, treatments that were sufficient to control the proliferation of monolayers are predicted to also be sufficient for spheroids. However, for slightly lower doses significant differences of the long-term evolution are seen (Figure . , bo om). Importantly, in that case the speed at which apoptotic cells are lysed has an influence. Those results suggest that if only geometrical effects are at play, the key factor to consider for the long-term treatment of spheroids is resistance acquisition, which can be assessed on monolayer cultures. Thus, if experimentally strong differences between monolayers and spheroids are observed, it would point out to other effects, for example ECM-mediated changes in gene expression. Our modeling framework can easily be used to study the impact of such potential effects. In addition, because the model is detailed at the molecular level, the insights it can provide are not necessary limited to the HeLa cell line. Indeed, studying the relationship between its parameters and the predicted behavior could improve our understanding of the molecular mechanisms that shape TRAIL resistance for other cell lines, provided those mechanisms rely on species and interactions that are represented in the model.

In this chapter, we have developped a multi-scale model to investigate TRAIL resistance in virtual experimental se ings that could bridge the gap between standard in-vitro se ings and in-vivo se ings. The first se ing corresponds to the long-term (several weeks) response of cancer cell populations grown in-vitro in standard conditions (e.g. monolayers) but repeatedly submi ed to TRAIL treatments. It is meant to bridge the temporal dimension of the gap that exists between standard cytotoxicity assays and in-vivo experiments; otherwise stated, it allows to study the contribution of time (and treatment repetitions) to in-vivo TRAIL resistance, independently of the resistance that is related to differences in the cells context between in-vitro and in-vivo se ings. The second virtual experimental se ing we have considered is the treatment of cells grown in-vitro as -dimensional spheroids (Figure . ). This se ing is complementary to the first as it allows to investigate the contextual dimension (i.e., differences in the micro-environment felt by cells) of the gap between standard in-vitro and in-vivo experiments.

Testable predictions about the long-term response of HeLa cells to repeated TRAIL treatments

A first important model prediction that was obtained by simulating repeated TRAIL treatments on HeLa cells is that killing efficiency is strongly reduced between the first and later treatments, this effect being present at all doses that enable killing (Figure . ). This prediction has implications for dosing strategies, as here the dose needed to balance cell prolif- it is likely that the strong increase of long-term killing efficiency (with respect to the initial killing efficiency) that we predict for HeLa cells also holds for other cancer cell lines.

Despite the resistance acquisition between initial and consecutive treatments, simulations also predicted that the best strategy to eliminate or control the size of a cell population using TRAIL is to maintain a constant stimulation. The obtained killing rate would be lower than the maximum instantaneous killing rate obtained when applying a single treatment on a 'naïve' population, yet it would be higher than the average killing rate obtained on the long-term when repeating 'pulse' treatments. This prediction can also be tested experimentally.

Testing this could be relevant for the therapeutical use of DRA's, since the frequency of administration and the stability of the DRA (that spans more than an order of magnitude between rhTRAIL and DR /DR antibodies) dictate the shape of the temporal drug profile sensed by the tumor. This would require a precise control on the TRAIL concentration in the cell micro-environment. This can be achieved with microfluidic devices, an approach that was successfully applied to investigate NF-κB signaling (Kellogg & Tay ). However, such technology has not yet been applied to investigate TRAIL resistance (but one study used such a device to study the efficacy of combinatorial treatment of TRAIL and Doxorubicin (Kim et al. )).

Potential impact of non-canonical signaling

Although we did not check whether this predicted property (i.e. that constant TRAIL stimulation the best strategy to limit cell proliferation on the long-term) is valid across a wide range of parameters, we believe it is likely the case and that this property is related to the model structure rather than specific parameter values. In fact, different but comparable model structures might also lead to the same behavior, given that it was also predicted by our -parameters phenomenological model (Figure . ). ) that is masked by cell death at high doses but dominant at low doses.

Recently, TRAIL resistance acquisition has been experimentally characterized in detail in another cell line (MCF A, (Flusberg et al. )). One important finding of this study is that a cell population can gain resistance to TRAIL in absence of cell killing: pre-treatment with TRAIL in presence of a caspase-inhibitor (hence preventing cell death even if caspaseis activated) led to reduced killing (as compared to naïve cells) when a normal TRAIL treatment was subsequently applied. This shows that at least for this cell line, resistance acquisition cannot be a ributed solely to selection effects, and that non-lethal TRAIL treatments can increase the resistance of individual cells. This is consistent with our simulation results (Figure . ), which also showed (for our reference model parameterization for HeLa cells) that the contribution of selection effects to resistance acquisition is not dominant. Their data also shows that this cell-level resistance increase occurs at the level of or upstream of caspase-, as downstream signaling was blocked. Again, in order to get deeper insights into the mechanisms at play, it will be interesting to investigate this in HeLa cells and quantitatively compare the results with model predictions.

Dissecting distinct types of contributions to TRAIL resistance

In agreement with experimental evidence (Flusberg et al. ), our simulation results show that even in absence of genetic mutations, the resistance to TRAIL of a cell population is a not a fixed, intrinsic trait but is dynamically shaped by the TRAIL stimulation profile it is submi ed to. Thus, to predict the outcome of repeated TRAIL treatments, knowing the killing efficiency of single treatments is not sufficient. In addition, the 'capacity' of the cell line to gain resistance when treatments are repeated should also be quantified. This capacity can be defined as the ratio between the 'long-term' killing efficiency (i.e. for a population that has been already treated several times) and the 'initial' killing efficiency (for a naïve popu-lation). By studying the relationship between the model parameters and the corresponding values of those two quantities, we found that different types of parameter changes affect them differently.

First, modifying the mean expression levels of pro-or anti-apoptotic proteins has an impact on both the initial and long-term killing efficiency, however, the ratio between the two quantities remain unchanged. Thus, the extent to which cell population can gain resistance is not controlled by the expression levels of those proteins.

Second, modifying the rates of targeted degradation for non-native species (such as activated caspase , etc…) had a very strong impact on both initial and long-term killing efficiency, with half-lives varying from to hours predicting a large spectrum of highly resistant to highly sensitive cell lines (Figure . ). ). More generally, it is becoming more and more appreciated that protein degradation often plays a key role in the quantitative functioning of signal transduction pathways (Loriaux & Hoffmann ).

It should be noted that in our model, most of the targeted degradation reactions are modeled with constant rate first order reactions. However in reality, many biochemical steps are probably involved, with the rate-limiting step being controlled by another protein (in fact, this is already the case for XIAP-mediated degradation of activated caspase-, as in the original model (Albeck, Burke, Aldridge, et al. ); also, as we assume that ).

While our model predicts a strong impact of non-native species degradation on TRAIL resistance, both initial and long-term killing are predicted to be affected equally, as their ra-tio remained constant. Thus, such degradation is not responsible for the gain of resistance arising when treatments are repeated. Searching for parameters that are implicated in resistance acquisition per se, we found that the 'dilution' rate, i.e. the rate at which stable proteins (all except Flip and Mcl-in our model) are 'degraded' because of cell growth, strongly impacts the ratio between initial and long-term killing efficiency without affecting much initial killing efficiency. Thus, our model predicts that resistance acquisition is highly dependent on the turnover of 'stable' proteins.

Those results can be interpreted as follows: upon the first TRAIL stimulation, targeted degradation of activated pro-apoptotic proteins, allows some cells to survive at the price of decreased levels of native pro-apoptotic proteins (and hence increased resistance). Recovery of those levels occurs slowly, at the rate of cell proliferation, because those proteins are mostly stable. This view can be contrasted by noting that the decrease in native proapoptotic proteins could be accompanied by a similar decrease of anti-apoptotic proteins that participated in the targeted degradation, such as for R/Flip, Bid/Mcl-and Bax/Bcl in our model. However, because some anti-apoptotic proteins exhibit a fast turnover, they rapidly recover their steady-state levels and hence there is a transient imbalance towards relatively more abundant anti-apoptotic proteins.

To formalize this heuristic explanation at a level that abstracts specific biochemical reactions, we have proposed a -parameters phenomenological model (Figure . ). In that diction that constant stimulation is the best strategy to limit long-term proliferation. This generic model provides a conceptual link between molecular mechanisms and the highlevel, population-scale dynamics of cell killing for any kind of cytotoxic drug and hence might be useful in the context of high-dimensional drug screening studies where the response of many cell lines to many drugs need to be interpreted. It would be also interesting to develop a variant of this model in which a pleiotropic effect of the drug is accounted for.

Relevance of non-genetic resistance to DRAs in-vivo

It is clear that for isogenic cancer cell lines cultivated in-vitro for timescales of several weeks or below, the contribution of mutational mechanisms to resistance acquisition is limited. Resistance acquisition rather arises from transcriptional and post-transcriptional mechanisms, such as in our model. Still, when longer time-scales are considered, it is possible that the contribution of mutational mechanisms to resistance acquisition starts to be dominant.

However, when month-long in-vitro selection of TRAIL resistant 'clones' is performed by cultivating cells in media with escalating TRAIL concentrations, the resistance of those clones is often found to be reversed following culture in TRAIL-free medium (for example, in a study using leukemia cells it was the case for one out of two of tested 'clones' (Cheng et al.

)). Thus, non-mutational resistance acquisition is also relevant for month-timescales.

In addition to the role of timescales, the difference in cellular context between in-vitro and in-vivo cancer cell populations might also influence to relative contribution of mutational versus non-mutational mechanisms to TRAIL resistance. Although in most in-vivo experimental studies of TRAIL resistance, the data and protocol is such that it is difficult to conclude about the presence or absence of resistance acquisition when treatment is repeated over time, we found strong indications for non-mutational resistance acquisition to DRA treatment in a recent study (Huet et al. ): in Figure of the corresponding paper, the first treatment allows a reduction of the tumor volume, while its not the case for the second treatment applied one week later (for both types of DRAs tested). Thus, nonmutational resistance acquisition to DRAs is probably relevant in-vivo. This means that in order to quantitatively understand the dynamics of TRAIL resistance for in-vivo tumors, it might be needed to combine intra-cellular models of extrinsic apoptosis such as ours (capable of predicting cell-intrinsic, non-mutational resistance acquisition) with agent-based models of tumor growth (in which the effect of cell-extrinsic factors can be modeled).

The role of cell-extrinsic factors in resistance acquisition

As a step towards quantitative, cell-based multi-scale models of TRAIL resistance for invivo tumors, we have simulated the response of multi-cellular spheroids to repeated TRAIL treatments. Our results suggest that limited TRAIL penetration within a tumor is probably not a major source of resistance (Figure . ), as even when assuming slow diffusion and fast degradation of TRAIL within the spheroid, no strong differences of long-term outcome of repeated treatments were predicted for treatment strategies that were efficient on 'monolayer' (in the sense that all cells face the same environment) populations.

Still, it is yet well accepted that the drug resistance of D vs D cultures of the same cells can be strongly different (Pampaloni et al.

; Yamada & Cukierman ). If those differences do not arise from limited drug penetration, it means that cell-level differences in drug resistance exist. Those differences will be caused by differences in gene expression profiles that would themselves be caused by differences in cell micro-environment.

In the specific case of TRAIL, the differential response of monolayers versus spheroids has almost not been investigated. However, a study of the response of mesothelioma cell lines to bortezomib (which by itself can induce apoptosis in some cancer cell lines) has been conducted (Barbone et al.

). They found that the same cells were much more resistant when grown as spheroids, and found important differences in the expression levels of many pro-and anti-apoptotic proteins, but could not make specific links between those changes and the change in resistance. Our modeling approach provides a valuable framework to quantitatively investigate this question, and it will be interesting to quantify experimentally the differential response of HeLa spheroids and monolayers to repeated TRAIL treatments.

Even if we manage to successfully map changes in gene expression between the two culture conditions with the corresponding response to TRAIL resistance, the cause of those gene expression changes will still be unknown. Understanding the underlying mechanisms would potentially reveal molecular targets to prevent pro-resistance gene expression changes to take place. Striking evidence for such changes was obtained for

Ewing's sarcoma cell lines (Merchant et al. ). In this study, tumor cells implanted and grown on mice (in absence of treatment) were then explanted and grown in-vitro. Their resistance was much higher than the parental cell line, showing that simply growing in an in-vivo context led to the development of TRAIL resistance in those cells, and that this resistance is not lost immediately after change of the environment (but it was lost after several passages in-vitro, showing it is does not involve mutational mechanisms). The underlying causes for this non-TRAIL mediated resistance acquisition are yet unknown.

Again, our modeling framework could be useful to investigate them. An early success of multi-cellular synthetic biology is the implementation of population control in liquid cultures of E. coli bacteria (You et al.

Chapter

). This system relies on a simple design: cells are implemented with the ability to sense the size of the population and to undergo cell death for too high population sizes. This was realized by using the AHL cellcell communication system from Vibrio fischeri (AHL is a small, diffusible molecule that can be synthized or sensed by cells) in combination with a gene coding for the toxic protein CcdB (i.e. causing cell death when present in sufficient amount in the cell), which production was placed under the positive control of AHL sensing. This system worked well in the sense that different target densities could be chosen by tuning AHL half-life and that those target densities were maintened for durations as long as hours.

Engineering spatial multi-cellular systems

In this example, cells of the population were cultured in liquid medium. Several others multi-cellular synthetic systems were demonstrated for populations growing in liquid cul-tures, for example, a prey-predator system (Balagaddé et al.

) and a system involving altruistic cell death (Tanouchi et al.

). However, in such cases the cell population do not exhibit any spatial organization, hence it does not really form a tissue. For many potential applications of multi-cellular synthetic biology (such as tissue engineering, synthetic organogenesis, engineered biomaterials…), cells of the population are expected to be spa- ), but the spatial configuration of cells is itself not controlled, and is mainly dictated by the endogenous growth program and the environment. In other terms, the size or the shape of the tissue are not subjected to self-organization (other than natural).

Engineered spatial systems with such self-organization capabilities still remain very limited. One notable exception is the obtention of ring-shaped pa erns of the cell density in bacteria by combining cell-cell communication with a genetic circuit controlling cell micromotility (Liu et al. ). A possible explanation for this limited success is the involvement of multiple spatial and temporal scales (cell growth, proliferation and death; diffusion; signal transduction; gene expression; etc.) that make it difficult to reason about the system.

In addition, the impact of the (unavoidable) cell-to-cell variability on system functioning is hard to predict and might be fundamentally different between non-spatial (e.g., populations grown in liquid culture) and spatial systems.

Cell-based multi-scale modeling is a priori well suited to address those challenges, as it allows the explicit representation of individual cells in space and renders it possible to account for both cell-intrinsic (such as stochastic gene expression) and cell-extrinsic (such as gradients of messenger molecules) sources of variability in cell state and fate. In this Chapter, I present the application of cell-based multi-scale modeling for guiding the development of an artificial yeast pa erning system.

Case study: implementing spatial pa erning in yeast

In principle, spatial pa erns in the configuration of cells in a population can be generated Therefore, they proposed a modular implementation of the pa erning system (Figure . ).

In addition to two modules for IP production and IP sensing, it contains two cell death modules connected to IP sensing meant to function either for too high IP sensing (highthreshold killing) or for too low IP sensing (low-threshold killing).

Context and objective of the collaboration

When we started collaborating with the Weiss lab on this project, a genetic implementation of this system was already constructed by Jing-Jing Sun (Sun ), building upon work from Ming-Tang Chen (Chen ). A description of the corresponding genetic circuits is shown in Figure . . Briefly, high-threshold killing is achieved via excision of the vital gene pkc (involved in cell wall integrity pathway) while low-threshold killing is achieved via the apoptotis-inducing protein Bax.

However, testing, tuning and debugging of this system is difficult. First, a brute-force exploration of the system behavior is not feasible given the high dimensionality of the space of possible experiments (three different knobs for tuning the IP production and the killing modules; also, the initial spatial configuration of cells is in principle critical in determining its evolution; and finally the dimension and physical properties of the solid medium in which cells are cultured a priori also impacts system functioning by modulating the dynamics of IP gradients). Moreover, as already mentioned, intuitive reasoning about the system behavior is also difficult because it involves multiple spatial and temporal scales.

Hence, our objective in this collaboration was to use bo om-up cell-based multi-scale modeling to help building a quantitative understanding of the system functioning. Such understanding is expected in turn to permit more rational and efficient model testing, tun-ing and debugging. Importantly, we aimed to exploit as much as possible the system modularity.

More precisely, we used liquid culture data on the IP-induced response of receiver cells . Existing experimental data characterizing system behavior

Here we describe the experimental data that we used either to estimate parameters of our models or to verify their prediction capabilities. First, in order to quantitatively characterize IP sensing module behavior for receiver cells grown in liquid culture, a high-throughput flow cytometry experiment was performed by our collaborators to measure IP-induced GFP expression dynamics with high dose and temporal resolution and at the single-cell level (Figure . ).

However, the final system is meant to function in solid cultures and in the presence of IP gradients. In addition, such gradients should exist at a sufficiently small scale, for example between close micro-colonies. To characterize the behavior of IP sensing for receiver cell micro-colonies in presence of IP gradients, an experimental setup making use of a cell printer (a technology that enables to deposit very few cells at defined locations on a surface with micrometric precision) was developped by our collaborators (Figure . ). As we will see, a differential response in space and time was obtained between micro-colonies, indicating that sufficient gradients of IP to allow system functioning can exist at this spatial scale.

Finally, the functioning of the high-threshold killing module was tested by printing patterning cells on a agar gel containing Dox (needed for its functioning) and IP (for activating it) (Figure . ). This data showed significant killing in the induced micro-colonies, albeit many hours after seeding cells. Interestingly, not all cells died, and the localization of dead or alive cells in the micro-colonies does not seem random: large clusters of dead cells are apparent. Because ultimately, the desired behavior of the system is the obtention of clear, dynamic, predictable spatial pa erns of dead or alive cells, it would be highly valuable to gain a quantitative understanding of the mechanisms driving the apparition of such death clusters.

IP induction levels

Induction durations

!"#$% !"#"$%"&'#"(()'
IP sensing module driving GFP expression *$+,("-#"(('&").(%"/'0"12.&3('/.)"' &")2.+)"'0.'45'

per IP dose X time point .

Modeling system behavior

In this section we describe the different computational models of the system that we developped. Their calibration and the assessment of their prediction capabilities will be presented afterwards.

A semi-mechanistic ODE model of IP sensing

IP sensing is central to the functioning of pa erning cells because both killing circuits rely on it. Hence, we focused first on obtaining a quantitative model of IP sensing that will be calibrated from the data characterizing receiver cells response to IP in liquid culture (Figure . ).

Because the precise mechanisms by which yeast cells expressing the (exogenous) AtCre receptor can activate the nuclear factor Skn in response to IP are yet not fully understood,

!"#$%&'&()*$++,)

all modules

Cell printer (10X)

(10X) T = 0h T = 14h T = 38h T = 0h T = 14h T = 38h
Cu -/ Dox -/ Gal -/ IP -

IP production if Cu HT killing if Dox LT killing if Gal

Cu -/ Gal - Dox + / IP + Performed by Brian Teague at MIT. Patterning cell micro-colonies were seeded using a cell printer on a agar gel. All gels also contained Phloxine B (a dye that accumulates in dead cells) and imaging was performed such that dead cells appear black. In a control experiment, the gel did not contain Dox, Cu, Gal or IP and the cells eventually formed large micro-colonies (right). To test the functioning of the high-treshold killing module, a gel containing Dox and IP was used. Strong killing occured between 14h and 38h and resulted in clusters of dead and alive cells.

we follow the approach adopted in (Palani & Sarkar

) by considering a semi-mechanistic ODE model of IP signaling, in which a single intermediate between receptor binding and transcription activation is represented (Figure . , top). The corresponding reactions are given in Table . . In addition to this structural choice, we also decided to normalize the expression levels of R, TF and (basal) GFP in order to minimize identifiability issues.

Semi-mechanistic ODE model of IP sensing (to compare to population data)

Model extension accounting for stochastic gene expression

(to compare to single-cell data)

!"# $# !"%$# &'# &' (# " &$))$*# +,'" # ,'" # …
σ uGF P → GF P GFP maturation kf R → ∅ R degradation γ dil + γ R → R R synthesis γ dil + γ R (*) T F → ∅ TF degradation γ dil + γ T F → T F TF synthesis γ dil + γ T F (*) uGF P → ∅ uGFP degradation γ dil + γ uGF P → uGF P uGFP synthesis (γ dil + γ uGF P + kf ) γ dil

Model extension accounting for stochastic gene expression

This model is deterministic and is meant to be compared to population averaged data. However, under this form it cannot capture cell-to-cell variability in the response to IP, while such variability might be important for the functioning of the full pa erning system. As we have shown in Chapter , combining this model with randomly distributed initial conditions to mimic cell-to-cell differences in protein level is not appropriate for timescales of the order of the cell cycle or longer. Therefore, we followed the same approach as for TRAIL-induced apoptosis by extending the model with stochastic gene expression ( see Figure . , bo om).

Importantly, as opposed to our model of TRAIL-induced apoptosis, here one the protein species regulates the transcription of another gene. We propose to model this transcriptional regulation by allowing the transcription rate of GFP (i.e., the stochastic rate at which mRNA are produced when the gene is active) to depend on the time-varying concentration of the transcription factor TF*. We assume that this dependence can be modeled with a Hill function. Note that under those assumptions, cells in which the promoter driving GFP expression is inactive would not directly produce mRNA when the transcription factor becomes activated, and that the transcription factor do not modulate the rates of promoter activity switches.

Cell-based multi-scale model of receiver cell micro-colonies

To model the behavior of receiver cell micro-colonies in presence of IP gradients, we proceeded essentially as for modeling TRAIL treatments on multi-cellular spheroids (Chapter ). More precisely, we integrated an instance of the single-cell IP sensing model into each cell of spatial, agent-based multi-cellular simulations in which cell growth, division and motion as well as diffusion of IP are represented (Figure . ). Note that here the cell-based paradigm is advantageous in the sense that the setup of the solid culture experiment can be precisely reproduced (seeding of individual cells at determined locations, geometry of the gel) and that a single-cell model of IP sensing calibrated from liquid culture experiments can be integrated in a straightforward manner.

Towards modeling pa erning cells behavior

Here we describe an extension of the above models in order to model the behavior of patterning cells instead of receiver cells. We decided to focus first on the high-threshold killing module, hence mimicking conditions in which IP production and low-threshold killing are not activated, as in the experiments described in Figure . . Then, a corresponding spatial cell-based multi-scale model is constructed as for modeling receiver cell micro-colonies. Importantly, in order to closely mimic the experimental reality, dead cells are not removed from the simulation and can still be passively displaced.

Single-cell model of IP sensing 3D diffusion in the gel

Agent-based model of receiver cell micro-colonies

One instance of the single-cell model per cell . Data analysis and parameter estimation

Analysis of liquid culture data

To allow comparison of our IP sensing models with the experimental data characterizing the response of receiver cells in liquid culture, we extracted from the raw FACS data singlecell distributions of measured fluorescence corresponding to each IP induction dose and induction duration (Figure . ). Those distributions revealed an important cell-to-cell variability in IP response. For example, with a dose of 0.18 µM and hours after induction, about half of the cells exhibit strong fluorescence while the other half have fluorescence levels close to uninduced levels. The shape of this bimodal distribution remains largely unchanged during the following hour.

From those distributions, the population level induction dynamics can be quantified by simple averaging (Figure . , bo om). Our strategy to calibrate our single-cell model of IP sensing is as follows: first, we fit the ODE model to the population level induction dynamics data to obtain a first estimate of the IP sensing 'signaling' parameters. Then, we ask whether the extended, single-cell model can predict the single-cell induction dynamics data under reasonable assumptions for the stochastic gene expression parameters.

Estimation of signaling parameters from population-level induction dynamics

Fi ing the IP sensing model on the population level induction curves identified a parameterization allowing a good agreement with the data (Figure . ). The corresponding parameters are given in Table . . During the fi ing procedure, two parameters of the model were not varied because they were known: the dilution rate was computed from the observed doubling time and the IP half-life was estimated to be hours in a previous study (Chen

& Weiss

). Of note, at very high induction (IP = 5 µM ), the experimentally observed population level response is slightly lower than the one observed for IP = 3.3 µM . Because this non-monotonic behavior very likely corresponds to an effect that cannot be captured by the structure of our model, it was excluded from the fi ing procedure.

Interestingly, the fi ed maturation rate for GFP is realistic, with τ m = ln(2)/kf = 12.5 minutes: for comparison, a variant that folds faster (sfGFP) was estimated to mature in .

minutes (Khmelinskii et al. ). This consistency suggests that overfi ing was avoided.

Of note, the IP receptor is estimated to be short-lived (half-life of . minutes). We could not find previous estimates of this parameter (remember that AtCre is not endogenously expressed by yeast cells). 4.7) was used to reproduce the population-level induction dynamics curves (data normalized to uninduced time-course, the 5 µM IP dose was excluded as explained in text). The Downhill simplex method was used (implementation from Numerical Recipes). Dilution rate and IP degradation rate were fixed based on previous estimations. Fitted parameter values are given in Table 4.2. In the plot, data points are squares and continuous lines correspond to simulated induction curves with the fitted parameters. Table 4.2: List of parameter values for the semi-mechanistic model of receiver cells behavior calibrated to liquid culture data. (*) Specific degradation rates that were estimated to be lower than 10 -5 hr -1 were set to zero for simplicity, as this did not affect the model behavior (dilution is then largely predominant).

Parameter symbol

Value Comment ;348/&4$4#%0'0&

Only the outer rim of cells is kept to avoid biases from piling up at colony centers presence of IP gradients. 4.2). The simulated spatiotemporal profiles for the three colonies shown in Figure 4.12 are compared to the corresponding data. Good agreement is obtained for D IP = 10 -4 mm 2 .s -1 . In simulations, for simplicity we used the IP concentration at the colony center as being representative of the whole colony response.

. Reproducing the spatio-temporal response of receiver cell

micro-colonies

In liquid culture, cell-level response to IP induction was found to be highly heterogenous, and our single-cell model of IP sensing was able to reproduce this variability (Figure . ). However, in FACS data no information about the population structure is available, while for micro-colonies grown on solid medium, close cells are more likely to be closely related.

Because new cells inherits molecular traits from their mother cell, such as protein concentrations, a certain spatial correlation of the fluorescence response within a micro-colony is expected (independently of IP gradients). This might have important implications for the functioning of the pa erning system. Therefore, we asked whether our cell-based multi-scale model of receiver micro-colonies could predict cell-to-cell variability in fluorescence and its spatial correlation. Hence, we reproduced in-silico the experiment and compared the resulting images to data (Figures . , . and . for colonies A, B and C respectively). The overall agreement appears good, and clustering of 'responders' or 'non-responders' seems to be qualitatively reproduced by the model. Of note, the growth of receiver micro-colonies (with respect to cell-level cell cycle durations and colony-level -dimensional growth) has not yet been precisely calibrated, but we don't expect strong qualitative differences.

Finally, we confirmed the visual impression of spatial correlation in the simulations by computing the Moran's I spatial correlation metric within micro-colonies (Figure . , values of I can range between -(negative spatial autocorrelation) and (positive spatial autocorrelation), indicates a random, non-correlated spatial pa ern). Unfortunately, we could not apply the same quantification of spatial correlation on the data because single-cell segmentation at this magnification was not reliable.

. Predicting death clusters in pa erning cell micro-colonies This carefully constructed cell-based multi-scale model of the behavior of receiver cells from the single-cell level up to spatially organized micro-colonies in presence of IP gradients establishes solid grounds on which to build a model of the full pa erning system. To illustrate this, we investigated under which conditions the clearly apparent clusters of dead/alive cells that were observed experimentally (Figure . ) can be reproduced by the extended model accounting for high-threshold killing (as described at the end of section . ). In each case, a negative control is obtained by shuffling the GFP levels among cells in order to evaluate the statistical significance of the results. For times below 4 hours, the cell number per colony is too low to conclude to a spatial correlation.

parameters, even when the overall killing is similar. Larger death clusters are obtained for higher tolerance to low Pkc levels, corresponding to longer delays between pkc deletion and actual cell death. Thus, while we already showed that cellular variability and memory in IP sensing can promote a spatial correlation within micro-colonies of the protein expression in response to IP, we identified for the high-threshold killing module an additional effect, e.g. a long delay between cell commitment to death and actual death, that promotes the apparition of clusters of dead/alive cells.

In principle, quantitative comparison of such simulations with corresponding data would allow us to infer parameter values for the high-threshold killing module. Note however that pa erning cells grow much slower than receiver cells. This might have quantitative consequences not only on the overall growth of micro-colonies, but also on the functioning of IP sensing and high-threshold killing because the corresponding slower dilution can impact both cellular variability and memory and the overall dynamics. Work is still ongoing to quantify and understand be er those growth differences such that they can be incorporated in the model.

. Discussion

In this Chapter, we illustrated the potential of our cell-based multi-scale modeling framework to build predictive models of engineered multi-cellular spatial systems via a rational, bo om-up approach. We have focused on an artificial sensing device implemented in yeast This approach was successful: a very small number of parameterization tests was then sufficient to obtain a good semi-quantitative agreement between the single-cell model and . ). To our knowledge, the dynamic, quantitative functioning of synthetic sensing circuits have rarely been characterized with such accuracy at the single-cell level. Those results adds to the potential of the modeling approach proposed in Chapter (which consists in extending ODE models of the protein-protein interactions involved in a given signaling pathway with parsimoniously parameterized stochastic gene expression models for all the proteins) to be broadly applicable.

In addition, the approach is here extended, as one of the protein modeled acts as a transcription factor and hence impacts the stochastic expression of another protein, while no transcription regulation was considered in Chapter . In the case of the IP sensing circuit in receiver cells, we have found that assuming a TF-dependent transcription rate (k on ) was sufficient to explain the data. We note that in mammalian cells, transcription factors have also been reported to be able to modulate transcriptional bursting (k on and k off ) (Senecal et al.

; Molina et al. ). It is likely that the effect of the transcription factor concentration on the values of k on , k off and k sm does not follow an universal rule but depends on the specific mode of action of this transcription factor and also on the location of the gene on the chromosome. When needed, specific dependencies other than the one used here could be introduced and tested.

Early on, synthetic biologists felt that having a reliable workflow to build quantitative models of the behavior of gene circuits was needed for allowing the successful design of circuits achieving complex functions (Purnick & Weiss ). One major hurdle towards this goal is the ubiquitous cell-to-cell variability in circuit behavior (Miller et al. ; Toni & Tidor ), which limits the potential of purely ODE-based strategies to characterize gene circuits. Our modeling approach could enable the systematic characterization of the quantitative, dynamic behavior of synthetic circuits with single-cell resolution. Other approaches have also been proposed to infer the rates of stochastic models of induced gene expression, for example based on moment closure techniques (Zechner et al.

). Although our approach is presumably less powerful, it constitutes a technically simple, approximation-free alternative to such approaches, which here was proven successful on a 'real-life' example.

Predicting circuit behavior in spatially organized cell populations from its characterization in liquid culture

The proposed approach enabled to quantitatively reproduce the response to IP of receiver cells grown in liquid culture. However, in order to be useful in guiding the design of the pa erning system, this model should be able to predict the behavior of receiver cells grown as spatially organized populations. Such predictive power is not guaranteed, as cells grown in solid medium might exhibit physiological differences compared to cells grown in liquid cultures, and such differences might in turn affect the functioning of the IP sensing circuit via context-dependency, a central issue for synthetic biology (Del Vecchio

). Whether

or not we can successfully transpose quantitative knowledge about a system obtained from liquid culture experiments to predict its functioning in solid culture is particularly important, because the former are relatively easy and rapid to perform as compared to the la er, hence greatly increasing the number of circuits that can be characterized in a given amount of time.

Here, we have found that a model of IP sensing calibrated from liquid culture data could predict the response of receiver cells micro-colonies grown in solid medium submi ed to IP gradients: indeed, simply reproducing the diffusion of IP in agar without changing IP sensing parameters led to an agreement with the data (Figure . and . ). Note that in this case, obtaining such agreement was not guaranteed, even if the physiology of the cells does not change between liquid and solid culture conditions. This is due to the fact that in the liquid culture data used to calibrate the model, cells were submi ed to IP levels that were monotonically (and exponentially) decreasing with time. On the contrary, in the solid culture experiment cells experienced IP induction profiles of different shapes, increasing first as IP diffuses to a micro-colony, and then decreasing as it diffuses away and get degraded.

In fact, an earlier version of the ODE model of IP sensing was able to explain the liquid culture data as precisely as the current one, but resulted in a signficantly bad agreement with the solid culture data (not shown). In general, performing de-induction dynamics in liquid might be useful to eliminate such incorrect model parameterizations.

While here our results show that cell physiology changes between liquid and solid culture conditions do not strongly impact the behavior of the IP sensing circuit, it is unlikely to be always the case. For example, changes in growth rate could significantly impact the quantitative behavior of a given circuit, via changes in dilution rates or global feedback on transcription (Klumpp et al. ). In such conditions, knowledge gained from liquid culture data could still be highly valuable in predicting the functioning of the circuit in solid culture, given that the model is sufficiently mechanistic and was cautiously calibrated to avoid overfi ing. Indeed, then in principle only a subset of the model parameters, associated with relevant differences in cell physiology, would have to be adapted in order to capture the circuit behavior in solid culture (for example the dilution rate).

Another central question for synthetic biology is the modularity of the synthetic genetic circuits (Del Vecchio ): does a circuit constructed from the combination of two smaller circuits behave as one would predict by composing two models characterizing their behavior in isolation ? How precisely the concept of modularity should be applied to enable synthetic biologists to build increasingly complex systems is still an open question (Neal et al.

).

Here, the pa erning system (Figure . ) rely on connecting the IP sensing module with two killing circuits, and the IP-responsive promoters (P TRSSRE and P SSRE ) are the genetic elements realizing that connection. In principle, if the IP sensing model calibrated from data on receiver cells faithfully represents the reality of the involved biochemical reactions, then a model of the high-threshold killing in which the signaling parameters (from IP binding to activation of the transcription factor), and the stochastic gene expression parameters for the receptor and the transcription factor remain the same as in the IP sensing model, should be successful in predicting the high-threshold killing circuit behavior; provided that additional parameters needed to describe TF-dependent expression of Cre, pkc gene deletion and the lethality associated to low Pkc levels are calibrated.

We have here tested this idea by constructing such extension of the IP sensing model and we have found that death clusters (within micro-colonies submi ed to IP) that are very similar to experimental observations can be obtained (Figures . and . ). We also obtained mechanistic insights about the system behavior, by pointing out the differential effect of Cre efficiency and Pkc -associated lethality on the shape of death clusters. Such insights are valuable for the tuning and debugging of the killing circuits in order to get a desired behavior.

However, we did not yet fully characterize the high-threshold killing behavior from data because the growth of pa erning cells is markedly slower than receiver cells. The causes for such difference and the expected impact on circuit functioning (which might be significant and non-trivial, although this problem has been studied mainly in bacteria (Klumpp &

Hwa

)) are still being investigated. Also, our models currently account for cell-to-cell heterogeneity in growth in a simple manner that might overlook potential coupling with the system circuits behavior. Such potential coupling is also under investigation.

Towards modeling the full pa erning system

While we focused first on characterizing the IP sensing module (because of its central role in the design), and then started to characterize the high-threshold killing module (because in principle a single IP-dependent killing circuit could generate spatio-temporal pa erns of dead/alive cells when combined with IP production and diffusion), the IP production module and the low-threshold killing module remain to be characterized in order to obtain a complete, bo om-up, model-based characterization of the pa erning system (Figure . ).

Regarding the IP production module, estimation of the IP production rate as a function of Cu induction ( Figure . ) has in fact already been realized in liquid cultures, via an experimental setup that utilize receiver cells as quantitative sensors of IP concentrations. However, the quality of this characterization remain to be tested in solid culture experiments. For example, pa erning cell micro-colonies, induced only for IP production, could be printed on a gel together with receiver cell micro-colonies positioned at different distances, the la er responsing to the spatio-temporal IP gradients created by the former, similarly to the exper-iment shown in Figure . .

Once such a complete model of the pa erning system is constructed, it can be used to predict the outcome of well-defined experiments (gel geometry, initial positions of the microcolonies, induction levels (Cu, Dox and Gal) for the different modules) involving pa erning cells. Iterative comparison between predicted and experimental results would allow refinment of the model until a deemed satisfying predictive power is reached. Then, extensive in-silico exploration of the system functioning can be performed to identify promising initial configurations and levels of induction in terms of the spatio-temporal pa erns of dead/alive cells that are predicted.

At this stage, it is still not clear what kind of pa erns can robustly be obtained with the current implementation of the system. It is possible that fundamental limitations apply, originating for example from the IP diffusion coefficient and stability, the delay between commitment and cell death, the overlap of those timescales with cellular growth and division, etc. Such limitations do not necessarily impact different implementations of the system equally (in other words, some implementations might function more robustly than others).

While here we used a bo om-up approach to gain a detailed, quantitative understanding of the behavior of a particular implementation of the system, it might be useful to use a top-down approach to derive simpler, more generic models of the system from its highlevel design (Figure . ). Such models can be used to broadly (but coarsely) explore the 'design space' of the system, and to provide hints on how robust system functioning could be implemented.

Chapter

Simulation of cell-based multi-scale models

The task of constructing and simulating the cell-based multi-scale models of TRAIL-induced apoptosis (Chapters and ) and of the yeast pa erning system (Chapter ) was challenging. Indeed, those models are high-dimensional and involve a large number of equations.

Moreover, their dynamics is described by a mixture of different mathematical formalisms: ordinary differential equations (for example for intra-cellular signaling reactions or the motion of cells in D), continuous time Markov chains (for example for stochastic gene expression), agent-based/population modeling and partial differential equations (for example for the diffusion of extra-cellular molecules). This renders the efficient description and simulation of those models a difficult and challenging problem. To our knowledge, no existing tool directly enable the description of simulation of such complex models.

The goal of this chapter is to describe methodological difficulties that appeared critical during this thesis, and to discuss potential solutions. We start by describing challenges relative to the choice of adequate simulation algorithms, considering the problem of intracellular models first and of the whole multi-cellular simulation afterwards. Finally, we also ).

However, as we argued in Chapter , the problem with this approach is that cell-to-cell variability cannot emerge from the cell itself and can arise only from external cues, in contradiction with experimental observations. A key concept in our modeling framework is to allow for cell-instrinsic sources of cell-to-cell variability by simulating stochastic gene expression for all the proteins that take place in the cellular decisions.

This choice poses a computational challenge regarding the simulation of cell state. The Gillespie algorithm requires to interpret all reactions as being stochastic and is highly inefficient when some species are abundant. Because of the marked separation in copy numbers between genes and mRNAs on the one side and proteins on the other, it appeared reasonable to model all reactions as deterministic except gene activity switches, transcription and mRNA degradation (see Chapter , Figure . ). This type of approximation is referred to as 'partitioning' and has been reviewed elsewhere (Haseltine & Rawlings ).

An alternative would have been to use an approximate stochastic simulation algorithm, such as tau-leaping (Cao et al.

). One advantage of our choice is that the approximation is done at the level of the mathematical model rather than its simulation, which might facilitate reproducibility of the results, and allows to use existing analytical results (Paszek ) to judge the validity of this approximation.

method ((Hairer et al. ), which is an eigth-order Runge-Ku a method. The other is an the Semi-Implicit Extrapolation (SIE) method ((Deuflhard )). It has been developped

for stiff problems, and requires that an expression for the jacobian is provided (although it can be an approximation).

Computational efficiency

Using our model of TRAIL-induced apoptosis (see Chapter ) as a real-life example, we benchmarked the computational efficiency of the mRNA storage simulation algorithm for the two ODE solvers described above (Figure . ). In each case, different storage durations (i.e. for how long the genes and mRNA are simulated alone and the corresponding mRNA trajectories stored) were tested.

The results strongly depended on the ODE solver method. For the Runge-Ku a based method (Dopr ), the storage duration has almost no impact on the computational efficiency. Moreover, this efficiency is very similar to the one obtained with the naïve algorithm (in which the ODE solver is called between each stochastic event).

However, the results are very different for the SIE method. First, the computational efficiency strongly depends on the chosen storage duration, up to an optimal storage duration of about minutes. For those storage durations, the gain compared to the naïve algorithm is large. Notably, the total gain of using a sophisticated solver and storing mRNA trajectories is more than -fold and is much larger than the gain of ODE solver change when the naïve algorithm is used.

In summary, the choice of the ODE solver and the simulation algorithm is critical for the computational efficiency. Sophisticated methods such as SIE, which are already faster than more common methods with the naïve simulation algorithm, might strongly benefit from the mRNA trajectory storage strategy, probably because they usually have a strong 'fixed' cost that is paid each time the solver is called. This speed-up has been instrumental for our investigation of TRAIL-induced apoptosis presented in Chapter and , as it allows to explore more regions of the parameter space, and to simulate more cells and longer.

A small software tool allowing the generation of c++ code for the combined simulation of signal transduction and stochastic gene expression from a programmatic description of the model reactions has been deposited on GitHub: http://github.com/fbertaux/FluctuProtST (it might be slower than the code used for hEARM in this thesis as it does not yet allow the use of the SIE solver, which requires the generation of the Jacobian function).

storage duration for mRNA trajectories (min)

Parameterization of stochastic protein turnover models from protein variability and mixing time

As we have seen in Chapter , the transcriptional bursting parameters (k on and k off ) are usually unknown, not measurable directly and difficult to infer. However, knowing their exact values might not be needed if we are only interested in reproducing the protein level fluctuations, as is the case in most applications.

For this reason, we have proposed to use the protein level coefficient of variation and half-autocorrelation time as two 'statistics' that quantitatively characterize the fluctuation phenotypes of a given protein. This would mean that when the impact of the fluctuations of a given protein on the overall response has to be studied, one could explore the twodimensional space CV × τ instead of the larger -dimensional space of stochastic protein turnover models (we excluded the mean expression level as it is usually available by standard population-level measurements and does not impact CV and τ ).

To enable this kind of parsimonious exploration, a procedure allowing to find a set of parameters that are compatible with a given CV and τ is needed. Relying on the analytical expressions for CV (given in Chapter ) and τ (given in Appendix ), we have developped such a procedure. Note that constraining CV and τ is a priori not sufficient to uniquely determine a single parameter set, as three degrees of freedom remain. For example, if one fixes the value of the rates γ p , γ m and γ g = k on + k off (an 'aggregate' switching rate for gene activity), the two parameters EG = k on k on +k off (mean time fraction the gene is on) and EM = EG k sm γ m (mean mRNA level) remain free and might be chosen to match the two constraints on CV and τ .

By examining the analytical expressions for CV (given in Chapter ) and τ (given in Appendix ), we noticed the following:

• The expression giving CV 2 can be decomposed in a 'gene' and 'mRNA' part:

CV 2 (γ g , γ m , γ p , EG, EM ) = CV 2 gene (γ g , γ m , γ p , EG, EM ) + CV 2 mRN A (γ m , γ p , EM )
• The first part CV 2 gene of this expression varies from from 0 + to +∞ when EG is varied (between 0 + and 1 -) and other parameters are fixed • The second part CV 2 mRN A depends only on γ m , γ p and EM , and varies from 0 + to +∞ when EM is varied and γ m , γ p fixed • The protein autocorrelation function can be expressed solely as a function of the three rates γ p , γ m and γ g and the 'gene contribution ratio to CV 2 ', that we call α = CV 2 gene CV 2

• Thus, when γ p , γ m and γ g are fixed, all the possible values for the mixing time τ can be enumerated by varying α between 0 + and 1 -. τ seems to vary monotically with α These findings led us to a procedure (described in Figure . ) that always find the couple (EG, EM ) (when it exists) compatible with given CV and τ constraints, when γ p , γ m and γ g are fixed. We implemented this procedure in python, and the code is available at http: //github.com/fbertaux/sgeLytics.

To illustrate the type of analysis that can be done with this procedure, we investigated how different parameterizations that share the same CV and τ occupy the parameter space (Figure . ). The slice γ mγ p shows that both mRNA and protein half-lives cannot be large together to meet the imposed constraints, probably because of the τ value, which increases with each half-life. The slice EG -EM show a region (large EG and large EM) for which no solution to the CV and τ constraints exist.

!"#$#%# & ' $#& ( $#& )# *+#,-#./0-## %#1#% (23 4& ' $#& ( $#& ) 5 # /36# #%#7#% (/8 4& ' $#& ( $#& ) 5#9# # :#;23)<-#;+<=>2+3# 4?@$?A5#-82;>;# # B+#;+<=>2+3#>+#>.-#%# C+3;>D/23>#)20-3#>.+;-# D/>-; # !"#$%&'()#*"+(#)%% ,-./&#/"#+'()0%&'()#*"+(#)% E# F+=36-6#;C/</D#+'>2(2G/>2+3#+3# E#>+#(-->#>.-#%#C+3;>D/23>## BH# I?J# 12% 13% !+('=>-#?A#;=C.#>./># !" K (LB:# M#4N#O#E5#P#!" K# !+('=>-#?@#;=C.#>./># !" K )-3-M#E#P#!" K#

Simulating regulated gene expression

While this is not the case for the canonical TRAIL-induced apoptosis pathway, signal transduction pathways often involve transcriptional regulation, via signaling reactions that change the activity and/or localization of transcription factors. As a result, the transcription of target genes for those transcription factors is altered (often positively). The products of the target genes can either participate in the same pathway, hence implementing a feedback in the pathway response; or they can be involved in other pathways. For example, in N F -κB signaling, the cytosolic/nuclear localization of a protein complex that can bind DNA and activate the transcription of target genes is controlled by other signaling proteins that respond to external stimuli such as the cytokine T N F -α. Among the target genes,

IκBα and A20 are involved in the signaling pathway responding to T N F -α (Tay et al.

).

In such context, if one is interested in modeling the expression of a gene that is targeted by a transcription factor, the rates of the two-state transcriptional bursting model cannot all be assumed to be constant: instead, at least one and maybe all of the rates k on , k off and k sm depends on the level of the transcription factor, depending on the molecular mechanisms by which the transcription factor influences transcription. For example, in our model of the This has implications for the computational simulation of the model. Because the propensities of some stochastic reactions now evolve continuously according to the differential dynamics of the continous variables (the protein levels), the Gillespie algorithm, which rely on a constant vector of propensities in-between two reactions, cannot be used anymore. An 'exact' simulation algorithm exists: it relies on integrating the propensities of stochastic reactions together with the differential dynamics of the continuous dynamics (Salis & Kaznessis ; Alfonsi et al. ). But this algorithm is difficult to implement and requires to solve an event detection problem to stop integrating propensities when one of the stochastic reactions has fired, which might have a high computational cost.

Alternatively, a pragmatic, approximate solution is to adapt the naïve algorithm described previously (or the storage algorithm) with a maximum time step and update the rates of the stochastic reactions according to the values of the continuous variables at each time step. In this case, the algorithm approximates those rates as stepwise functions, with the size of the steps depending on ) the chosen maximum time step and ) the rate of changes of the continuous variables involved in the stochastic rate functions. Thus, the accuracy of the approximation can be tested by varying the maximum time step and evaluate whether the simulation results change significantly or not.

We used this approximate algorithm for the IP sensing simulations (See Chapter ) and found that for time steps below minutes, the results did not depend on the time step any-more (not shown). Thus, this appears to be a viable solution. However, it might not be computationally optimal, and further investigation is required to compare the computational efficiency with the 'exact' algorithm. Still, the computational efficiency of the approximate algorithm was sufficient for our needs, as we were able to simulate the full multi-scale model of spatially-organized receiver cell micro-colonies (Chapter , for example Figure . ) in reasonable computational times on a standard computer.

. Simulating spatially organized cell populations Because many of the computational challenges that arise when simulating cell populations with an agent-based framework have been described previously, notably by previous members of our group (Hoehme ; Jagiella ; Weens ), we will here simply highlight points that appeared critical for the present work.

Simulating cell motion

The spatial arrangement of multi-cellular assemblies such as spheroids or micro-colonies is constantly changing because of cell growth, division, death, migration through physical interactions between cells in contact as well as between cells and their environment. To simulate this, we have used in Chapters and an existing la ice-free approach in which cells are assumed to behave as visco-elastic spheres (Drasdo & Höhme ; Drasdo et al.

).

In this approach, an equation of motion relates at each instant velocities of individual cells with cell-cell contact forces ( . ), which values depend on the geometry of those contacts.

Thus, the positions of all cells evolve according to a set of 3N coupled differential equations of the form

dx i dt = F x i (x 1 , y 1 , z 1 , ..., x N , y N , z N ), …, dz i dt = F z i (x 1 , y 1 , z 1 , ..., x N , y N , z N )
where N is the number of cells. However, solving this system with a standard ODE solver is challenging because ) N is not constant and can be very large compared to traditional ODE systems and ) the right hand side is generally costly to evaluate, as it requires that all cell-cell contacts are known.

The first point was originally adressed by using a simple first-order, explicit scheme to solve cell motion, with time-step control based for example on a maximum distance that cells are allowed to travel in a single time-step (Weens , used here in Chapter ). (Of note, in a large population of cells where division takes place, this time-step choice strategy is not ideal as there will always be a division event imposing a 'worst-case' time-step.

More sophisticated strategies that prevent local constraints caused for example by division to have a global impact on the computational efficiency of the simulation would be beneficial.) Alternatively, we also adapted the implementation of existing ODE solvers to solve the cell motion equation in synchrony with cell growth as well as intra-cellular dynamics (more details are given in a following section).

The second point is critical for computational efficiency, as using a naïve contact detection algorithm (in O(N 2 ) where N is the cell number) will prevent the simulation of large cell populations (~ -K cells on standard machines) in reasonable computational times. A efficient algorithm is to partition space into cuboïdal 'contact' voxels which size is chosen such that cells in contact have necessarily their centers in neighbor voxels (here neighbors in the sense of sharing one or more vertices). Then, contact detection is linear in N (provided their is a limit to the number of cells that can belong to the same voxel): in each step, each cell is mapped to its containing voxel, and then contacts are tested only for cells in either the same voxel or neighbor voxels. More details about such algorithm can be found elsewhere (Drasdo et al.

).

An additional difficulty arises when cell-cell friction is accounted for. In that case, terms involving the velocities of cells in contact appears in the equation of motion ( . ). As a result, the velocity of each cell cannot be extracted directly, instead a non-diagonal linear system A V = B has to be solved. Given the size of the system, the associated computational cost can be very high and adds up to the contact detection cost in the total cost of the derivatives computation for the ODE system of cell positions. An approximate method such as the conjugated gradient is likely a good choice to optimize this cost. More details on this can be found elsewhere (Weens ).

Another difficulty also arises when a stochastic term is used to represent micro-motility.

If micro-motility is modeled as brownian motion, the resulting equations are stochastic differential equations. Depending on the numerical scheme that is used, this might cause problems, notably when a time step control procedure is used to limit the maximum distance traveled by a cell. In this thesis, we have decided to model micro-motility (see Chapter ) as a persistent random walk instead. Briefly, the migration force (direction and magnitude) is stochastically changing at a given rate but is constant between two consecutive changes (the 'persistence' time between changes is exponentially distributed and in average is equal to the inverse of the change rate). Oppositely, for the brownian motion model, the migration force is changing within all potential time windows. When a change occurs, a new migration direction and magnitude is sampled into distributions. This facilitates the numerical resolution of the motion equations.

Simulating the diffusion of extra-cellular molecules

To simulate the diffusion of TRAIL within spheroids (Chapter , see Figure . ) or of IP in the solid medium on which yeast micro-colonies were grown (Chapter , see Figure . ), we have taken a simple approach: we used a finite differences, Euler (first-order, explicit) numerical scheme on regular cuboïdal grid. This scheme is robust (stability is ensured when the CF L criterion is met and accuracy is usually satisfying under those conditions) and easy to implement.

However, this method can be computationally burdensome when a fine spatial resolution (relative to the domain size) is required. In our case, we usually used a spatial resolution (i.e., the length of voxels composing the discretization grid) of about a cell diameter.

We believed it was not meaningful to go below that value as the surface of our cell agents is assumed to be homogenous and hence was not spatially resolved.

If needed for more complex models (for example in which cell surface receptors are not distributed homogenously on the surface, etc…), more sophisticated methods such as finite element methods with adaptative meshes can be much more powerful, but they are more difficult to implement and existing libraries are not straightforward to handle. In our group, a cell-based multi-scale modeling software currently in development (TiSim) will have such capabilities.

The global synchrony problem

Most existing simulation algorithms for cell-based multi-scale models rely on the composition of sub-models representing different processes occuring at different scales (Sü erlin et al.

). For example, in our simulations of TRAIL treatments on multi-cellular spheroids (Chapter ), TRAIL diffusion, cell motion and intra-cellular apoptosis signaling were not simulated in a completely simultaneous fashion.

Rather, a global time step was defined, for which those processes were simulated sequentially, and shared variables between processes (for example TRAIL concentration from diffusion and TRAIL concentration as seen by the intra-cellular model) are updated only discretely between each global steps. As a result, simulated dynamics are only approximating the true theoretical dynamics, and the global time step should be taken small enough such that this approximation remains acceptable (we used a global time step of minutes for spheroid TRAIL treatment simulations).

In some circumstances, the timescales at which shared variables evolve can be such that a very small global time step should be used to ensure sufficient accuracy, which can result in excessively long computation times. For example, when a antibody with very high affinity to some cell surface receptors is released in a tissue, the spatial spread of the antibody would be governed by a 'saturate then diffuse' mechanism that will require a global time step smaller than the receptor binding timescale to be correctly reproduced (e.g. molecules that should be bound at the first cell layer will be incorrectly moved to the next cell layer during the step computing diffusion). To address such issues, one might want to construct an 'asynchrony error' measure such that it can be controlled via an adaptive global time step procedure. But rigorously choosing such measure is not straigthforward and might depend on the specifics of the model.

A radical alternative would be to abandon the 'composition of sub-models' simulation paradigm and unify the continuous dynamics of the whole model into a single set of equations that will be fed to a single solver. After all, mathematically, the equations representing protein-protein reactions, cell motion due to cell-cell contacts, the evolution of TRAIL concentration in a voxel because of diffusion, etc… are all differential equations. In such case, the simultaneity of all continuous processes is automatically accounted for.

One difficulty with such approach is that the overall system of differential equations changes each time a cell agent is created or destroyed. We adapted the code of the DOPR ODE solver from Numerical Recipes (a 5 th order Runge-Ku a based scheme with adaptative time step) and integrated it into multi-agent simulation code to allow for changing the differential equations each time such event occurs. We tried this approach for simple multi-scale models of the yeast pa erning system and obtained comparable results to the standard 'composition of sub-models' algorithm. We also rely on a similar approach for our tool allowing the automated simulation code generation from a programmatic model description (see next section).

Another potential problem with such approach is that the ODE solver will be fed with a system of equations of an unusually large size for which it might not be adapted. Moreover, we considered first only models with deterministic dynamics (except for agent events that can be stochastic), preventing for example to have stochastic gene expression modeled into each cell. Still, in principle this 'unified and synchronous simulation across agents' paradigm can be extended to models in which the resulting equations are hybrid stochastic/deterministic or stochastic differential equations, provided a suitable solver (for arbitrary hybrid models or arbitrary stochastic differential equations) is available.

Numerical convergence tests

Here we present numerical tests related to the choice of timesteps to simulate the repeated application of TRAIL on multicellular spheroids (Chapter ). With this approach, we were able to obtain functional simulation programs for the spe-cific models we were working on. However, this required a lot of time and efforts: designing, testing and debugging such very large (> K lines for the original code simulating cell-cell physics) and heterogenous code was extremely tedious. Moreover, the obtained simulation programs lack flexibility and re-usability: slight changes in the model assumptions would at best require error-prone manual modifications at many locations in the code, or in many cases involve significant re-writing if the structure of the simulation algorithm is impacted by the change. In addition to model analysis and simulation, the task of simply describing the model in an unambiguous mathematical form can be extremely tedious because of the large number of variables, parameters and equations. This is already true for single formalism models, such as large ODE models of signal transduction (for which rule-based modeling can be an effective solution, Hlavacek et al.

), but this problem is even more stringent in the context of cell-based multi-scale models that combine several mathematical formalisms and hierarchical levels.

To solve some of those issues, we developped a tool that enables the automated generation of simulation code from a programmatic description of the corresponding mathematical model (Figure . ). ). In this thesis, we initially reasoned that extending ODE-based models of biochemical pathways with stochastic protein fluctuation models for all its proteins could greatly extend the prediction scope of those models (Chapter ). More precisely, the temporal scope should be extended compared to approaches that account for cell-to-cell variability in protein levels via static distribution in initial conditions (Spencer et al.

; Gaudet et al. ; Kallenberger et al. ).

Importantly, we have hypothesized that a precise knowledge of all the biochemical parameters involved in the expression of a given protein is not needed, provided that the amplitude and speed of the resulting protein fluctuations are sufficiently accurate. We have proposed an approach allowing such systematic inclusion of stochastic protein fluctuation models (Chapter ) despite incomplete knowledge and have demonstrated its relevance on two examples. In the case of a large model of the TRAIL-induced apoptosis pathway, it enabled to correctly predict subtle measurements of how the cellular state of sister cells diverge from one another as time goes by because of stochastic fluctuations in the level of signaling proteins (Chapter ). In the case of a synthetic sensing circuit implemented in yeast cells, our approach permi ed to reproduce its induction dynamics with high accuracy, temporal scope and at the single-cell level (Chapter ).

Models of biochemical reactions that systematically combine stochastic gene expression with protein-protein reaction kinetics at the scale of a full signal transduction pathway (typically for more than proteins) are rare. Perhaps the closest example to our work is the model of N F -κB used to investigate the digital nature of this pathway activation upon stimulation (Tay et al.

; Kellogg & Tay ). In this work, stochastic gene expression is modeled for several genes which products feedback on signaling dynamics. However, for several other proteins, their cell-to-cell variability is modeled via static distributions in initial conditions, hence se ing a limit to the temporal scope of the model of the order of the mixing time for those proteins, as was the case for the ODE-based model of TRAIL-induced apoptosis that we extended in Chapter . ) or resistance to a drug (Nevozhay et al.

Linking extended biochemical pathway kinetic models with cellular decisions to simulate cell population dynamics

). Also, in E.coli, fluctuations of a metabolic enzyme in conditions such that its level is highly limiting were shown to generate fluctuations in cell growth, as measured by cross-correlation analysis between expression level and instantaneous growth rate in single cells (Kiviet et al. ).

What In the context of TRAIL-induced apoptosis, we have found that limited drug penetration within tumor spheroids can exacerbate non-genetic resistance acquisition as predicted for cell populations facing the same drug concentration but is unlikely to cause by itself long-term resistance to repeated treatments. We have focused on a single dimension of heterogeneity in cell micro-environment (drug concentration) but many other environmental factors might play a role, such as mechanical cues that can impact gene expression of many genes (that might include genes involved in the response to TRAIL) via mechanotransduction. For example, it has been found that E-cadherin expression, which is known to be mechanically regulated, correlates with the expression of TRAIL death receptors (Lu et al. ).

In the context of the engineering of synthetic spatial systems, we have shown that our approach could enable the obtention of highly realistic, single-cell resolved predictions of the spatio-temporal response of micro-colonies bearing a sensing circuit and responding to dynamic gradients of messenger molecules (Chapter ). More traditional modeling approaches used by synthetic biologists (such as ODE models for gene circuit dynamics that are extended into PDEs to represent space (Tabor et al.

) or integrated into spatial cellbased models (Gorochowski et al.

)) would have failed to reproduce the spatial heterogeneity in the micro-colonies response.

.

Limitations and perspectives

Limitations of the single-cell models

One limitation of the single-cell models we used in this thesis (for TRAIL-induced apoptosis and sensing in yeast) is that neither cell size nor progression to the cell cycle is represented.

As a result, when we model molecule copy numbers, it does not account for the fact that copy numbers usually increase with cell size (Marguerat & Bähler ). Thus our model represents 'effective' copy numbers rather than real copy numbers, and as a consequence when modeling division we do not split randomly the content into the two daughter cells but rather copy the state. This was justified in the case of TRAIL-induced apoptosis as ) death time distributions were not dependent on cell cycle stage and ) recently divided sister cells had very high correlation of death times, showing that spli ing noise does not have a strong contribution (Spencer et al.

). However, it is possible to explicit the dependency with growth and the effect of spli ing, as was done for a small network (Volfson et al.

).

It might be required in conditions where growth can be highly variable between cells and is regulated by one or several molecular factors represented in the model.

Another limitation to those models is that despite their level of detail, many steps are not represented or highly simplified. For example, this is the case for receptor clustering (not modeled) and caspase-processing at the DISC (highly simplified) in the kinetic model of TRAIL-induced apoptosis we used (Ho & Harrington ; Dickens et al. ; Neumann et al. ). This might or might not be a problem depending on the purpose of the model and the capacity of the simplified representation used for a given set of reactions to realistically capture the input/output relationship of this set of reactions.

Another concern that apply to our modeling of TRAIL-induced apoptosis relates to the relevance of simulating cell populations for many generations with a model that represents only a single pathway with molecular detail: as time goes by, the risk that cross-talks with other pathways starts interfering non-trivially with the modeled pathway to determine cell phenotypes increases. Our response to this concern is that it is anyway necessary to take this risk in order to push our understanding forward, but that it might be alleviated by adequate modeling approaches.

Combining mechanistic and phenomenological models into coarsegrained whole-cell models Indeed, while efforts spent into building whole-cell models that aim to represent as much molecular detail as possible are promising (Karr et al. ), the field of multi-cellular systems biology that address problems such as human diseases or animal development should not wait for the obtention of such models for human or animal cells before a empting to gain a molecular-level, quantitative understanding of the behavior of complex populations evolving in complex environments (such as developping tumors treated by anti-cancer drugs).

This imply to make choices regarding which aspects of cellular environmental states should be understood with molecular resolution and which aspects can be represented phenomenologically via abstract sub-models. The resulting models will thus necessarily mix different levels of detail/abstraction. For example, in our simulations of repeated TRAIL treatments on multi-cellular spheroids, the decision to enter apoptosis is modeled with molecular details while cell growth and division are modeled in a very phenomenological and simple manner. Such models will maybe fall short on explaining some observations, motivating either an improvement of the abstract sub-models (for example, accounting for memory effects in the cell cycle duration by using stochastic processes instead of sampling a cell cycle duration in the same distribution at each division) or the replacement of abstract/phenomenological sub-models with mechanistic, molecular level representations.

While such models might seem un-elegant and somehow messy, they will probably be necessary for realizing the potential of multi-cellular systems biology to quantitatively explain complex biological phenomena with molecular detail.

Recently, a coarse-grained whole-cell model of E. coli has been proposed (Weisse et al.

). This model has been constructed to account for three fundamental cellular trade-offs that link gene expression, growth and metabolism: making proteins costs energy, mRNAs compete for translation, and the total proteome is constant (in steady-state growth). It correctly predicts empirical relationships between growth rate, nutrient levels and the ribosomal protein fraction. Thus, this model could be used as a 'chassis' into which mechanistic, molecular-level models are integrated in order to realistically simulate how a specific process interacts with the global physiology of cells. However, because this model abstracts progression within the cell cycle and cell division, a 'single-cell' version should be constructed first to enable the meaningful integration with single-cell models of biochemical reaction pathways.

Using synthetic circuits to study the population consequences of cellular decisions

Building quantitative models in order to realistically simulate the dynamics of cell populations with molecular detail is still very difficult and is currently rarely a empted. Much more progress is needed to establish methods and concepts that are robust enough to make such modeling rewarding in terms of biological insights. We believe that multi-cellular synthetic biology can play a decisive role towards this goal. This could be realized by the combined experimental and theoretical investigation of how synthetic gene circuits (that are designed, well-characterized and tunable) can interact with, and eventually control, the dynamics of cell populations. The gain for systems biology research could be double: in addition to establishing modeling methods and concepts, it is likely to provide insights on the 'design principles' by which natural gene circuits implement various, complex and robust functions at the population level.

  which reaction-diffusion equations shed light on the mechanisms of pattern formation), the cell-based multi-scale modeling paradigm (Thorne et al.

  tain cases without requiring cell destruction, therefore enabling live-cell observation. Key elements in the emergence of single-cell experimental biology are genetic engineering, i.e. the ability to manipulate the genome of cells; the discovery of fluorescent proteins (Chalfie et al. ) and their utilization as tags for measuring the level or activity of any protein of interest (Crivat & Taraska ; Regot et al. ); flow cytometry and live-cell microscopy (Selimkhanov et al. ); the development of microfluidic cell culture systems (Mehling & Tay ), that enable very precise control on the micro-environment in which cells are grown; and more recently, 'omics' techniques that allow large-scale, close to genome-wide measurements (Wang & Bodovi ): single-cell mass spectrometry (Bendall et al.

  tions. Historical examples are the variability in the number of viral particules released by infected bacteria, that can be mainly a ributed to stochastic production of viral particules in the infected cell (Delbrück ); or the variability in how much and how fast cells produce a protein in response to promoter activation, to which instrinsic stochasticity of transcription and translation largely contributes (McAdams & Arkin ). Other examples include the contribution of a triple positive feedback loop in N F -κB signaling to digital pathway activation (Tay et al.

  tool to infer parameters of a protein fluctuation model from the protein variability and mixing time, and tools that allow the automated generation of simulation code from a programmatic description of a model). Those tools could be useful to other researchers interested in the simulation of cell-based multi-scale models and are available on GitHub http://github.com/fbertaux. . Outline This document is organized as follows. Chapter describes the construction and validation of a single-cell model of TRAIL-induced apoptosis. The model is constructed from a previous ODE model using a generic and principled approach to systematically account for protein level fluctuations. Model predictions are compared to existing data on cell fate variability, cell fate correlation between sister cells, and transient resistance gain after treatment.In Chapter , the obtained single-cell model of TRAIL-induced apoptosis is integrated in a cell-based multi-scale model to predict the long-term response of cancer cell populations to repeated TRAIL treatments in two virtual experimental in-vitro se ings: 'monolayer' and spheroid cultures. In the spheroid se ing, the fate of cells are coupled by reduced diffusion and increased degradation of TRAIL within the spheroid. In the 'monolayer' se ing, no direct or indirect cell-cell couplings are assumed to play a role. We use first the 'monolayer' se ing to investigate the role of time in TRAIL resistance. We vary the treatment parameters (dose and period) as well as parameters relating to the molecular implementation of TRAILinduced apoptosis (mimicking cell lines with different resistance properties). The results provide insights on the mechanisms driving long-term resistance that we formalize in a simple phenomenological model. Finally, we use the spheroid se ing to investigate the additional impact of limited TRAIL penetration into the long-term resistance to repeated treatments.

  the importance of protein fluctuations in TRAIL-induced apoptosis TRAIL: an inducer of extrinsic apoptosis selective towards cancer cells Apoptosis is a form of controlled cell death playing a crucial role in the development and homeostasis of multi-cellular organisms. Many human diseases are linked to anomalies in the regulation of apoptosis. Apoptosis can be triggered either from inside the cell (intrinsic apoptosis) or by external death ligands, i.e. proteins that are able to bind certain type of membrane-bound proteins called death receptors (extrinsic apoptosis).Several types of death ligands and death receptors are known. Among them, TNF (Tumor Necrosis Factor) and its receptors (TNFRs) were the first to be characterized. TNF plays a central role in the regulation of inflammation and many successful drugs target its receptors. Other death ligands are CD L (and its receptor CD ) and TRAIL (and its receptorsDRs) (Ashkenazi & Dixit).When TRAIL was discovered in the mid-'s (Wiley et al. ; Pi i et al. ), it raised high interest because it was observed to induce apoptosis preferentially in cancer cells while being harmless to healthy cells (Walczak et al. ; Ashkenazi et al.

Figure 2 . 1 :

 21 Figure 2.1: TRAIL-induced apoptosis. A simplified view of the TRAIL-induced apoptosis pathway is shown.The displayed interactions are abtracted from an established kinetic model(Albeck, Burke, Spencer, et al. 2008).

  proach to account for protein fluctuations in signal transduction pathways Extended vision of signal transduction pathways Protein synthesis and degradation are subjected to noise, resulting in fluctuations of protein concentrations in individual cells and in cell-to-cell variability at the population level (Raser & O'Shea ). Such variability could have consequences on signal transduction: aside of conventional epigenetic differences (Rando & Verstrepen ), unequal access to ligand molecules or simply noise in signaling reactions, it often contributes importantly to heterogeneous behavior within an isogenic population (Spencer et al. ; Colman-Lerner et al. ; Geva-Zatorsky et al.

Figure 2 . 3 :

 23 Figure 2.3: Accounting for stochastic protein turnover in signal transduction pathways. Scheme of the modeling approach. Protein-protein interactions mediating signal transduction (signal transduction layer) are modeled by ordinary differential equations. In parallel, promoter activity changes, mRNA production and degradation (gene expression layer) are seen as stochastic events and generate fluctuations in mRNA levels. This impacts the synthesis rates of the corresponding proteins. Together with protein degradation, it generates fluctuations in protein levels (here shown in absence of transduction). Only a fragment of the extrinsic apoptosis pathway is shown. Deterministic/stochastic interpretation of chemical reactions is represented with black/red arrows respectively.

Figure 2 . 4 :

 24 Figure 2.4: Protein level auto-correlation function shape depends on model parameters. Examples of autocorrelation curves for three different parameterizations of the stochastic protein turnover model. In each case, the corresponding mixing time τ , i.e. the time at which auto-correlation equals 0.5, is indicated.

Figure 2 . 5 :

 25 Figure 2.5: Fluctuations of long-lived proteins are less dependent on transcriptional bursting kinetics than short-lived proteins. Comparison of protein level coefficient of variation (A) and mixing time (C) as a function of transcriptional bursting rates for two settings: a short-lived protein and mRNA (half-lives of 2 and 1 hours, resp.) and a long-lived protein and mRNA (27 and 9 hours, resp.). Combinations of Ton and Toff values ranging from 0.1 to 5 hours and 0.1 to 10 hours respectively were tested (T on = 1 k off and T off = 1 kon are mean ON and OFF time of the gene). In each case, the remaining rate parameter k sm was chosen such that E[m] = 17. (B and D) Representation of the range of values obtained for all models tested in (A) and (C).

Figure 2 . 6 :

 26 Figure 2.6: Behavior of a standard stochastic protein turnover model. Promoter switching rates respect typical ranges observed in (Suter et al. 2011) and lead to a protein level CV of 0.25. Left plots show three representative single-cell time courses of protein and mRNA levels. Histogram at the right displays the corresponding distribution of protein level obtained when simulating a large number of cells for a long duration, corresponding to a snapshot of the cell-to-cell variability expected in a population.

  Model predicts transient memory in cell sensitivity to TRAIL and CHX Using live-cell microscopy, Spencer et al. investigated the fate of hundreds of cells after exposure to TRAIL and CHX ( ng/mL and . µg/mL, Figure . -A). All cells undergo MOMP with a highly variable delay (from to hours, Figure . -B). To study cell fate inheritance, the authors also recorded hours before treatment to identify sister cells (Figure . -A). They were found to have highly correlated MOMP times (correlation coefficient close to for recently divided cells, about . for older sisters -Figure . -B, black curve). Here, the MOMP time distribution provides a quantification of the cell fate variability, while MOMP time correlations between sister cells also give information on the transient memory in cell state. Within our framework, in-silico reproduction of those experiments is straightforward (Figures . -D and . -C), enabling us to investigate possible origins of transient cell fate inheritance.

Figure 2 . 8 :

 28 Figure 2.8: Cell fate variability in TRAIL-induced apoptosis. (A-C) Cell fate variability experiments performed in (Spencer et al. 2009). (A) HeLa cell populations were treated with either 10 ng/mL of TRAIL and 2.5 µg/mL of cycloheximide (CHX) or 250 ng/mL of TRAIL alone. Cells were tracked during 8 hours by live-cell microscopy and MOMP time was detected via mitochondrial release of a fluorescent reporter. (B-C) Histograms of MOMP times and surviving fractions observed for treatment with (B) TRAIL and CHX or (C) TRAIL alone. (D-F) In-silico reproduction of those experiments with our "fitted" model (i.e. the parameterization in the explored parameter space region giving the best agreement for cell fate variability data, see detailed description in Appendix 2). (D) Simulations (see Appendix 3 for details). (E-F) Results for the (E) TRAIL and CHX or (F) TRAIL alone treatments. For the latter case, representative model trajectories are given in Figure 2.14.

Figure 2 . 9 :

 29 Figure 2.9: Transient cell fate inheritance in TRAIL-induced apoptosis. (A-B) Experiments measuring correlation of MOMP times between sister cells performed in (Spencer et al. 2009). (A) HeLa cells were recorded from 20 hours before treatment as in Figure 2.8-A. Sister cells were identified to permit comparison of their fate. (B) Quantification of cell fate inheritance was realized by computing the correlation between sister cells MOMP time as a function of the duration between division and MOMP (averaged between sisters). (C-D) In-silico reproduction of those experiments with the model of Figure 2.8. (C) Description. (D) Quantification of cell fate inheritance was applied to simulation results as in (B). See Appendix 3 for details.

  tained (Figure . -E). Further analysis revealed that such agreement requires FLIP and Mclprotein half-life to be short and to fall within a narrow range (between . and . hours, Figure . -C). This model prediction is consistent with previous measurements in HeLa cells ( and minutes for FLIP short isoform and Mcl-respectively, (Nijhawan et al. ; Poukkula et al. )). In contrast, FLIP and Mcl-mRNA half-life and promoter switching rates are not strongly constrained, probably because their influence on cell fate is limited by the rapid protein level decrease caused by synthesis blockade. We then asked whether our extended model also capture transient cell fate inheritance (Figure . -B). It is the case: fi ed models accurately predict the MOMP time correlation between sister cells (Figure . ,D,black curve and Figure . -C). Of note, assuming standard promoter switching rates for FLIP and Mcl-(but accounting for their short mRNA and protein half-life -this parameterization will later be referred as the "non-fi ed" model) already provides a good agreement for both MOMP time distribution and MOMP time correlation between sister cells (Figure . , upper plots, +CHX data). This non-trivial result shows that the speed at which the sensitivity to TRAIL and CHX fluctuates in single cells is well captured and thus suggests that our generic approach permits to describe fluctuations of protein levels with sufficient accuracy. . Model explains fractional killing and predicts transient cell fate inheritance Spencer and colleagues repeated this experiment but treated cells with TRAIL alone ( ng/mL). In this condition, an important fraction of cells died fast (MOMP in ~ hours) but % were still alive after hours (Figure . -C), illustrating the fractional killing property.

  ->MOMP) avg (h)> MOMP time (hrs) Frequency In uence of Flip and Mcl1 protein turnover model rates on Model-Data agreement in the + CHX condition

Figure 2 . 10 :Figure 2

 2102 Figure 2.10: Stochastic protein turnover models capture fluctuations of cell sensitivity to TRAIL and CHX. (A) Best agreement found between model and data for MOMP times distribution in the +CHX condition. Obtained for FLIP and Mcl-1 model rates such that protein/mRNA half-life and mean ON/OFF promoter activity duration equaled 0.4/1 and 1.9/3.1 hours respectively. See Appendix 3 for the description of how model data agreement was quantified. (B) Best agreement found between model and data for MOMP time correlation between sisters in the +CHX condition. Obtained for FLIP and Mcl-1 model rates such that protein/mRNA half-life and mean ON/OFF promoter activity duration equaled 0.3/1 and 0.35/24 hours respectively. (C) Influence of FLIP and Mcl-1 model rates on Model-Data agreement in the +CHX condition.For each parameter, we plot the model to data distance corresponding to the best model when all other three parameters are varied.

  Figure2.11: The "non-fitted" model quantitatively predicts TRAIL+CHX single-cell data and lead to fractional killing and reversible resistance for TRAIL alone treatments. In the non-fitted model, FLIP and Mcl-1 promoter switching rates are standard (Ton=0.1 hours and Toff=2.6 hours) but the short half-life of their mRNA and protein is accounted for (2 hours and 0.5 hours respectively). We reproduce here for this model all the results presented in Figures 2.8, 2.9, 2.16 and 2.17 for the "fitted" model. Quantitative agreement is obtained for TRAIL+CHX single cell data from(Spencer et al. 2009) (MOMP time distribution and sister cell MOMP time correlations). No quantitative agreement is obtained in the case of TRAIL alone treatments, but the existence of fractional killing and reversible resistance is nevertheless predicted. Note that because fewer cells were simulated compared to main text figures (5.10 4 instead of 10 5 for sister cell experiments), sister correlation curves appears slightly noisier.

Figure 2 . 12 :

 212 Figure 2.12: Model fitting to cell fate variability data predicts large, rare fluctuations of FLIP / Mcl-1 and transient cell fate inheritance. (A-C) Agreement between model prediction and experimental data for (A) death (i.e. MOMP) time distribution, (B) surviving fraction after 8 hours, and (C) both together, for treatment by TRAIL alone (250 ng/mL), as a function of FLIP / Mcl-1 promoter switching times (other parameters as in TableS1). (D) Representative protein level fluctuations of Mcl-1 described by a stochastic protein turnover model allowing good agreement for both MOMP time distribution and surviving fraction. This model has been used for Figures 2.8, 2.9, 2.16 and 2.17. (E) Model-data agreement for MOMP time correlation between sister cells. For (A), (B), (C) and (E), agreement quality increases from red to green. The quantification algorithm is detailed in Appendix 3.

Figure 2 . 14 :

 214 Figure 2.14: Representative single-cell trajectories before and after TRAIL treatment for the "fitted" model. Trajectories for two dying and two surviving cells (after 12 hours of TRAIL treatment) are shown. T-marked arrows denote the time of TRAIL addition (250 ng/mL), D-marked arrows denote the time of death commitment (MOMP). mRNA (lower left of each panel) and native form protein levels (upper left of each panel) are shown for pro-caspase 8, Bid and Mcl-1. Levels of activated caspase 8, truncated Bid and activated caspase 3 are also shown (upper right of each panel), as well as the ratio of released Smac and of cleaved PARP (lower right of each panel).

Figure 2

 2 Figure 2.15: Large, rare fluctuations of Mcl-1 alone are sufficient to explain cell fate variability and transient inheritance in both conditions. While FLIP and Mcl-1 protein and mRNA half-lives were the same as for the "fitted" model (0.4 and 1.0 hours respectively), only the Mcl-1 promoter was assumed to have low switching rates (Ton and Toff are 16 and 24 hours resp.). The switching rates of the FLIP promoter were assumed to be standard (Ton=0.1 hours and Toff=2.6 hours). All the results presented in the main text for the "fitted" model are reproduced here. Note that because fewer cells were simulated compared to main text figures (5.10 4 instead of 10 5 for sister cell experiments), sister correlation curves appears slightly noisier.

  Figure . -A,B). Future survivors strongly stood out by their Mcl-protein level and gene activity state (Figure . -B). FLIP also appeared to play an important role in determining cell decision, and smaller but significant effect was also seen for Bid, Bax, Bcl and XIAP. Although it is a good predictor of cell fate, initial Mcl-gene activity status does not completely determine survival: neither all Mcl-"ON" cells survived nor all Mcl-"OFF" cells died. Thus, preexisting differences in protein levels and promoter activities are major determinant of cell fate but stochastic events in gene expression occurring during signal transduction also play a role. While timing of death for cells treated with TRAIL and CHX appeared to be multifactorial (Spencer et al.), our results suggest that cell survival is predominantly

Figure 2 . 16 :

 216 Figure 2.16: Reversible resistance in repeated TRAIL treatments. (A) Schematic description of the 'repeated TRAIL' experiments performed in [14,15] to characterize reversible resistance in HeLa cells. (B) In-silico reproduction of these experiments with our model (details in Appendix 3). (C) Resistance gains in surviving cells relative to naïve cells as a function of time between the two TRAIL treatments. Data are shown for experimental observations (Flusberg et al. 2013; Flusberg & Sorger 2013) and model predictions (our study).Comparison of resistance gain experimental measurements between the two TRAIL doses (500 and 50 ng/mL for(Flusberg et al. 2013) and(Flusberg & Sorger 2013) respectively) should be done with care, as the measurement and quantification method differed.

Figure 2 . 17 :

 217 Figure2.17: Molecular determinants of cell fate and resistance to repeated TRAIL treatments. (A) Cartoon illustrating that the determinants of cell fate and resistance can be studied by analyzing the overrepresentation of protein levels in 'future survivors' (cells that will still be alive after treatment) at the time of treatment, and in surviving cells at day X, respectively. (B) Cell fate determinants analysis: overrepresentation (compared to initial population) of protein level (blue) and promoter activity (red) at the time of treatment in 'future survivors'. Asterisks mark differences that passed a 5% significance test. (C) Resistance determinants analysis: over-representation of protein levels in surviving cells at day X. Inset illustrates the recovery kinetics expected from protein turnover only (i.e. in absence of significant selection effect or residual signaling activity). Therefore, deviation from such kinetics indicates the presence of a selection effect or residual signaling activity.

  Figure.-C, inset -see the death receptor (R), pro-caspase ,Bar and Bid). Deviation from such kinetics indicates either the persistence of signaling reactions that continue to consume proteins (it is the case for Apaf, pro-caspase and XIAP, as further analysis confirmed) or is a consequence of important selection. Indeed, while Mcl-and FLIP should recover normal levels in a few hours in absence of selection because of their high turnover rate, Mcl-levels (but not FLIP levels) are still strongly higher than in naive cells one day after TRAIL treatment, consistently with previous observations on the relative selection strength that operated on them. Together, those results indicate that recovery phenotypes in surviving cells result from a complex interplay of three distinct effects: selection during apoptosis, transcriptional noise and protein turnover as a driving force tending to reset protein levels to their initial, pre-stimulus distribution, and long-term residual signaling activity. This explains why it is difficult to understand the recovery process and justifies the use of modeling to disentangle the various contributions. . Discussion Rehm et al. and Spencer et al. (Rehm et al. ; Spencer et al. ) made two insightful observations about TRAIL-induced apoptosis. First, recently born sister cells died almost

  of the Mcl-protein has the potential to rescue cells from TRAIL apoptotic signaling. In this context, solely accounting for protein fluctuations within the TRAIL apoptosis pathway predicts the fractional killing property (Figures . and .

  ). While our results challenge current opinion on the role of survival pathways in TRAILinduced apoptosis, they are consistent with observations made on wild type HeLa cells that neither blocking NF-κB response nor inhibiting the Akt pathway do significantly change the surviving cell fraction after TRAIL treatment (Braeuer et al. ; Lalaoui et al. ). The pivotal role for Mcl-in TRAIL-induced apoptosis predicted by our model is consistent with the recent finding that Mcl-silencing by shRNA in HeLa cells completely sensitize cells to TRAIL (Lemke, Karstedt, Abd El Hay, et al. ). While moderate fluctuations of Mcl-levels were sufficient to obtain fractional killing, a quantitative agreement with the Spencer et al. ( ) single-cell data (MOMP time distribution and surviving fraction) required large and rare Mcl-fluctuations, caused by rare switches between long periods of gene activity or inactivity. Interestingly, in that case, the observed rapid loss of MOMP time correlation between sister cells quantitatively emerged from model simulations. However, this result was obtained for a single TRAIL dose. Given the predicted importance of Mcl-in determining cell fate, it would be valuable to investigate model ability to reproduce surviving fraction, MOMP time distribution and sister cell fate correlation for lower and higher TRAIL doses. FLIP is often mentioned as a key factor in cell resistance to TRAIL (Shirley & Micheau ), but in our model FLIP has less impact on cell survival than Mcl-. Consistently, Lemke et al. ( ) silencing experiments demonstrated a dominant role for Mcl-and a synergy with FLIP. However, our model might under-estimate the role of FLIP: the representation of DISC-related events in EARM is simple and thus does not account for recent biological findings, including the stoichiometry between its components (Dickens et al. ; Schleich et al.

Figure 2 .

 2 Figure 2.19: Impact on resistance gain in one-day survivors of the presence/absence of caspase feedback loop and the degradation of non-native ('active') forms. In-silico repeated TRAIL experiment (as in Figure 2.16) was repeated for variants of the "fitted" model regarding presence/absence of the caspase feedback loop and the default active forms half-life. Resistance gain in one-day survivors is shown. Simulations were repeated 4 times with 10 4 cells, error bars indicate standard deviation of estimated resistance gain between replicates.

  towards a multi-scale understanding of TRAIL resistance In the previous Chapter, we developped a single-cell model of TRAIL-induced apoptosis that accounts for ) protein-protein signaling reactions linking TRAIL exposure to commitment to apoptosis, ) stochastic gene expression for the proteins involved in this signaling and ) protein degradation. An important result that was obtained is that under parsimonious and realistic assumptions for parameter values, fractional killing and transient resistance acquisition readily emerged from model simulations. Those two properties relating to TRAIL resistance are observed in-vitro for many different cancer cell lines. A priori, such resistance can have implications for the therapeutical usage of TRAILlike drugs (i.e. death receptor agonists, abbreviated DRAs), notably regarding the choice of the DRA (potency and pharmaco-kinetic properties) and the treatment strategy (dosing, scheduling and co-drugging for decreasing resistance). Also, in addition to the mechanisms of TRAIL resistance implicated in in-vitro experiments, other factors could impact the response of in-vivo tumors to TRAIL treatments, such as limited and heterogenous drug exposure, context-dependent changes in cell state affecting TRAIL sensitivity, influence of other cells such as endothelial or immune cells, etc. (Pampaloni et al. ; Yamada & Cukierman ).

(

  Figure3.1: Modeling populations of single cells submitted to TRAIL treatments. Our single-cell model of TRAIL-induced apoptosis is simulated in each cell of a population. In addition, cell proliferation is modeled by dividing cells after a random, normally distributed cell cycle duration. At division cell state is duplicated in the two daughter cells. Such population can be submitted to arbitrary time-varying TRAIL concentrations (for example generated by bolus administration at discrete times), leading to single-cell resolved population dynamics driven by cell proliferation and TRAIL-induced apoptosis.

  our single-cell model of TRAIL-induced apoptosis, because this parameterization readily predicts the phenomena of fractional killing and transient resistance acquisition. This allow us to investigate the long-term consequences of those phenomena without relying on too specific, yet not validated assumptions (such as the predicted Mcl-bimodal distribution and slow transcriptional bursting). Hence, we assumed that Mcl-and Flip transcriptional bursting rates have standard values (T on = 0.1 hours, T off = 2.6 hours, see Figure .

  cell populations treated with various TRAIL doses for hours resembles experimental observation. Of note, the dose-response curves for Flip shRNA and Flip + Mcl-shRNAs qualitatively differed from observations, as expected because the representation of Flip in the model is simplistic (not shown).

Figure 3 . 2 :

 32 Figure3.2: Predicted killing at 24 hours as a function of TRAIL dose resembles experimental observations on HeLa cells. Simulated dose-response curve (the principle of such simulations is detailed in Figure3.3) for a 24 hours end-point to reproduce experiments from(Lemke, Karstedt, Abd El Hay, et al. 2014). Inaddition, to allow comparison with the Mcl-1 shRNA silencing experiment performed by Lemke et al, we also performed simulations in which Mcl-1 synthesis rate was reduced by a factor 10. Markers (circles for WT, squares for Mcl-1 shRNA knock-out) indicate simulated (red) or experimental (green) cell count/viability 24 hours after treatment. Solid lines indicate corresponding sigmoid fits. Because the TRAIL formulation (different strategies are used to enhance the cross-linking capabilities of the original TRAIL protein) used by(Lemke, Karstedt, Abd El Hay, et al. 2014, isoleucine zipper-tagger TRAIL) differs from the one used in the Sorger group ('SuperKiller' TRAIL), we introduced a scaling factor that was chosen such that sigmoid fits of the simulated and observed dose-response curve (in absence of silencing) agree the most (here, one unit of 'Sorger group' TRAIL corresponds to 17.7 units of 'Lemke' TRAIL).

  Figure 3.3: Simulating the TRAIL dose-response of populations of single cells. Populations of growing cells were submitted to different TRAIL concentrations when they reached a size of 5000 cells. The resulting cell number curves are shown (two independent simulations per dose were performed). Note that the TRAIL concentration sensed by each population decreases exponentially with time because TRAIL is unstable (9 hours half-life). From such curves, one can extract a standard cytotoxicity dose-response curve by dviding the number of cells after 72 hours by the cell number expected in absence of TRAIL. Black thin lines indicate standard deviations for the two replicates. The horizontal lines indicate either the no killing level or the killing level required to balance growth.

  such non-genetic transient resistance acquisition was observed experimentally (Flusberg et al. ) and is an emerging property of our model of TRAIL-induced apoptosis (see Chapter ). By performing longer in-silico experiments, in which TRAIL treatments are repeated every hours, we can investigate the long-term consequences of non-genetic resistance acquisition as predicted by our model for timescales that are close to in-vivo timescales for TRAIL-based anti-tumor assays. From such simulations, the temporal evolution of cell number curves informs about the long-term cytotoxic efficiency of the repeated TRAIL treatments (Figure . ).

  TRAIL treatment (days)

Figure 3 . 4 :

 34 Figure3.4: Cell populations repeatedly treated with TRAIL acquire resistance. The long-term evolution of population sizes when TRAIL treatment is repeated every 72 hours is shown. Simulations in which population size reached 15000 cells were stopped. Although the highest doses allows for strong reduction of the population size, resistance acquisition is clearly visible: later treatments are always less efficient than the initial treatment. Treatment times are indicated with vertical black lines. A log-scale is used such that constant killing efficiency among each consecutive treatment translates into linear progression between the peaks. Note that as expected, variability and noise in those curves increases as cell number decreases.

  strong impact on the long-term outcome of repeated TRAIL treatments. But what are the mechanisms driving this resistance acquisition? Two distinct types of contribution to resistance acquisition could be conceptually distinguished.First, the resistance to TRAIL of individual cells could increase because of molecular changes induced by previous TRAIL treatments. Note that the notion of "TRAIL resistance" for individual cells is probabilistic: for example, a cell having a probability of % to survive the next treatment (because of its current molecular state) could indeed survive; and the probability to survive another treatment could have rised to % because of TRAIL-induced changes. We term such 'adaptation' of individual cells to subsequent treatments cell-level resistance acquisition.

FigureFigure 3 . 6 :

 36 Figure . we show the temporal evolution of population averages of mRNA levels for cells repeatedly treated with different doses of TRAIL. Those curves clearly indicate selection effects, with, as expected, transcripts coding for anti-apoptotic proteins (Mcl-and XIAP

Figure 3 . 7 :

 37 Figure 3.7: Moderate selection effects revealed by changes in average mRNA levels in the population.Evolution of average mRNA levels in the population for long-term repeated (every 3 days) TRAIL treatments, for 3 different TRAIL doses (in ng/mL). Only R, pC8, Flip, Bid, Mcl-1 and XIAP are shown. vertical thin black lines indicate the times of TRAIL treatments. Note that the y-axis scale only range from 14 to 20 mRNA copies. In absence of TRAIL, an average of 17 is observed as expected (horizontal dashed line). Deviations above/below 17 indicate positive/negative selection of cells with respect to that mRNA (and hence the protein it encodes).

Figure 3 . 8 :Figure 3 . 9 :

 3839 Figure 3.8: Evolution of average protein levels in the population suggests that cell-level, TRAIL signaling related degradation plays an important role. Evolution of average native protein levels in the population for the same in-silico experiments as Figure 3.7. Only R, pC8, Flip, Bid, Mcl-1 and XIAP are shown. Vertical thin black lines indicate the times of TRAIL treatments. Expected levels in absence of TRAIL are indicated by the horizontal dashed line. Deviations above/below those levels is caused either by selection effects or cell-level TRAIL signaling related degradation.

  Figure3.10: Increasing non-native forms degradation equally decreases both initial and long-term TRAIL resistance. In all plots, different colors correspond to different values for the aggregate parameter describing the degradation of non-native forms. Top: cell number curves for 1000 ng/mL TRAIL treatments repeated every 3 days. In each case, two independent replicates are shown. From those curves, the initial and longterm killing efficiencies can be computed. Bottom left: Long-term versus initial killing efficiency. Bottom right: ratio between the initial and long-term killing efficiencies ("resistance factor") versus the initial killing efficiency. The grey dashed lines highlight the apparent linear relationship between initial and long-term killing efficiencies.

Figure 3 . 11 :

 311 Figure3.11: Turnover rate of stable proteins governs resistance acquisition. Same analysis as in Figure3.10 except that the turnover rate of stable proteins is varied (i.e., dilution rate). The fold change corresponding to the reference model is indicated in the colorbar (1 corresponds to a half-life of 27 hours, 10 to a half-life of 2.7 hours, etc.). In each case, the corresponding protein synthesis rates are scaled by the same factor such that the mean expression level is unchanged. Grey dashed lines indicate the linear relationship between initial and long-term killing efficiency that was valid for variation in non-native forms degradation (see Figure3.10). It does not hold here, instead, resistance acquisition is strongly affected while initial resistance does not change much.

Figure 3

 3 Figure 3.12: A 3-parameters phenomenological model to describe the population-level killing dynamics in presence or absence of death signals. A single variable represents the instantaneous ability of the population to transform death signals into cell death (transmission capacity, noted G). G is simply the death rate that is obtained in presence of death signals. Resistance acquisition is modeled via the parameter k, which represents the cost of transmission in terms of transmission capacity: in presence of death signals, killing occurs but is consumed. If k = 0, there is no resistance acquisition. The 'synthesis' of G and its normal (i.e. in absence of death signals) 'degradation' are modeled with the two parameters σ and γ respectively.

  Figure3.13: Resistance factor dependency on treatment period of the mechanistic model can be reproduced with the 3-parameters phenomenological model. For the mechanistic model, repeated treatments with 1000 ng/mL TRAIL were performed for different treatment periods. In each case, the resistance ratio (killing efficiency of late treatments divided by killing efficiency of the first treatment, where killing efficiency is the ratio between cell number at the time of next treatment and cell number at treatment time) was computed. Parameters of the phenomenological model were adjusted to reproduce the same dependency, assuming that each treatment leads to 3 hours of death signal presence. Those parameters are σ = 1, k = 2 and γ = ln 2 55 hrs -1 . For the mechanistic model, increase of the resistance factor at low TRAIL periods is expected, as delays between TRAIL exposure and cell death result in death caused by the first treatment being attributed to the second treatment, etc.

First

  , to describe the motion of individual cells, we model the physical forces exerted on them and solve the corresponding equations of motion. Three types of forces are considered: adhesive/repulsive forces for cells in contact, micro-motility forces to mimic the ability of cells to explore their environment by active migration, and friction forces to account for the viscosity of the extra-cellular matrix.

Figure 3 .

 3 Figure 3.14: Multi-scale model of multi-cellular spheroids treated with TRAIL. Cells are physically represented as growing spheres.Cell motion is governed by adhesive/repulsive forces for cell-cell contacts, micro-motility and a friction term. Inertia is neglected because motion is friction-dominated. The rate at which cell radius grow is chosen to match a normally distributed cell-cycle duration as for previous section. The diffusion and degradation of TRAIL is also solved. Diffusion is assumed to be small when cell density is high, degradation is assumed to be high when cell density is high. In each cell, an instance of the singlecell model of TRAIL-induced apoptosis is simulated and governs cell decision to enter apoptotis. Apoptotic cells have their radius decreasing at a constant rate.

  described the mathematical equations describing the model but did not completely describe how those equations are numerically solved. Briefly, the three processes (cell growth and motion, extracellular TRAIL diffusion and degradation, intracellular TRAIL signaling) are simulated asynchronously with a global timestep of minutes, and each process has its own internal timestep (the timestep for cell growth and motion is computed adaptively based on maximum displacement, the timestep for extracellular TRAIL diffusion and degradation is computed from the CFL criterion, the timestep for intra-cellular TRAIL signaling is controlled by the ODE solver based on absolute and relative error tolerance parameters). The numerical aspects related to the simulation of such models are discussed in more details in Chapter (section . ). In particular, tests of numerical convergence justifying the choice of the global timestep and the diffusion and degradation internal timestep are given (Figure.).

  Figure . . It displays different time points and cellular properties per time point: TRAIL local concentration, active caspase-levels and normal/apoptotic cell status. It is interesting to see how the 'non-noisy' cell-extrinsic information (TRAIL concentration) is transformed in a more noisy pa ern of caspase-activation and an even more noisy pa ern of cell fate decision because of cell-intrinsic heterogeneity and stochasticity.

Figure 3 .Figure 3

 33 Figure 3.15: Repeated TRAIL treatments on multi-cellular spheroids. Visualization of spheroid evolution under repeated TRAIL treatments (125 ng/mL every 3 days, starts when spheroid size reaches 8000 cells). Six time points are shown. In each case, full spheroid (upper-left), TRAIL penetration (upper-right, dark blue: 0 ng/mL, red: 100 ng/mL), caspase-8 activation (bottom-left, dark blue: 0 copies/cell, red: 400 copies/cell) and apoptotic cells (bottom-right, in red) are displayed. Radius reduction rate for apoptotic cells is set to 1.0 µm per hour.

  If it is experimentally invalidated, it would indicate that effects not captured by the model have an important role. A likely candidate would be a pleiotropic effect(Hart et al.) of TRAIL, which might have a positive influence on cell proliferation (Ehrhardt et al.) and/or the level or activity of anti-apoptotic proteins (through transcriptional, posttranscriptional or post-translational regulation) (Falschlehner et al.

  Mcland Flip complexes with tBid and R* are degraded at the same rate as Mcl-and Flip, those two proteins are effectively targeting tBid and R* for degradation). Explicitely modeling all the proteins involved in the targeted degradation of the core apoptosis signaling might be needed to gain a more accurate and complete picture of the regulation of TRAIL-induced apoptosis, and hence of TRAIL resistance acquisition. Indeed, a study has found that the half-life of Flip and Mcl-was markedly changed in lung cancer cell lines that acquired TRAIL resistance through weeks-long treatments with escalating TRAIL doses (Wang et al.

  model, both selection effects and cell-level increase of resistance in surviving cells are abstracted by a single parameter: the 'cost' that translating death signals into cell killing has in terms of the capacity to translate subsequent death signals. Interestingly, this model captures essential features of the mechanistic TRAIL-induced apoptosis model, such as the pre-

  towards spatial self-organization of cell populations Multi-cellular synthetic biology aims to extend the scope of synthetic biology by implementing tissue-level functions in living systems. This task is highly challenging because of the multi-scale gap between what should be modified (genes in cells) and what should be obtained (a desired behavior of the whole cell population).

"

  Quorum sensing" cell-cell communication (when cells can both synthetize and sense a messenger molecule) is in principle useful for implementing spatial organization, as it provides cells information about their social context. An early demonstration of such potential was brought by Basu & Weiss. They engineered E. coli to be either AHL "senders" (synthetizing and releasing AHL) or "receivers" (expressing fluorescent proteins in response to AHL) and obtained with this system geometric pa erns of fluorescence, which shapes depend on the spatial configuration of sender cell colonies. In that example as well as several others (Chen & Weiss ; Tabor et al. ; Payne et al. ), control over the population spatial organization is only exerted on the phenotypes of the cells (expression/secretion of certain molecules, etc.

  from simple local rules driving cell fate. A famous example is the mathematical object called the 'Game of Life' (Conway ). In the 'Game of Life', a D square la ice composed of cells that are either 'occupied' (i.e. alive) or 'un-occupied' evolves with time in a discrete and deterministic fashion: alive cells can die either from loneliness or over-crowding, and un-occupied cells can become alive if they have an adequate number of alive neighbors (to represent cell division). The fact that such simple rules can generate rich and complex dynamic pa erns depending on the initial configuration fascinated many scientists.Drawing inspiration from the 'Game of Life' rules, the Weiss lab at MIT proposed to construct a pa erning system in yeast. The high-level design for this system is illustrated in Figure . . It relies on combining density-sensing (the ability for cells to sense the local cell density) with cell death programs such that individual cells would die if they are either isolated or overcrowded while growing normally for intermediate local cell densities. A natural choice to implement density-sensing in yeast cells is the previously developped artificial cell-cell communication system based on the plant cytokine IP (Chen & Weiss).

Figure 4 . 1 :Figure 4 . 2 :

 4142 Figure 4.1: High-level design of a patterning system inspired from Conway 'Game of Life'. Engineered cells would be able to sense their local cell density and to die if it is either too low or too high (left). Only intermediate densities would allow the cell to grow. The functioning of this design for a given 2D spatial configuration of cells is illustrated (right). The color of each cell represents the local density it is sensing. Group of cells that would die because of 'overcrowding' or 'isolation' are indicated. Other cells would grow normally.

  (i.e. bearing only IP sensing module driving GFP expression, see Figure.) to construct a quantitative, single-cell model of IP sensing. Then, we integrated this model into a spatial cell-based multi-scale model reproducing solid culture experiments characterizing the response of receiver cell micro-colonies to gradients of IP. We used this model to ) infer the diffusion coefficient of IP in the medium and ) assess the ability of the knowledge gained from liquid experiments to predict the behavior of the same cells cultured as spatially organized populations (thereby assessing transposability from liquid to solid). Finally, we parsimoniously modified the model of IP sensing to model the behavior of the high threshold killing in response to IP. We used it to predict the evolution of pa erning cell micro-colonies for which high threshold killing was activated in response to IP, and compared the results with corresponding experiments.

Figure 4 . 3 :

 43 Figure 4.3: Genetic implementation of the yeast patterning system. Receiver cells, that only bear the IP sensing module driving GFP expression, are also shown. For patterning cells, the distinct modules (IP production, IP sensing, high-and low-threshold killing) are highlighted. In both cases, constitutive promoters are indicated in green. Other promoters are either inducible by copper (blue), galactose (brown) or by IP via the IP sensing module (orange). For the high-threshold killing module, strict repression of Cre synthesis is achieved by TetR and is relieved only when Dox is added to the medium. Therefore, in principle, in patterning cells IP production, high-threshold killing and low-threshold killing are active only if respectively copper, Dox or galactose are present in the medium. High-threshold killing is achieved by Cre-mediated excision of the essential gene pkc1 (which has been flanked with loxP sequences). Low-threshold killing is achieved by Bax-mediated apoptosis induction (in presence of IP, this induction is prevented via the production of the Bax antagonist Bcl-2). In both strains, endogenous genes sln1 and ssk1 have been deleted to minimize non IP-induced activation of Skn7. Those circuits were constructed by Jing-Jing Sun (Sun 2014) and used parts developped earlier by Ming-Tang Chen (Chen 2008).

Figure 4 . 4 :

 44 Figure 4.4: Experimental characterization of the IP sensing module behavior in liquid culture. Performed by Brian Teague at MIT. Receiver cells, a yeast strain in which an IP sensing module (see Figure 4.3) driving GFP expression has been integrated in the genome, were grown in liquid cultures containing various IP concentrations and their fluorescence was analyzed by flow cytometry at various time points after induction. This experiment was repeated for two different initial cell densities (OD = 0.05 and OD = 0.2).

Figure 4 . 5 :

 45 Figure 4.5: Experimental characterization of the IP sensing module behavior in solid medium and in presence of IP gradients. Performed by Brian Teague at MIT. A line of receiver cell micro-colonies was seeded using a cell printer (~5-20 cells per colony initially) on a polyacrylamide gel and a drop containing IP was deposited on the first colony, hence initiating IP gradients in the medium. The colonies were grown for 24 hours in a microscope chamber and brightfield and fluorescence images were acquired every hour.

Figure 4 . 6 :

 46 Figure 4.6: Experimental characterization of the high-threshold killing module behavior in solid medium.Performed by Brian Teague at MIT. Patterning cell micro-colonies were seeded using a cell printer on a agar gel. All gels also contained Phloxine B (a dye that accumulates in dead cells) and imaging was performed such that dead cells appear black. In a control experiment, the gel did not contain Dox, Cu, Gal or IP and the cells eventually formed large micro-colonies (right). To test the functioning of the high-treshold killing module, a gel containing Dox and IP was used. Strong killing occured between 14h and 38h and resulted in clusters of dead and alive cells.

Figure 4 . 7 :

 47 Figure 4.7: Modeling the IP sensing module behavior in cells. Top: diagram of a semi-mechanistic kinetic model to describe the behavior of the IP sensing module in cells. It abstracts IP signaling with 17 reactions(15 rate parameters) involving 8 species: IP, its receptor and the complex with IP, an inactive transcription factor, its complex with the IP-bound receptor, the activated transcription factor and finally GFP (for which a maturation step is considered). Corresponding reactions are given in Table4.1. They are assumed to follow mass-action kinetics. Bottom: extension of this model accounting for stochastic gene expression. Stochastic transcriptional bursting models are driving the synthesis of R, TF and uGFP as in Chapter 2 for the TRAILinduced apoptosis model. In addition, to account for TF* -mediated synthesis of uGFP, the transcription rate of uGFP mRNA is made dependent on TF* level.

  GFP degradation γ dil + γ GF P IP : R → ∅ IP:R degradation γ dil + γ IP :R IP : R : T F → ∅ IP:R:TF degradation γ dil + γ IP :R:T F T F * → ∅ TF* degradation γ dil + γ T F *

First, the

  single-cell model IP sensing is modified according to the molecular implementation of the high-threshold killing module (Figure . ): GFP is replaced by the recombinase Cre, a protein fluctuation model is added to describe the vital Pkc protein, a stochastic reaction describing Cre-mediated deletion of the pkc gene is added, and finally cell death is triggered when Pkc protein level goes below a lethal threshold.

Figure 4 . 8 :

 48 Figure 4.8: Cell-based multi-scale model of receiver cell micro-colonies submitted to IP gradients.Individual cells are modeled as growing and interacting spheres (similarly to the modeling of tumor cell spheroids in Chapter 3, except that parameters are adapted for yeast cells (E = 1 kP a, no cell-cell adhesion, R birth = 4 µm, we did not implement budding as it did not appear necessary at this stage, and motion is limited in 2D). Each cell is equipped with an instance of the single-cell model of IP sensing described above. In addition, 3D diffusion of the IP molecules deposited in the gel is simulated and the resulting spatio-temporal profiles of IP are read by cells based on their localization.

Figure 4 . 9 :

 49 Figure 4.9: Quantification of IP induction dynamics in liquid culture. Analysis of the experimental data characterizing the IP sensing module behavior in liquid culture (see Figure 4.4). Gating was performed based on FSC-H and SSC-H to remove debris or cell aggregates. Resulting distributions of cell fluorescence for various IP induction levels and for different times after inductions are shown. Corresponding population average time-courses (normalized to the uninduced time-course) is also shown.

  Figure 4.10: Estimating IP sensing signaling parameters. The semi-mechanistic ODE model of IP sensing (see Figure 4.7) was used to reproduce the population-level induction dynamics curves (data normalized to uninduced time-course, the 5 µM IP dose was excluded as explained in text). The Downhill simplex method was used (implementation from Numerical Recipes). Dilution rate and IP degradation rate were fixed based on previous estimations. Fitted parameter values are given in Table4.2. In the plot, data points are squares and continuous lines correspond to simulated induction curves with the fitted parameters.

Figure 4 . 11 :

 411 Figure 4.11: IP sensing model accounting for stochastic gene expression captures observed single-cell induction dynamics.In-silico reproduction of the single-cell induction dynamics and comparison to data. 10 thousands cells per induction dose were simulated for the same duration as in experiments and their fluorescence was stored at the same timepoints as in the experiment. The simulated data was scaled such that the mode of the uninduced distribution correspond to the experimental one, allowing direct comparison between model and data. Then the same binning was used to construct histograms. Corresponding parameters are given in Table4.3.

Figure 4 . 12 :

 412 Figure 4.12: Quantification of induction profiles for individual micro-colonies. Analysis of the experimental data characterizing the IP sensing module behavior in solid culture (see Figure 4.5). The fluorescence images were used for segmenting colonies. For quantification of induction dynamics per colony, only the outer rim of each colony was considered to avoid biases introduced by piling up at colony centers. Background fluorescence was substracted such that the measurement reflects fold-change of cell fluorescence with respect to un-induced colonies. Three colony induction profiles are shown and their locations indicated by letters (A is the first colony on which IP was deposited).

  Figure 4.13: IP sensing model calibrated from liquid culture data is consistent with solid culture data and provides an estimation of IP diffusion coefficient. In-silico reproduction of the solid culture data. Three different values of the IP diffusion coefficient are tested. The parameter values for the ODE model describing IP sensing are the one estimated from liquid culture data (Figure 4.10 and Table4.2). The simulated spatiotemporal profiles for the three colonies shown in Figure4.12 are compared to the corresponding data. Good agreement is obtained for D IP = 10 -4 mm 2 .s -1 . In simulations, for simplicity we used the IP concentration at the colony center as being representative of the whole colony response.

  Among the few new parameters needed to describe high-threshold killing, two are a priori of central importance: the efficiency of Cre-mediated excision of the pkc gene, and the tolerance of cells to low levels of the Pkc protein. Intuitively, the former would govern the efficiency of the commitment to cell death, while the la er would somehow set the delay between commitment and death: after excision, Pkc protein levels would decrease because of degradation and dilution, but if cells are tolerant to low levels of Pkc , a commi ed cell might survive a long time after deletion and might even divide.Results of simulations in which those two parameters were varied are shown in Figure.. First, as expected, overall killing efficiency increases with both increasing pkc deletion rate (higher commitment rate) and reduced tolerance for low Pkc protein levels. Second, clearly distinct shapes of death/alive clusters are obtained for different values of those

Figure 4 .

 4 Figure 4.14: In-silico reproduction of the spatio-temporal response of receiver cell micro-colonies. The cell-based multi-scale model was used to mimick the experiment. Images of predicted single-cell fluorescence in growing micro-colonies are shown and compared to data (colony A here). The gray colorscale used to display simulation results was chosen according to the experimental one to allow a fair comparison. Clustering of 'responders' or 'non-responders' seems to be qualitatively reproduced by the model. A slightly decreased cell-cycle duration (2.2 hours in average) was used as compared to the doubling time of receiver cells measured in liquid culture (1.73 hours). Adapting the model to precisely reproduce colony growth dynamics would require single-cell tracking to quantify cell-level growth and a quantification of the 3-dimensional shape of micro-colonies driven by piling up.

Figure 4 .

 4 Figure 4.15: Same as Figure 4.14 but for colony B.

Figure 4 .

 4 Figure 4.16: Same as Figure 4.14 but for colony C.

Figure 4 . 17 :

 417 Figure4.17: Predicted individual cell fluorescence within micro-colonies is spatially correlated. Temporal evolution of spatial correlation of simulated cellular GFP levels within micro-colonies. The Moran's I metric was used, with the inverse of cell-cell distances as spatial weights. In each case, a negative control is obtained by shuffling the GFP levels among cells in order to evaluate the statistical significance of the results. For times below 4 hours, the cell number per colony is too low to conclude to a spatial correlation.

(

  IP sensing, Figure . (top), (Chen & Weiss )) which is a central element of a pa erning system currently developped in the Weiss Lab at MIT (Figure . (bo om)). Towards model-based, quantitative, single-cell level characterization of gene circuits behavior using flow cytometry A first achievement was the construction of a single-cell model of the IP sensing and its calibration from flow cytometry data quantifying IP induction dynamics in receiver cells. Importantly, instead of comparing directly predictions of the single-cell model with observed induction distributions, an intermediate step in which an ODE version of the IP sensing model (Figure . , top) was fi ed to population averaged induction dynamics allowed to obtain a first estimation of IP signaling parameters (Figure . and Table . ). Once this ODE model was calibrated, we applied a parsimonious strategy for choosing the stochastic gene expression parameters of the single-cell model, similarly to what was done for modeling TRAIL-induced apoptosis (Chapter ).

  PKC1 gene deletion rate (hr -1 per Cre unit) Lethal threshold of PKC1 protein level (rel. to avg. expression)

Figure 4 .

 4 Figure 4.18: Predicting the shape of death clusters in micro-colonies submitted to high-threshold killing.An in-silico reproduction of the experiment presented in Figure4.6 was performed with the cell-based multiscale model for the patterning cells (in which only the high-threshold killing module is activated). Different pairs of value for the two main parameters of high-threshold killing (Cre-mediated pkc1 deletion rate and lethal threshold of Pkc1 protein levels) were tested. In each case, an image of a micro-colony 15 hours after seeding (and IP induction as it is present in the gel) is shown. Overall killing increases with both parameters, but for similar overall killing the clustering effect increases when the Pkc1 lethal threshold gets lower.

  discuss the more high-level problem of how to optimize the 'workflow' of computational biologists, notably in terms of code re-usability and safety regarding potential programming errors, and propose a tool for the automated generation of simulation code from a programmatic model description (programmatic means taking the form of a computer program, such as in PySB,(Lopez et al. )).. Simulating combined signal transduction and stochastic gene expression Cellular decisions are governed via biochemical pathways involving protein-protein reactions. Those pathways are traditionally modeled via ordinary differential equations. Therefore, an effort has been spent by the community of cell-based modeling towards the development of frameworks that enable the simulation of ODEs into each cell of a multi-cellular simulation, in such a way that the state variables of those ODEs can be linked to the extracellular part of the simulation, hence allowing a coupling between cells and their environment (Sü erlin et al. ; Starruß et al.

Figure 5 . 3 :

 53 Figure5.3: Schematic of a procedure allowing the parameterization of stochastic protein turnover models from constraints on the protein level fluctuations (i.e. the protein CV and the protein mixing time τ ) with additional constraints on the three rates γ p , γ m and γ g = k on + k off . The procedure is capable of determining whether the rate constraints are compatible with the τ constraint, such that the user can change them until they are compatible.

Figure 5 . 4 :

 54 Figure 5.4: Exploration of parameter sets sharing the same protein fluctuation statistics CV and τ . More than 10 5 sets of rate constraints (γ g , γ m , γ p ) were sampled in a log-space hybercube according to the Sobolev algorithm. In each case,the parameterization procedure described in Figure 5.3 was used together with the fluctuation constraints CV = 0.5 and τ = 35 hrs. The sets for which a solution exist are displayed in blue in different slices of the parameter space. In the slices involving only γ g , γ m or γ p , the sets for which no solution exist are displayed in grey.

IP

  sensing circuit implemented in yeasts (Chapter , Figure., bo om), we have modeled the effect of the transcription factor on its target gene by assuming that k sm follows a Hill function of the transcription factor concentration.

First, the globalFigure 5 . 5 :

 55 Figure 5.5: Numerical tests for choosing timesteps. Top: influence of the global timestep on cell number temporal evolution (left) and TRAIL gradient (right, characterized via the maximum and minimum TRAIL concentration felt by cells). Bottom: influence of the internal timestep for TRAIL diffusion and degradation, chosen as a fraction of the CFL timestep. To facilitate the comparison between simulations, we turned off stochastic gene expression in cells (replaced by equivalent constant protein synthesis).

!

  Figure5.6: Principle of the automated code generation procedure. A model is constructed by writing a small python program using functions of a model description module. Those functions allow to define agents and their properties and dynamical rules (see Figure5.7 for a more detailed example). From such python-encoded description of a model, c++ simulation code corresponding to that particular model can be generated using an other python module. The generated c++ code is ready to compile and can be easily modified to define custom state initialization and model state outputs.

Figure 5 . 7 :

 57 Figure 5.7: Example of programmatic model description in python. The model represents a population of yeast cells bearing the high-threshold killing module and the IP production module (see Figure 4.3 in Chapter 4) that grow in liquid culture. Events (deterministic or stochastic) and differential dynamics are defined using dedicated functions.

Figure 5 . 8 :

 58 Figure 5.8: Example of c++ simulation code generation from a programmatic model description. Portions of the c++ code generated from the python model description shown earlier (Figure 5.7) are shown.

  What can we gain by extending the temporal scope of of biochemical pathway kinetic models by accounting for protein level fluctuations? The main motivation is in fact to enable the simulation of the dynamics of population of cells submi ed to given environmental conditions with a level of detail that encompass the molecular scale. Indeed, the behavior of biochemical pathways ultimately governs cell growth, division, death and other cellular phenotypes. In such cases, fluctuations and memory in cells internal state arising from stochastic gene expression are expected to impact population dynamics, because selection effects will apply as soon as different cellular states confer different cell fitness in a given environment. Therefore, integrating extended biochemical pathway kinetic models into cellbased population dynamics models by linking biochemical variables with cellular decisions allows to investigate complex effects such as non-genetic selection operating on fluctuating cellular states. This was a key question we wanted to adress in the context of TRAIL-induced apoptosis (Chapter and ), as such non-genetic selection could be implicated in transient resistance acquisition observed after treatment(Flusberg et al. ). It could also have consequences for the quantitative functioning of the yeast pa erning system studied in Chapter .Since several years, the topic of non-genetic selection operating on fluctuating cellular states raised high interest. Theoretical studies considering a single gene which product impact fitness (Sato & Kaneko ; Mora & Walczak ) showed how the steady-state distribution of that product at steady-state is changed by selection effects. Experimental demonstration of such effects was obtained using synthetic circuits that exhibited random switching between low and high expression for proteins conferring higher metabolic capac-ity(Acar et al.

  is the contribution of selection effects into transient resistance acquisition of cancer cells to TRAIL? In Chapter , we integrated our TRAIL-induced apoptosis single-cell model into multi-cellular simulations to investigate the long-term response of populations repeatedly treated by TRAIL. We have found a strong increase of resistance in the long-term when comparing killing efficiencies of later treatments with first treatments. Interestingly, a detailed analysis revealed that transient resistance acquisition as predicted by our model involves both selection and 'adaptation' effects. Such 'adaptation' was caused by a cell-level decrease of pro-apoptotic protein levels relatively to anti-apoptotic protein levels, which was mainly mediated by targeted degradation of activated proteins and differential constitutive turnover between pro-and anti-apoptotic proteins. Thus, selection effects are always present but might be minor contributors to non-genetic resistance acquisition. It should be noted that this does not mean that modeling fluctuations was not necessary: 'static' cell-tocell variability models would have grossly over-estimated selection effects, hence masking the contribution of cell-level resistance acquisition. To our knowledge, this mechanistic explanation for transient resistance acquisition is novel, highlighting the potential of our modeling approach to provide original insights. Another important modeling work that relates to the question of non-genetic resistance acquisition in cancer cell populations has been recently published (Chisholm et al. ). It considers the emergence of drug tolerance in certain type of cancer cell populations treated by a chemotherapeutic agent, which has been experimentally linked to reversible epigenetic changes (Sharma et al. ). The model accounts for two phenomenological, stochastically fluctuating cell traits: their survival potential (i.e., drug resistance) and their proliferation potential. Therefore, under drug exposure, cells with low survival potential are eliminated while other can survive. In addition to this selection effect, the model also accounts for cell-level, stress-induced adaptation of the proliferation level: in presence of drug, cells decrease their proliferation. Hence, cell population dynamics is here governed by an interplay between phenotypic fluctuations, drug-mediated selection on those fluctuations, and druginduced cellular adaptation. This situation is very similar to our simulations of repeated TRAIL treatments. However, a main difference is that our model is 'mechanistic' in the sense that cell response to TRAIL is governed by existing biochemical reactions between proteins. Hence, cell-level and selection-based resistance acquisition is an emerging property of the model and does not result from modeling choices. 'Mechanistic' cell population dynamics models with molecular details are rare but might become more common in the next years. In a recent study (Shokhirev et al. ), a multi-scale model of how N F -κB, cell cycle and apoptosis molecular pathways drives the population dynamics of lymphocytes was constructed and compared to experimental data. Because they relied on ODE models, they modeled protein level variability with static distributions. To solve the inherent problem with that (as already discussed, for example overestimation of selection effects), they decided to remix this 'extrinsic' noise by resampling protein levels from the corresponding distributions each time a cell divides. Our approach provides a less artificial and more mathematically-grounded solution to this problem, and hence is less likely to introduce biases and allows to use temporal decorrelation data to finely calibrate the model. From stochastic gene expression in cells to spatially-organized populations While dissecting the quantitative functioning of cellular pathways in single cells requires precise control of the environment (Spiller et al. ), cells in a population are usually not sensing an homogeneous environment. Moreover, they often themselves create environmental heterogeneity by the way they grow and arrange in space or by adopting different phenotypes (Snijder et al. ; Hirata et al. ). For example, this is the case for tumor spheroids treated by TRAIL (as we studied in Chapter ) or for cell populations able to send and sense messenger molecules (as for micro-colonies of engineered yeast cells studied in Chapter ).In this thesis, we used spatial cell-based modeling to investigate those situations. Spatial cell-based modeling has been increasingly used during the last years to investigate the emerging behavior that results from environment-cell and cell-cell interactions (Graner & Glazier ; Drasdo et al. ; Galle et al. ; Walker et al. ; Leeuwen et al. ; Ramis-Conde & Drasdo ; Jang et al. ; Macklin et al. ; Schlüter et al. ). Most used relatively abstract rules for cellular decisions, while some modeled molecular species (Walker et al. ; Ramis-Conde & Drasdo ; Schlüter et al. ). However, spatial cell-based model that include stochastic gene expression in biochemical pathways are rare (to the best of our knowledge, this work represents the first a empt).

  5 pC8 R * + k 6 pC8 : R * k 14 pC8 C6 + k 15 pC8 : 8 pC3 C8 + k 9 pC3 : C8 k 69 pC3 Apop + k 70 pC3 : 11 pC6 C3 + k 12 pC6 : C3 dP ARP dt = k P ARP sp mRN A P ARPγ P ARP p P ARPk 17 P ARP C3 + k 18 P ARP : C3 df lip dt = k f lip sp mRN A f lipγ f lip p f lipk 20 f lip R * + k 21 f lip : R * dBar dt = k Bar sp mRN A Barγ Bar p Bark 22 Bar C8 + k 23 Bar : C8 dXIAP dt = k XIAP sp mRN A XIAPγ XIAP p XIAPk 24 C3 XIAP + k 25 C3 : XIAP + k 26 C3 : XIAPk 50 XIAP Apop + k 51 XIAP : Apop k 62 smac XIAP + k 63 smac : XIAP dBid dt = k Bid sp mRN A Bidγ Bid p Bidk 27 Bid C8 + k 28 Bid : C8 dBax dt = k Bax sp mRN A Baxγ Bax p Baxk 30 Bax tBid + k 31 Bax : tBid dP ore dt = k P ore sp mRN A P oreγ P ore p P orek 39 P ore Bax * m,4 + k 40 P ore : Bax * m,4 dsmac m dt = k smac m sp mRN A smac mγ smac m p smac mk 55 smac m P ore * + k 56 smac m : P ore * k 52 CyC m P ore * + k 53 CyC m : P ore * k 42 M cl1 tBid + k 43 M cl1 : tBiddBcl2 dt = k Bcl2 sp mRN A Bcl2γ Bcl2 p Bcl2k 44 Bcl2 Bax * m + k 45 Bcl2 : Bax * m k 46 Bcl2 Bax * m,2 + k 47 Bcl2 : Bax * m,2k 48 Bcl2 Bax * mk 1 T RAIL R + k 2 T RAIL : R + k 4 R *dT RAIL : R dt = -γ T RAIL:R p T RAIL : R + k 1 T RAIL Rk 2 T RAIL : R k 3 T RAIL : R dR * dt = -γ R * p R * + k 3 T RAIL : Rk 4 R * k 5 pC8 R * + k 6 pC8 : R * + k 7 pC8 : R *k 20 f lip R * + k 21 f lip : R * dpC8 : R * dt = -γ pC8:R * p pC8 : R * + k 5 pC8 R *k 6 pC8 : R * k 7 pC8 : R * dC8 dt = -γ C8 p C8 + k 7 pC8 : R *k 8 pC3 C8 + k 9 pC3 : C8 + k 10 pC3 : C8 + k 16 pC8 : C6k 22 Bar C8 + k 23 Bar : C8k 27 Bid C8 + k 28 Bid : C8 + k 29 Bid : C8 dpC3 : C8 dt = -γ pC3:C8 p pC3 : C8 + k 8 pC3 C8k 9 pC3 : C8 k 10 pC3 : C8 dC3 dt = -γ C3 p C3 + k 10 pC3 : C8k 11 pC6 C3 + k 12 pC6 : C3 + k 13 pC6 : C3 k 17 P ARP C3 + k 18 P ARP : C3 + k 19 P ARP : C3k 24 C3 XIAP + k 25 C3 : XIAP + k 71 pC3 : Apop dpC6 : C3 dt = -γ pC6:C3 p pC6 : C3 + k 11 pC6 C3k 12 pC6 : C3 k 13 pC6 : C3 dC6 dt = -γ C6 p C6 + k 13 pC6 : C3k 14 pC8 C6

+ k 67 +:

 67 Apaf * pC9k 68 Apop k 69 pC3 Apop + k 70 pC3 : Apop + k 71 pC3 : Apop dXIAP : Apop dt = -γ XIAP :Apop p XIAP : Apop + k 50 XIAP Apopk 51 XIAP : Apop dCyC m : P ore * dt = -γ CyC m :P ore * p CyC m : P ore * + k 52 CyC m P ore *k 53 CyC m : P ore * k 54 CyC m : P ore * dCyC : r dt = -γ CyC:r p CyC : r + k 54 CyC m : P ore *k 58 CyC : r + k 59 CyC dsmac m : P ore * dt = -γ smac m :P ore * p smac m : P ore * + k 55 smac m P ore *k 56 smac m : P ore * k 57 smac m : P ore * dsmac : r dt = -γ smac:r p smac : r + k 57 smac m : P ore *k 60 smac : r + k 61 smac dCyC dt = -γ CyC p CyC + k 58 CyC : rk 59 CyC k 64 Apaf CyC + k 65 Apaf : CyC + k 66 Apaf : CyC dsmac dt = -γ smac p smac + k 60 smac : rk 61 smac k 62 smac XIAP + k 63 smac : XIAP dsmac : XIAP dt = -γ smac:XIAP p smac : XIAP + k 62 smac XIAPk 63 smac : XIAP dApaf : CyC dt = -γ Apaf :CyC p Apaf : CyC + k 64 Apaf CyCk 65 Apaf : CyC k 66 Apaf : CyC dApaf * dt = -γ Apaf * p Apaf * + k 66 Apaf : CyCk 67 Apaf * pC9 Apop + k 69 pC3 Apopk 70 pC3 : Apop k 71 pC3 : Apop

cell survive TRAIL exposure: molecular events counteracting apoptotic signaling

	What could be the molecular mechanisms that promote cell survival in absence of genetic
	defects in the TRAIL apoptosis machinery? This question was at the center of extensive
	research efforts. Many mechanisms of various nature were reported. Some are well under-
	stood but many are poorly characterized (Falschlehner et al.	).

killing is obtained (Song et al.

;

Flusberg et al. )

. This reversible resistance (or transient resistance acquisition as we prefer to call it) could also have important implications for therapy, notably for the scheduling of drug applications. Together, those observations demonstrate the existence of non-genetic factors in TRAIL-induced apoptosis that cause ) cell fate variability among naïve (i.e. previously unexposed to TRAIL) cell populations and

) transient resistance acquisition for the surviving population after TRAIL treatment.

How can a

Those mechanisms can be classified into three categories. First, competitive binding of antiapoptotic proteins can prevent activated death signaling proteins to activate downstream targets. This is the case for example for decoy receptors towards TRAIL (Sheridan et al.

), for FLIP towards DISC (Shirley & Micheau ), for XIAP towards caspase-, for Mcl-towards Bid (Youle & Strasser ) and for Bcl towards Bax (Oltvai et al.

Description, characterization and simulation of the stochastic protein turnover model. It is

	Stochastic protein turnover model		protein			protein
	Protein	γ p	Numerical		mRNA			mRNA
	k sp		simulations		promoter			promoter
		γ m			Fluctuations in			Steady-state population
		mRNA			single cells			variability
	k sm			protein	Autocorrelation Prot = f (k on ,k o ,k sm , m ,k sp , p )	protein	Mean Prot = f (k on ,k o ,k sm , m ,k sp , p ) Var Prot = f (k on ,k o ,k sm , m ,k sp , p )
		k on	Analytical results	mRNA	Autocorrelation mRNA = f (k on ,k o ,k sm , m )		mRNA	Mean mRNA = f (k on ,k o ,k sm , m ) Var mRNA = f (k on ,k o ,k sm , m )
	"ON"	promoter k o	"OFF"	promoter	Autocorrelation Gene = f (k on ,k o )		promoter	Var Gene = f (k on ,k o ) Mean Gene = f (k on ,k o )
	Figure 2.2: defined by six rates governing model reactions. Gene activity switches, mRNA production and degradation
	(red arrows) are stochastic reactions. Protein synthesis and degradation reactions (black arrows) are deter-
	1 ministic. Numerical simulations can be used to generate the corresponding temporal fluctuations in single . cells. When a population of cells is simulated, the cell-to-cell variability can be studied. After some time, k sm This popular model (Raj et al. ; Paszek ; Shahrezaei & Swain cell-to-cell variability reaches a steady state. Analytical calculations provide expressions characterizing ; Suter et al. ; Dar et al. ; Singh et al. ; Viñuelas et al. ; Sherman & Cohen the steady-state variability (moments of the steady-state distributions), but also fluctuations (autocorrelation ) has functions).
	many names, such as the central dogma model of gene expression or the two-state transcriptional
	bursting model. Here we will call it protein fluctuation model or alternatively stochastic protein
	turnover model to highlight the importance of the protein half-life in shaping fluctuations (as
	will be shown later).
	Of note, the variant we consider does not account for noise in translation and protein
	degradation. This choice could be justified by the high copy number in which proteins are
	generally found in human cells (Schwanhäusser et al.	). It facilitates the numerical
	simulation of the model without complicating much its analytical treatment. The validity
	of this approximation is discussed in more details elsewhere (Paszek	).

  ). It informs about the typical amplitude of the protein fluctuations and depends on the rate parameters k on , k off , k sm , γ m and γ p . Note that it does not depend on k sp : because protein translation and degradation are deterministic, k sp simply acts as a scaling factor that sets the mean protein level.

	The protein level steady-state distribution does not contain information about the speed
	at which it is visited in individual cells through fluctuations. A simple argument to this
	claim is that scaling all rate parameters by the same factor would leave the steady-state
	distribution unchanged as it amounts to rescale time.
	A natural way to characterize this speed of fluctuations is to rely on the auto-correlation
	function of the protein level. This function measures how much information (quantified by
	the correlation) about the initial protein level remains after a given time. In other words, it
	gives the time evolution of the memory of single cells with respect to their protein level. Here,
	we provide an exact derivation of the auto-correlation function for the stochastic protein
	turnover model. This derivation and the resulting expression are presented in Appendix .
	Three different autocorrelation curves corresponding to different model parameteriza-
	tions are shown in Figure . . To characterize this function with a single number, one can
	extract the time needed to lose half of the correlation with initial protein level. This value
	is called the mixing time and has already been measured experimentally for endogeneous
	proteins in human cells (Sigal et al.	).
	In summary, we propose here to use the two real-valued statistics CV (of the protein
	level distribution) and mixing time τ (of the protein level) to characterize the amplitude and
	speed of protein fluctuations generated by a given stochastic protein turnover model. A more
	detailed discussion on the parameterization of this model from constraints on CV and τ is
	given in Chapter (see . , section Parameterization of stochastic protein turnover models from
	protein variability and mixing time).

Native proteins in the kinetic model (17) stable proteins: equip with a standard stochastic protein turnover model short-lived proteins: Flip and Mcl1

		γ p	such that half-life = 27 hrs (dilution)	γ p	such that half-life ranges between 0.1 and 1.5 hours
		γ m	such that half-life = 9 hrs (typical mRNA stability) (Schwanhausser et al., 2011)	γ m	such that half-life ranges between 1 and 4 hours
	standard constraints	k on + k Mean protein level Mean mRNA	= 10.4 hrs -1 (typical switching speed) (Suter et al., 2011) (Schwanhausser et al., 2011) (typical value) = 17	k on level Mean mRNA k	Wide range explored (Schwanhausser et al., 2011) (typical value) = 17 Wide range explored	+
	particular constraints	level Variance of	as in (Spencer et al., 2009)
		protein level		Mean protein
				level
		All rates determined
	lives (Vucic et al.	). To account for this fact, we used a unique half-life ( hours) for
	all non-native species, except for the mitochondrial pores ( . hours, same value as in the
	original model) and the complexes involving the short-lived protein FLIP and Mcl-(the
	half-life of the native form was used, based on the assumption that FLIP and Mcl-are
	targeted for degradation in a similar fashion alone or with their binding partner).

All rates determined particular constraints + explore reasonable ranges one standard constraint to keep space to explore tractable + as in

  While we cannot exclude that model parameterizations allowing a quantitative agreement exist, it might be needed to include additional mechanisms such as survival pathways induction to explain the observed sustained resistance gain after one week when treating cells with a high TRAIL dose(Flusberg et al.

		al.	; Flusberg &
	Sorger	) is only qualitative.

). It should also be noted that the two experimental measurements of resistance gain in one week survivors were obtained with different methods. It would thus be interesting to confirm the sustained resistance observation for high TRAIL dose in the video-microscopy experimental se ing, which allows a more direct and precise measurement of changes in the overall resistance of populations (provided that a sufficient number of cells are analyzed).

Impact on short-term model behavior of the presence/absence of caspase feedback loop and the degradation of non-native ('active') forms.

  Model-Data agreement is shown for MOMP time distributions, surviving fractions and sisters correlation of MOMP time in both treatment conditions for model variants when the caspase feedback loop is either present/absent and the default active forms half-life is 15, 5 or 2 hours.

	With feedback loop	Without feedback loop
	Default active forms half-life = 15 hours Default active forms half-life = 5 hours Default active forms half-life = 2 hours Figure 2.18: 2 0 5 0 0.2 0.4 0.6 0.8 1 0 2 4 6 0 0.05 0.1 0.15 0.2 Data , -CHX Model , -CHX 8 Data , +CHX Model , +CHX 0 5 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 0 0.05 0.1 0.15 0.2 Data , -CHX Model , -CHX Data , +CHX Model , +CHX 0 5 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 0 0.05 0.1 0.15 0.2 Data , -CHX Model , -CHX Data , +CHX Model , +CHX With feedback loop 10 15 Data , -CHX Data , +CHX Model , -CHX 20 Model , +CHX -CHX 0 0.2 0.4 0.6 0.8 1 10 15 20 Data , -CHX Data , +CHX Model , -CHX Model , +CHX -CHX 0 0.2 0.4 0.6 0.8 1 10 15 20 Data , -CHX Data , +CHX Model , -CHX Model , +CHX -CHX 0 0.2 0.4 0.6 0.8 1 With feedback loop Non-native forms default half-life +CHX Data Model 0 2 4 0 0.05 0.1 0.15 0.2 0 0 0.2 0.4 0.6 0.8 1 +CHX Data Model 0 2 4 0 0.05 0.1 0.15 0.2 0 0 0.2 0.4 0.6 0.8 1 +CHX Data Model 0 2 4 0 0.05 0.1 0.15 0.2 0 0 0.2 0.4 0.6 0.8 1 Without feedback loop 6 8 Data , -CHX Model , -CHX Data , +CHX Model , +CHX -CHX 0 0.2 0.4 0.6 0.8 1 5 10 15 Data , -CHX Data , +CHX Model , -CHX Model , +CHX 6 8 Data , -CHX Model , -CHX Data , +CHX Model , +CHX -CHX 0 0.2 0.4 0.6 0.8 1 5 10 15 Data , -CHX Data , +CHX Model , -CHX Model , +CHX 6 8 Data , -CHX Model , -CHX Data , +CHX Model , +CHX -CHX 0 0.2 0.4 0.6 0.8 1 5 10 15 Data , -CHX Data , +CHX Model , -CHX Model , +CHX Without feedback loop Default active forms half-life = 15 hours Default active forms half-life = 5 hours Default active forms half-life = 2 hours Default active forms half-life 5 10 15 27 0 0.5 1 1.5 2 2.5 3 Experimental Measurement No Feedback Loop With Feedback Loop Resistance gain in one-day survivors Resistance gain in one-day survivors	+CHX +CHX +CHX	Data Model Data Model Data Model

  sistance in such se ings and how it relates to the in-vitro resistance of the same cell lines remains limited. Indeed, those studies mostly focused on demonstrating qualitatively the anti-tumor activity of a new DRA or a new combination of a DRA with another agent.

	Ashkenazi et al.	; Chinnaiyan et al.	; Fulda et al.	; Merchant et al.	;
	DeRosier et al.	; Lee et al.	; Huet et al.	), our understanding of TRAIL re-
	Mathematical modeling and numerical simulations are potentially very helpful to build
	such quantitative understanding because it allows to translate previous knowledge and hy-
	potheses on the mechanisms at play into predictions that can be compared with experimen-
	tal observations. However, no mathematical models enabling the investigation of TRAIL
	resistance over long time scales and for cell populations exhibiting spatially heterogenous
	micro-environments have yet been developped. A potential explanation for this is that so
	far no model able to quantitatively account for the observed dynamics of cell-to-cell variabil-

cacy. However, despite numerous in-vivo studies on DRA treatments

(Walczak et al. ;

ity was available. Moreover, studying spatial and long-term effects requires the integration of quantitative intra-cellular models into quantitative cell-based models, which represents a challenge.

In this Chapter, we integrate our single-cell model of TRAIL-induced apoptosis into multi-cellular simulations to investigate TRAIL resistance in two virtual experimental settings that could bridge the gap between standard in-vitro and in-vivo se ings. The first se ing corresponds to the long-term in-vitro, monolayer culture (several weeks) of cell populations repeatedly treated by TRAIL. The second se ing extends the first but instead cells are grown as D multi-cellular spheroids, allowing us to investigate the additional role of spatial heterogeneity in TRAIL exposure. Although spheroids are smaller than in-vitro tumors, their size is higher or comparable to the typical distance between two vessels within in-vivo tumors (around -micrometers, (Primeau et al. ; Baker et al. )), making it relevant to study heterogeneity in drug spatial distribution. Those two se ings are intermediate between standard assays of TRAIL cytotoxicity and pre-clinical in-vivo assays in mice. Importantly, they are experimentally relatively easy to perform, and much more amenable than in-vivo experiments. Therefore, our modeling predictions are readily testable experimentally.

Table 3 .

 3 1: List of parameter values for computing the motion of cells.

	Parameter	Symbol	Value	Unit
	Young modulus	E		P a
	Poisson ratio	ν	.	none
	Adhesivity	σ	10 -4	J.m -2
	Friction coefficient	µ	10 7	P a.s.m -1
	Micro-motility persistence	τ M M		min
	time
	Micro-motility average	F M M avg	10 -10	N
	force
	Micro-motility force	CV F M M	.	none
	variability
	Maximum displacement	∆l max	.	µm
	per motion step

Table 3 .

 3 2: List of parameter values for simulating cell growth and division.

	Parameter	Symbol	Value	Unit
	Average cell cycle	CC avg		hours
	duration
	Cell cycle duration	CV CC	.	none
	variability
	Cell radius at birth	R birth	.	µm
	Cell radius at division	R div	2 1/3 R birth	µm
	Daughter cells distance	d div	2 R birth	µm
	at division

Table 3 .

 3 3: List of parameter values for simulating TRAIL penetration in the spheroid.

	Parameter	Symbol	Value	Unit
	TRAIL diffusion in	D medium	50	µm 2 .s -1
	absence of cells
	TRAIL diffusion in	D cells	0.2	µm 2 .s -1
	presence of cells
	TRAIL half-life in	HL medium T RAIL		hours
	absence of cells

Table 3

 3

	Parameter	Symbol	Value	Unit
	Radius reduction rate	γ lysis	varied ( or )	µm.hr -1
	of dying cells
	Radius under which
	cells are removed

.4: List of parameter values for simulating TRAIL-induced apoptosis and lysis of dying cells.

Table 4 .

 4 1: Reactions of the IP sensing semi-mechanistic ODE model. Note that because IP is supposed to be in large excess in medium, IP binding does not reduce extra-cellular IP concentration. (*): such that steady-state level (in absence of IP) is 1. (**): such that GFP steady-state level (in absence of IP) is 1.

	Reaction	Description	Rate
	IP + R → IP : R + IP	IP binding	kb 1
	IP : R → R	IP unbinding	ku 1
	IP : R + T F → IP : R : T F TF binding	kb 2
	IP : R : T F → IP : R + T F TF unbinding	ku 2
	IP : R : T F → IP : R+T F * TF activation	kc 1
	T F * → T F * + uGF P	TF-dependent
		uGFP synthesis

  Computational efficiency of the mRNA trajectories storage algorithm as a function of the storage time for two ODE solver methods: Dopr853 and SIE. The absolute and relative error tolerances were 10 -3 in both cases. For each method, comparison to the performance of the naïve algorithm is provided. In each case, 100 cells have been simulated for 8 hours after a treatment of 250 ng/mL TRAIL. A loglog scale is used.

					mRNA storage, DOPR853 method
					naïve algorithm, DOPR853 method
		10 3			mRNA storage, SIE method naïve algorithm, SIE method
	computation time per cell (ms)	10 2
		10 -1	10 0	10 1	10 2
	Figure 5.2:

  R * + k 20 f lip R *k 21 f lip : R * + k 41 P ore : Bax * m,4k 52 CyC m P ore * + k 53 CyC m : P ore * + k 54 CyC m : P ore * k 55 smac m P ore * + k 56 smac m : P ore * + k 57 smac m : P ore *

	dcP ARP dt df lip : R * dt dBax * m,4 = -γ cP ARP p = -γ dt = -γ Bax * cP ARP + k 19 P ARP : C3 m,4 p Bax * m,4 + k 37 Bax * m,2 Bax * m,2 -k 38 Bax * m,4 f lip:R * p -k 39 P ore Bax * m,4 + k 40 P ore : Bax * m,4 f lip : dBar : C8 dt = -γ Bar:C8 p -k 48 Bcl2 Bax * m,4 + k 49 Bcl2 : Bax * m,4 Bar : C8 + k 22 Bar C8 -k 23 Bar : C8 dC3 : XIAP dt = -γ C3:XIAP p C3 : XIAP + k 24 C3 XIAP -k 25 C3 : XIAP dP ore : Bax * m,4 dt = -γ P ore:Bax * m,4 p P ore : Bax * m,4 + k 39 P ore Bax * m,4 -k 40 P ore : Bax * m,4
		-k 26 C3 : XIAP -k 41 P ore : Bax * m,4
	dC3 : ub dt dBid : C8 dt dtBid dt dP ore * = -γ C3:ub p = -γ Bid:C8 C3 : ub + k 26 C3 : XIAP p Bid : C8 + k 27 Bid C8 -k 28 Bid : C8 -k 29 Bid : C8 = -γ tBid dt = -γ P ore * p P ore * dM cl1 : tBid = -γ M cl1:tBid p M cl1 : tBid + k 42 M cl1 tBid -k 43 M cl1 : tBid dt dBcl2 : Bax * m = -γ Bcl2:Bax * m p Bcl2 : Bax * m + k 44 Bcl2 Bax * m -k 45 Bcl2 : Bax * m dt dBcl2 : Bax * m,2 = -γ Bcl2:Bax * m,2 p Bcl2 : Bax * m,2 + k 46 Bcl2 Bax * m,2 -k 47 Bcl2 : Bax * m,2 dt dBcl2 : Bax * m,4 dt = -γ Bcl2:Bax * m,4 p Bcl2 : Bax * m,4 + k 48 Bcl2 Bax * m,4 -k 49 Bcl2 : Bax * m,4
	+ k 34 Bax * m dApop dt = -γ Apop
	dBax * m dt	= -γ	Bax * m p	Bax * m + k 33 Bax * -k 34 Bax * m
		-2k 35 Bax * m Bax * m + 2k 36 Bax * m,2
		-k 44 Bcl2 Bax * m + k 45 Bcl2 : Bax * m
	dBax * m,2 dt	= -γ	Bax * m,2 p	Bax * m,2 + k 35 Bax * m Bax * m -k 36 Bax * m,2
		-2k 37 Bax * m,2 Bax * m,2 + 2k 38 Bax * m,4
		-k 46 Bcl2 Bax * m,2 + k 47 Bcl2 : Bax * m,2
		+ k 15 pC8 : C6 + k 16 pC8 : C6
	dpC8 : C6 dt	= -γ	pC8:C6 p

pC8 : C6 + k 14 pC8 C6k 15 pC8 : C6 k 16 pC8 : C6 dP ARP : C3 dt = -γ P ARP :C3 p P ARP : C3 + k 17 P ARP C3k 18 P ARP : C3 k 19 P ARP : C3 p tBid + k 29 Bid : C8k 30 Bax tBid + k 31 Bax : tBid + k 32 Bax : tBid k 42 M cl1 tBid + k 43 M cl1 : tBid dBax : tBid dt = -γ Bax:tBid p Bax : tBid + k 30 Bax tBidk 31 Bax : tBid k 32 Bax : tBid dBax * dt = -γ Bax * p Bax * + k 32 Bax : tBidk 33 Bax * p Apopk 50 XIAP Apop + k 51 XIAP : Apop

Acknowledgements

Introduction

Acknowledgements

DATA 10 1

Estimation of stochastic gene expression parameters using single-cell level induction dynamics

We then asked whether the single-cell IP induction dynamics (Figure . ), which revealed an important cell-to-cell variability, can be captured with our single-cell extension of the IP sensing model (Figure . ). As for modeling protein fluctuations in TRAIL-induced apoptosis (Chapter ), we followed a parsimonious parameterization approach. First, none of the parameters estimated previously from the population averaged data and that still appear in the single-cell model were changed. Moreover, we used the same transcriptional bursting parameters for the three proteins considered in the model (R, TF and GFP) and started with reference values obtained from a previous study of transcriptional bursting in budding yeast (Brown et al. ).

In a first step, we manually tuned those parameters in order to get a distribution of GFP in absence of IP that resemble observations. Then, we also manually adjusted the three parameters of TF* -mediated transcription of GFP by comparing the simulated induction dynamics single-cell distributions with observations. Strikingly, we very rapidly obtained a good semi-quantitative agreement (Figure . ). The model correctly predicts detailed features of the induction distributions present in the data, such as the straight, non-smooth angle at the mode of the distribution observed for inducing with IP = 0.12 µM , and the bimodal distributions for doses ranging from 0.18 µM to 0.41 µM . Thus, our single-cell model of IP sensing is able to reproduce observed single-cell IP induction dynamics with a good accuracy. Note that although at this stage we were satisfied with this result, an even be er agreement could probably be obtained by using an automated fi ing procedure.

Analysis of solid culture data

In order to assess if our model-based characterization of IP sensing obtained from liquid culture data can predict the observed behavior of receiver cell micro-colonies in presence of IP gradients (Figure . ), we first had to quantify this response in way that allows meaningful comparison with model predictions. To this aim, we developped an image analysis pipeline to compute a fluorescence temporal profile for each micro-colony (Figure . ). Of note, because piling up occured at the colony centers and inter-colony growth was variable, we segmented the outer rim of each colony for which a single layer of cells contribute to the measured fluorescence. We show the resulting profiles for three representative microcolonies (Figure . ). As could be expected if IP diffusion is sufficiently slow, clear differences of fluorescence profiles between the colonies are obtained, consistently with their distance to the location at which IP was deposited: close colonies responded faster and to a higher extent than distant colonies.

Estimation of the IP diffusion coefficient

Under the assumption that the physiology of receiver cells does not change much between liquid and solid culture conditions (at least with respect to the functioning of the IP sensing module), an in-silico reproduction of this experiment, using our model of IP sensing calibrated on liquid culture data and in which only the IP diffusion coefficient parameter is varied, should be able to predict this data. Indeed, the geometry of the gel is known, the stability of IP is probably not much affected, and the exact amount of IP deposited is known.

We tested a few values of the IP diffusion coefficient and found a good agreement between model predictions and data (Figure . ) for D IP = 10 -4 mm 2 .s -1 . This value is in the range of what could be expected from previous knowledge (see for example (Brown & Johnsen )). Thus, our model-based characterization of IP sensing derived from liquid culture data is quantitatively consistent with the behavior of receiver cell micro-colonies in

Simulation algorithms

To simulate such hybrid stochastic/deterministic models, a straightforward, 'naïve' simulation algorithm can be proposed:

Data: vector of cell state (gene activities and mRNA levels as integers, protein levels as floats) at time Note that this algorithm leverages the fact that in our case, the rate of stochastic reactions depends only on the discrete variables (gene activity status and mRNA levels), allowing the use of the Gillespie algorithm for the stochastic reactions. More generally, for a given partition of a reaction scheme into stochastic and deterministic reactions (the species changed by the stochastic reactions are then the discrete variables, and the rest the continuous variables), this algorithm is valid only if the rates of stochastic reactions do not depend on continuous variables. Otherwise stated, the evolution of stochastic variables is independent from the evolution of deterministic variables (but the inverse dependence is possible). This algorithm calls the ODE solver between each stochastic event, which might be costly if the ODE solver has some constant component in its execution time. Thus, we propose a potentially improved algorithm, that utilizes the independence of stochastic variables evolution: we simulate many steps of gene activity switches, mRNA production and mRNA degradation and we store the resulting mRNA trajectories in a look-up table before calling an ODE solver, which will read mRNA levels in that table as time-dependent parameters of the differential equations, as described in Figure . .

Model equations

Here we give all the differential equations governing the evolution of protein levels. The mRNA levels are governed by stochastic reactions as described earlier.

Appendix : Simulation of the single-cell model of TRAIL-induced apoptosis and comparison with experimental data

Simulating TRAIL-induced apoptosis with stochastic protein turnover

In general, the TRAIL signaling protein-protein reactions are taking place concurrently with stochastic protein turnover. When the noise in signaling reactions is neglected due to high protein copy number, those reactions can be simulated using ODEs. However, rates of protein synthesis are in our model stochastic, as they follow mRNA fluctuations. Promoter activity and mRNA fluctuations were simulated using an implementation of the Gillespie algorithm in C++ (Numerical Recipes). Messenger RNA trajectories were computed and stored in advance because protein levels do not affect the rates of promoter state switches, mRNA production and degradation (see Chapter ). The ODEs governing evolution of all protein levels were then simulated using the Semi-Implicit Extrapolation method implemented in C++ ((Deuflhard ), Numerical Recipes). This method was significantly faster than a more standard Runge-Ku a method (Dormand-Prince, C++, Numerical Recipes) but gave identical results (see Chapter , section Simulating combined signal transduction and stochastic gene expression).

In-silico sister cells experiment

To sample the state (promoter activity, mRNA and protein levels) of the mother cells, all stochastic protein turnover models were simulated during days (starting for the mean state) for each of the ⁴ (or ⁵ for Figures . , . b) and -g) in this reference), we applied a sampling algorithm to approximately reproduce those distributions. The overall impact on correlation curves was generally low compared to results obtained with assuming a fully uniform distribution of division times in the pre-stimulus recording interval. MOMP was considered to have occurred when half of mitochondrial Smac has been released.

In-silico repeated TRAIL experiment

A naïve population of ⁴ cells was obtained as in the sister cells experiment. Each cell was assigned a random time of next division. To account for the fact that the distribution of next division times is not uniform in growing cell populations, we used a distribution obtained by simulating simple growth. New cells were a ributed a next division time according to a cell cycle duration normally distributed with hours mean and hours standard deviation. Cells in which cPARP levels exceeded ⁵ were considered dead as in (Gaudet et al. ). To closely mimic the experimental protocol used in (Flusberg et al.

), we accounted for the effect of passing cells by checking population size each day and if needed, removing randomly cells until ⁴ were left. Resistance gain is computed as RG = Resistance survivors Resistance naive cells where Resistance = AliveCellN umber (treatment+8hours) AliveCellN umber (treatment) ,

similarly to (Flusberg et al. ).

Quantification of model-data agreement

Comparison between the different types of data

To permit comparison between the different types of data, in each case a threshold for the cost defining agreement/disagreement was manually set by visual comparison of experi- Appendix : List of software tools used in this thesis