Acknowledgment

First and foremost, I want to thank God for blessing me with this life and for being there in good and bad times.

I would like to thank my supervisor, Dr. Carito Guziolowksi, who I met for the first time when I arrived to Nantes by train at 11 o clock at night. She was there to pick me up at this late hour. She won my heart right there by this simple act of kindness. I knew at that moment that I am going to learn science as well as humanity from her. She has played a significant role in my life in these 3 years of PhD. Her constant support and guidance has helped me to understand my thesis subject properly. I would also like to thank my supervisor, Dr. Jeremie Bourdon, who I get to know in the mid of my PhD. He has provided me with immense support and guidance during my thesis work. I have also learned a great deal of politeness from him. I am also grateful to my parents, for not following the typical traditions and helping me to get the education, which allowed me to write this thesis. My mother, Rubeena, encouraged me to be an independent woman. My father, Razzaq, appreciated my individuality instead of making me follow the crowd. My brothers, Adnan and Munsfeen, provided constant loving support even during the times when we disagreed. I would also like to thank my special friend, Roland, with whom I enjoyed an 80km cycling trip and 1000m hike (record of my life, almost died). He also annoyed me a lot with his ambition of teaching me how to use articles. He taught me the importance of persistence and confidence.

I am thankful to my friends, Julia, Abhilash, and Khaoula for taking care of me when I got sick which happened quite a lot during first two years of my PhD. Thanks to Julien for taking care of my plants and listening to me during stressful times. Thanks to Simon for trying to teach me how to swim which I still cannot do. I am also grateful to Bodgan, Konstatin, Sylvain, Deep, Tahir, Sebestien and (annoying) Saman for great moments and discussions.

I would like to end this section with this beautiful quote: "Your task is not to seek for love, but merely to seek and find all the barriers within yourself that you have built against it." Rumi.

Introduction

Systems biology is the study of biological functions and mechanisms consisting of DNA, RNA, and proteins in a systematic way. It is also concerned with understanding the complex and dynamic biological organization within the living organism in order to explain experimentally observed behaviors and do future predictions [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF]. Cells are the fundamental units of a living organism and share the same building blocks consisting of DNA, RNA and proteins with different complexity and structure. The central dogma of molecular biology explains the continuous reproduction of cells and flow of genetic information. It states that the blueprint of the DNA is carried by the mRNA (messenger RNA) to make a "protein" as a functional product, which performs the cellular work. Figure 1.1 explains the processes involved in the central dogma of biology: (1) Replication, (2) Transcription, and (3) Translation. Replication is the process of generating new DNA as functional product. Transcription produces new RNA as a result. Translation results in a protein as final product. Even after translation, proteins have to go through many modifications in order to become fully functional [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF].

These modifications involve the folding into three dimensional structure and moving to a specific location in the cell. Each cell contains different types of proteins and each protein can be expressed differently under different circumstances. Proteins play an integral role in all cellular functions of living organisms, i.e., creating a signaling cascade and regulating various biological mechanisms. Proteins are made up of the primary sequence of 20 amino acids. Usually proteins interact with other proteins to create complexes of varying sizes. These complexes work with other proteins or complexes to create pathways or modules to carry out tasks in a cell such as signal transduction, metabolism, duplication, DNA transcription, and DNA damage repair, etc. Understanding the interactions among proteins and complexes is crucial to fully identify and characterize the structure and functions of the cell machinery.

There have been a lot of advances in experimental, analytical and computational techniques in the last couple of decades. These advances helped to generate large amounts of proteomics data, such as the yeast two-hybrid system and mass spectrometry. Large scale protein to protein interaction databases contain a large number of experimentally verified protein to protein interactions. These advances have made things less clear and at the same time they helped to understand the protein functions, interactions, and structure at varying levels. For example, linear dynamics were used to model signal transduction inside cells; now it seems that signals follow the more complex nonlinear dynamics. Understanding the mechanisms and interactions among proteins helps to unravel how information propagates within cells in different diseases such as cancer [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF][START_REF] Duan | The roles of post-translational modifications in the context of protein interaction networks[END_REF]. In this thesis, I am focused on inferring and understanding the protein signaling networks in four breast cancer cell lines (BT20, BT549, MCF7, UACC812).

Motivation

Protein signaling networks are not static in nature since they respond to stimuli and perturbation. They constitute complex regulatory systems controlled by crosstalk and feedback mechanisms. These networks are often regulated in diseases. Discovering the precise mecha-1.1. MOTIVATION nisms of signal transduction may provide a better fundamental understanding of disease behavior. For instance, a main difficulty in cancer treatment is that different signaling networks reveal that cell populations specialize upon treatment and therefore patient responses may be heterogeneous. Computational models of signaling control for different patient groups could guide cancer research towards a better drug targeting system. In this work, we propose a methodological framework to discriminate among the regulatory mechanisms of four breast cancer cell lines by building predictive computational models.

Boolean network (BN) modeling is a simple yet powerful framework to study biological models such as signaling pathways or regulatory networks. BNs are based on qualitative approaches, allowing to model large scale biological networks [OPS + 15]. Part of the research about this paradigm focuses on the topology of the network, searching for interesting characteristics such as cycles and hubs. Some researchers focus on identifying the influences among the components of the network, knowing which gene or protein activates or inhibits the others. There is also research about studying the dynamics of the model, how the combinations of inherent influences make the system evolve and how this evolution changes under different conditions or situations. While many BN approaches exist to model biological systems, they focus mainly on system properties, and few exist to integrate experimental data in them. In this thesis, we use the caspo time series (caspo-ts) method to learn cell line specific BNs by integrating protein signaling networks with experimental data. The caspo-ts method uses Answer Set Programming (ASP) and Model Checking techniques to solve the combinatorial optimization problem of enumerating a family of Boolean networks (BNs) optimally explaining time-series data [OPS + 16]. Figure 1.2 shows the overall process of caspo-ts, a publicly available software at [START_REF] Razzaq | Caspo-ts[END_REF].

This thesis is focused on applying and ameliorating the caspo-ts method to large scale experimental data of four breast cancer cell lines to identify cell line specific Boolean models.

More precisely, I apply the caspo-ts method to multi-perturbation time series data of four breast cancer cell lines (BT20, BT549, MCF7, UACC812) along with the traditional breast signaling network. I study how these cell lines derive different signaling behaviors under the same conditions. This behavior is important to study to unravel the remarkable heterogeneity among breast cancer types. It may help to design effective therapeutic strategies by having a better understanding of the underlying behavior of the system.

Contributions

1.2.1 The caspo-ts system applied to breast cancer Protein signaling networks are static views of dynamic processes where proteins go through many biochemical modifications such as ubiquitination and phosphorylation to propagate signals that regulate cells and can act as feedback systems. Understanding the precise mechanisms underlying protein interactions can elucidate how signaling and cell cycle progression occur within cells in different diseases, such as cancer. This knowledge may guide better drug designs.

In this work, I focused on computational identification of BNs representing protein signaling behavior using the caspo-ts [OPS + 15] method. I used four breast cancer cell lines of the HPN-DREAM challenge dataset [HHC + 16, HNJC + 17]. This dataset contains multiperturbation time series data of four breast cancer cell lines (BT20, BT549, MCF7, UACC812). I thoroughly constructed and refined the prior knowledge network (PKN) from public databases, such as Reactom [WDD + 14a], to cover the maximum number of proteins existing in the HPN-DREAM challenge dataset. I modeled a family of cell line specific BNs for the four breast cancer cell lines given this PKN.

My key findings suggest that this method is capable of constructing cell line specific Boolean models, which is extremely valuable given the heterogeneity of breast cancer due to many genetic modifications. An algorithm is implemented to analyze these cell line specific BNs to 1.2. CONTRIBUTIONS study similarity among these breast cancer cell lines. I have highlighted the common and distinct behaviors.

Further, I validated the inferred BNs using the testing dataset provided by the HPN-DREAM challenge. This dataset was not included while inferring BNs. My models have a Root Mean Square Error (RMSE) of 0.31 with respect to the traces for the testing data, providing an optimal fit to the testing data. Furthermore, I also validated the cell line specific Boolean models by comparing them to the canonical mTOR pathway. The obtained results are comparable to the top performing teams of the HPN-DREAM challenge. In addition, this approach can also be used as a complementary method to identify erroneous experiments. Related to this work, I have published an article in Plos Computational Biology journal [RPS + 18].

Extension of the caspo-ts system Diverse Boolean models

Since the ASP solver uses a backtracking algorithm to exhaustively generate BNs, it can lead to a situation where successive BNs share very similar properties. This can be problematic specially in the case of a large solution space where discovering or analyzing all BNs becomes computationally hard. To resolve this issue, a diverse enumeration scheme has been introduced. This feature has been implemented in caspo-ts and allows it to break up the clusters of similar BNs, hence generating diverse BNs. I refer to the modified caspo-ts as caspo-ts D .

I have demonstrated the results of the proposed approach on two different benchmark scenarios in systems biology: (1) an artificial dataset to model TCR signaling and (2) the HPN-DREAM challenge dataset to model breast cancer cell lines.

Results suggests substantial improvements of caspo-ts D in solution quality by discovering more signaling behaviors than caspo-ts. Moreover, caspo-ts D is able to find BNs in cases where caspo-ts is unable to find any. Related to this work, I have published an article in the Computational Methods in Systems Biology conference [RKR + 18].

Parallel model checking

The ASP part of the caspo-ts system over-approximates BNs. This over-approximation removes a large set of BNs that have no reachable traces, reducing the number of invalid BNs.

However, over-approximated reachability does not guarantee to reproduce all time series traces.

Hence, at the final step the caspo-ts system uses a model checker to check exact reachability of all (binarized) traces existing in the experimental data by the given BN. This is the most time consuming part of the caspo-ts system. The verification of this reachability is a PSPACE-hard problem and the computation time for checking reachability is highly variable depending on the BN under verification. It can take from an hour to months. Moreover, caspo-ts specifies all properties in one big specification, which can be slow to verify, especially in the case of large scale networks. I have improved this step, by splitting up the specification, to reduce the computational time of true positive BN detection. Results suggest substantial improvements in computation time.

Organization of the thesis

This thesis consists of six chapters.

Chapter 1 "Introduction" provides a broad overview and motivation behind this thesis. It precisely specifies the aims and contributions of this work.

Chapter 2 "Background and Related Work" provides background introduction on breast cancer, phosphoproteomic data, signaling networks, and computational modeling. We discuss four different computational methods to model signaling networks: ordinary differential equations, genetic algorithms, dynamic bayesian networks, and integer linear programming. We state the advantages and disadvantages of these methods as compared to the modeling approach (answer set programming) used in this work.

Chapter 3 "Answer Set Programming" introduces the modeling approach used in this thesis.

We define basic notations and explain them with examples. Then we discuss related work in the context of modeling signaling networks using answer set programming. We also discuss related work in the context of learning diverse solutions using answer set programming.

Chapter 4 "Computational Discovery of Dynamic Cell Line Specific Boolean Networks from Multiplex Time-Course Data" discusses the first contribution of this thesis. We describe in detail the caspo-ts system and its application on a real case study (HPN-DREAM Challenge).

We build signaling networks given phosphoproteomic data and prior knowledge networks. We used multiple criteria to evaluate the learned networks. We also highlight different characteristics (such as enumeration order, false positive rate) of the caspo-ts method, which are resolved in the following chapter.

Chapter 5 "Computing Diverse Boolean Networks from Phosphoproteomic Time Series Data" introduces the improved version of the caspo-ts methods. This chapter is related to the second contribution of the method. Here, we describe the new enumeration criteria to sample the solution space of the caspo-ts method. We also give a new algorithm to improve the computational time of verifying a solution using model checker.

Chapter 6 "Conclusions and Future Work" provides the summary and the future perspectives of this work.

Background and Related Work 2.1 Breast cancer

Breast Cancer is a remarkably complex and heterogeneous disease. It develops as a result of uncontrolled growth of abnormal cells due to genetic mutations. Mutation causes cells to multiply and divide chaotically. Mostly, this leads to the situation where multiple copies of abnormal cells give birth to a tumor. There exists a high diversity within the same type of breast tumors in terms of genomic alterations [RvJH + 13]. Moreover, there is also substantial difference within tumor bearing patients. This high level of diversity poses a great challenge for cancer therapy, and demands diverse clinical features [START_REF] Peri | Human protein reference database as a discovery resource for proteomics[END_REF].

Breast cancer types

Breast cancer can develop in different sites of the breast: the lobules, the ducts or the tissue in between them. The milk producing section of the breast is called lobules. Milk travels through the pathway from lobules to the nipple using ducts (see Figure 2.1). Depending on the area, breast cancer can be divided into two subtypes: (1) Non-invasive and (2) Invasive. With the non-invasive type, the cancer is completely confined to the ducts or lobules and does not spread to the surrounding connective tissues. The non-invasive breast cancer is further subdivided into two types: (1) ductal carcinoma in situ and (2) lobular carcinoma in situ. With the invasive type, the cancer breaks through the lobular and duct wall, and spreads into the surrounding connective tissues. There are further categories of the invasive breast cancer; the most common ones are: (1) invasive ductal carcinoma and (2) invasive lobular carcinoma [SDS + 10].

Cell lines

Cell culture is the process in which cells are cultivated in an artificial environment. Cells can be isolated through a variety of means. They can be derived from the living tissue and disaggregated by mechanical or enzymatic means. They can also be taken from an already established cell line. Cells are maintained and proliferated under controlled conditions, known as the primary culture stage. At this stage, cells are usually subcultured to provide more room for cultivation. After subculturing, cells capable of proliferation are selected and as with the outgrowth from the primary culture stage, give rise to a cell line. Cell lines developed through a primary culture usually have a genetically determined life span. However, some cell lines are transformed into immortal cells and can grow indefinitely under optimal conditions. This can happen because of variety of reasons, e.g., chemical induction or spontaneous occurrence [START_REF] Freshney | Culture of cells for tissue engineering[END_REF].

Cell lines provide a simple and powerful model for studying and analyzing breast cancer biology. They serve as an unlimited source of homogeneous self-replicating material, and are easy to handle, and replace in case of contamination. However, they pose a number of challenges too. Cell lines are susceptible to genetic and epigenetic drifts during their culture [START_REF] Birgersdotter | Gene expression perturbation in vitro-a growing case for three-dimensional (3d) culture systems[END_REF]. If a cell line is stored for a long time, a subpopulation can arise and cause phenotypic drift resulting in various clones within it [START_REF] Dai | Breast cancer cell line classification and its relevance with breast tumor subtyping[END_REF].

Despite of their challenging nature, cell lines have been widely used as an experimental model for in vitro studies [NCF + 06, BFA + 10, PZS +

Molecular classification of breast Cancer cell lines used in this thesis

As described above, breast cancer cell lines are categorized according to the three important receptors: ER (estrogen receptor), PR (progesterone receptor), and HER2 (human epithelial receptor 2) [START_REF] Dai | Breast cancer cell line classification and its relevance with breast tumor subtyping[END_REF]. In Table 2.1, we show the molecular classification of the four cell lines used in this study. In the following, we discuss some history and research work related to these breast cancer cell lines.

BT20

In 1958, the BT20 cancer cell line was isolated from a 74 year old woman. It is the first human breast cancer cell line and was established by Lasfargues and Ozzello [START_REF] Etienne | Cultivation of human breast carcinomas[END_REF]. It depicts the invasive breast cancer type. Even though it was the first breast cancer cell line, it has not been widely used to study breast cancer. Here, we cite some of the work which has been done using this cell line alongside other cell lines.

In the past, there was a lack of experimental models to study biological properties of human breast cancer. For this, they transplanted the BT20 cell line into nude athymic mice. Results

showed the development of tumors at the injection site. Their observations showed that nude athymic mice are suitable for studying human breast cancer [OSM + 74].

This cell line together with other breast and ovarian cancer cell lines have been used to study the HER2/neu assay sensitivity. According to the results, the BT20 cell line did not contain an over-expression of HER2 [RJC + 02].

A group of 41 breast cancer cell lines (including BT20) was studied to identify BRCA1

MCF7

In 1970, the MCF7 cancer cell line was developed using pleural effusion from a 69 year old female metastatic breast cancer patient. It was established by Soule and colleagues at the Michigan Cancer Foundation [SVL + 73]. MCF7 is the most studied breast cancer cell line.

It mimics several invasive breast cancers which expresses the estrogen receptor (ER). It is widely used to study estrogen response and resistance both in vivo and in vitro, as this cell line maintains a substantial level of ER, considering ER maintenance is not a trivial task. There have been around 25,000 publications related to this cell line till now.

MCF7 played a major role in the development of antibodies for ER positive tumors, since Green et al. developed the first monoclonal antibody to ER [START_REF] Greene | Monoclonal antibodies to human estrogen receptor[END_REF]. These antibodies led to the identification of cDNA clones which express ER mRNA. This aided in the cloning and sequencing of the ESR1 gene [WGG + 85]. These antibodies also helped in measuring ER levels in human breast cancer, thus serving as a guide on the use of hormone therapy for ER positive tumors. This cell line also expresses progesterone, glucocorticoid, and androgen receptors [START_REF] Horwitz | Mcf-7: a human breast cancer cell line with estrogen, androgen, progesterone, and glucocorticoid receptors[END_REF]. Hence, it is used as a valuable model to study other hormone response pathways.

An area of research focuses on estrogen based stimulation of MCF7 cells. Earlier studies were focused on how estrogen regulates the growth factor signaling [START_REF] Kent | Mechanisms of endocrine resistance in breast cancer[END_REF]. Recent studies have revealed that estrogen represses and induces a large number of genes simultaneously. This generates a complex network of alterations, which coordinate to change growth [CPM + 12].

MCF7 cells have also been used to study hormone resistance. Some researchers have developed hormone resistant variants of MCF7 cells by either chronic exposure to anti estrogens or estrogen withdrawal. Initially, this slowed down the cell growth, but finally growth resumed.

It has been shown that estrogen deprived cells are highly sensitive to estrogen stimulation and express substantial levels of ER [JSB + 98]. These studies have identified transcriptomic and epigenetic alterations, which eventually lead to changes in growth factor signaling. Then prominent signaling pathways were pharmacologically targeted in these cell lines to study the response to targeted therapies. Results showed that different cell lines have different sensitivities. This suggests that the heterogeneity of this disease can be studied through these triple negative cell lines and can help with designing effective preclinical treatments [LBC + 11].

In [GMN + 12], a group of 25 triple negative breast cancer cell lines were analyzed to identify similarities between cell line and triple negative breast cancer. These cell lines were studied on three molecular levels: genomic, transcriptomic, and epigenomic. They categorized cell lines in three groups. Two groups consisted of ER negative cell lines, while one encompassed of three ER negative and all ER positive cell lines The first two groups agreed with the existing knowledge, while one group was not consistent with existing studies. They further extended the characterization of these breast cancer cell lines. This provided valuable knowledge about suitability of a particular cell line for modeling different features of the breast cancer disease [GMN + 12].

Recently, a new classification criteria was proposed to categorize breast cancer cell lines [SSK + 17]. The relationship between RNA, TP53 mutation status and protein expression was also quantified. Their analysis revealed heterogeneity within cultures of established cell lines.

They compared their finding with other studies, to help guide the selection of cell line models for in vivo and in vitro studies [SSK + 17].

UACC812

In 1988, the UACC812 cell line was isolated from a 42 year old women. This They suggested to completely block the HER network and inhibit ER to develop an optimal therapy [WMG + 11].

Phosphoproteomics experimental dataset

Cancer is both a genomic and a proteomic disease at a functional level. A genetic defect may ultimately lead to alterations in a protein signaling network Since 2011, an annual meeting is held to provide a platform for scientists and researchers for exchanging ideas related to RPPA technology.

RPPA is a high-throughput technology that performs quantitative measurements of hun-dreds of proteins in biological samples (clinical and preclinical). It consists of micro-blots of protein lysates from multiple samples of cell lines or tissues on a single array where each sample is described by at least one spot. Each array is incubated by one antibody to identify the corresponding protein expression across multiple samples simultaneously. Each array can contain thousands of samples. For high-throughput measurement of many proteins, multiplexing on multiple arrays of the same set of lysates is performed using different antibodies (see Figure 2.

3) [START_REF] Huang | Reverse phase protein arrays in signaling pathways: a data integration perspective[END_REF].

Computational modeling of protein signaling networks

Network modeling has been widely used for studying phosphoproteomic data, yielding important insights into protein interactions, functions, and evolution. Figure 2.4 shows the work-flow of the computational modeling. Computational modeling bridges the gap between traditional biology and high-throughput datasets. A model is constructed using different formalisms like mathematical modeling, stochastic search methods, bayesian networks, integer linear programming, answer set programming, and then integrated with the experimental data.

Afterwards, models can be validated through testing data. These testing data can be generated by traditional biology experiments. If a model's predictions coincide with the testing data, then it can be used to predict novel events. These novel events can be further validated through experiments before suggesting novel biological insights. Network modeling is an important part of this thesis. In the following, we discuss the general frameworks that can be used for network discovery.

Ordinary differential equation models

Ordinary differential equations (ODEs) based approaches have been widely used for modeling complex and dynamic signaling networks. ODEs represent the interaction among various CHAPTER 2. BACKGROUND AND RELATED WORK molecules (such as proteins or genes). In ODEs, the various biological species (such as proteins) are represented as variables, where each variable has an equation reflecting its dynamic evolution over time. Here, we describe two types of ODEs:

1. Law of Mass Action, 2. Hill Function.

Here, we show how to formulate biochemical reactions (such as phosphorylation) using above mentioned equations.

Law of mass action

This law states that the rate of reaction is proportional to the chance of a collision of reactants, which in turn, is proportional to the concentrations of the participating molecules to the power of the molecularity (such as the number in which they enter the specific reaction)

[JYL + 17]. For example, we can define the molecular interaction between A, B and C by the following equation:

A + B k + ---- k - C, (2.1)
where A and B are reactants, C is the product, and k + and k -are the forward and the reverse kinetic rate constants . The following ODEs can be used to derive the change in the concentration of A, B, and C over time:

d[A] dt = d[B] dt = k + [C] -k -[A][B], (2.2)
d[C] dt = k + [A][B] -k -[C]. (2.3)

Hill function

In the context of protein signaling networks, Hill functions can be used to represent the state of proteins such as activation or inhibition. In the following, Hill functions can be used to describe the change in protein expression over time. The Hill function is associated with each protein involved in the signaling network.

dy dt = M i=1 f + (x i ) + N j=1 f -(x j ) -y * d y , (2.4
)

f+ - (x) = k x y x + - nx H + - nx x + x + - nx , (2.5)
where y denotes the concentration of activated protein, x i represents ith protein (i = 1, 2, ..., M ), which activates protein y , and x j is the jth protein (j = 1, 2, ..., N ), which inhibits protein y. f+ -(x), represents the activating or inhibiting profile induced by protein x, respectively; k x y denotes the activation or inhibition rate; H x represents the microscopic dissociation constant and n x is the Hill coefficient; d y denotes the degradation rate of protein y [JYL + 17].

Iadevaia et al.

In [ILM + 10], a computational approach is proposed to model the IGFR signaling network in the MDA-MB231 breast cancer cell line, using a set of ODEs with the law of mass action.

The IGFR signaling network representing activation and inhibition interactions, was mapped into a set of 77 chemical reactions. For example, the following equation was used to represent the activation of protein PIK31 by IRS1,

P IK31 + IRS1 * k + ---- k - P IK31 * , (2.6)
where k + and k -are the forward and the reverse kinetic rate constants. These set of chemical reactions were mapped into a system of 127 ODEs.

d(pB) dt = K b * B, (2.7)
where p represents the phosphorylation of the protein.

The authors also took into account the trafficking of proteins between the cell membrane and cytoplasm, since HER3 can take hours or days to transfer to cell membrane. This slow transfer can be caused by the feedback loops involved in the pathway or DNA reprogramming (cells rewrite their DNA to survive). Two versions of the model were created to account for the transfer of protein between cell membrane and cytoplasm. DNA reprogramming was modeled by a long delay, and feedback was modeled as a reaction rate following the law of mass action. They implemented the process of protein transfer between cell membrane and cytoplasm by defining two compartments. One compartment contained HER2 and HER3, and other contained internal HER3. Transfer between these two compartments was modeled using diffusion differential equations. Parameter estimation was performed using simple search algorithm, allowing the models to fit experimental data. In this work, they described the detailed analysis of one model which has a feedback loop mechanism to show the slow transfer of HER3 between two compartments. The authors applied Gefitinib on this model to study its effects. For this, the model was augmented with a chemical reaction (which follows the kinetic law of mass action), representing that Gefitinib binds to HER2 and its dimers. Results showed pHER3 and AKT recovered their levels even after persistent inhibition of HER2, and HER3 transferred from cytoplasm to the cell membrane. After they inhibited the trafficking between compartments and applied Gefitinib. Results showed pHER3 and AKT did not recover their levels, and HER3

did not transfer to the cell membrane. Another simulation was performed where AKT was activated externally. Results showed decrease in pHER3 level and HER3 is transferred from the cell membrane to the cell. They kept the parameters where simulation results agreed with the experimental data. Finally, they proposed and simulated a treatment scheme where AKT was activated for two days, which caused drainage in HER3 from the cell membrane, and then Gefitinib was applied which inhibited the HER2/3 and reduced transferring. Their simulation results suggest that the proposed scheme works better than the traditional scheme where Gefitinib was applied alone. Moreover, low level of AKT was maintained, eventually leading to the apoptosis of cancer cells. Then this mathematical model was trained against experimental data (imaging data) to build a cancer cell line specific pathway. To build the liver cell specific pathway, the generic pathway was trained against the cues response signaling dataset. This model was represented with 52

ODEs (following the law of mass action) and 83 parameters. For parameter estimation, an objective function was defined to minimize the distance between the simulation and experimental results. The genetic search algorithm was adopted as a two stage approach for optimal parameter estimation. During the first stage, the algorithm was repeated 50 times to select the best parameters with minimal error between simulation and experimental results. Then the model was simplified by removing links (containing non-best parameters). During the second stage, best parameters were refined by repeating the search algorithm for a simplified model.

To perform validation, a leave-one out approach was used. Simulation predictions agreed with experimental results. They analyzed the effects of 27 kinase inhibitors and identified 6 key inhibitors. Out of these, the PF02341066 kinase inhibitor was identified as a proper inhibitor for suppressing cancer while avoiding damage to the liver cells. They also analyzed the effect of combined inhibitions by selecting 4 out of 6 effective key kinase inhibitors. Their findings revealed the threshold for each kinase inhibitor in a combination to predict the expected combined effect.

Conclusion

Mathematical modeling offers useful and powerful tools to model small-scale networks.

Models elucidated using mathematical modeling are complex and require explicit specifications of kinetic parameters of the system; parameter estimation becomes computationally intensive as networks grow larger [WMG08, SK13, MTH +

Stochastic models

Stochastic modeling frameworks are used to model signaling networks accurately, by offering a strategy to cope with the uncertainty and noise inherent in biological processes. Different stochastic methods can be used to model signaling networks such as genetic algorithms.

Genetic algorithms

A genetic algorithm is a bio-inspired search method, which follows the principal of natural selection and evolution. It starts with the belief that the offspring inherit properties of their parents. If parents have the best fitness then their offspring will be better than parents. This is an iterative process, and finally the fittest individuals will be found. The key findings suggest that the trained networks have fewer interactions as compared to the protein signaling network, but has higher false positive, false negative rate. New links or interactions were suggested which improved the fitness of BNs to the experimental data. The resulting models were more predictive than the protein signaling network due to the elimination of nonfunctional links.

Mishra et al.

A hybrid approach based on genetic algorithm and ODEs was proposed by Mishra et al.

[MBC + 15] to model cell fate decisions and cancer signaling pathways. The authors started with building an apoptosis network (N 1 ) representing the signaling activities inside a normal cell. The network N 1 consisted of 22 proteins, and was constructed by exploring literature knowledge. Twenty two ODEs were formulated for each protein to simulate the N 1 network.

These ODEs were solved using ode45 (Matlab library function), to infer time series for each protein. Some parameter values were inferred directly from the literature, while others were selected randomly. From the N 1 network, the cancer network (N especially with their inherent ability to deal with noise. However, quality of the solution in case of genetic algorithms is highly dependent on the initial population, well written objective function, and the number of iterations. Similar to our method (caspo-ts), the genetic algorithm based method generates a family of solutions. In comparison to caspo-ts, stochastic search methods cannot generate a complete set of solutions, hence they cannot guarantee a global optimal solution.

Dynamic bayesian networks

A bayesian network is a directed acyclic graph G, representing a probabilistic relationship between random variables X i where i = 1 . . . n. These random variables represent nodes and the probabilistic dependencies among them is represented by edges. The edges are described by the joint probability distribution P (X 1 , . . . , X n ). The joint probability distribution of the graph has the following general form:

P (X 1 , . . . , X n ) = Π N i=1 P (X i = x i ||X j = x j , . . . , X j+p = x j+p )
where

X j is a parent of X i in G [FLNP00].
Dynamic bayesian networks [DGH92, DGH + 95, MR02] are an extension of bayesian networks in which a directed graph represents how random variables X i evolve over time. Dynamic bayesian networks can incorporate feedback loops. In the context of modeling time series data for protein signaling networks, proteins are represented by random variables X i , time is represented by T j where j = 1 . . . n. The dynamic bayesian networks described here follow the assumption that time can flow forward only, therefore the value of the variable X i at t + 1 is only dependent on value of the variable X i at t [FMR98, Hus03]. The dynamic bayesian networks can be referred to as static bayesian networks unrolled through time where each variable is represented at multiple time points [HLM + 12].

There are two steps to model bayesian networks: (1) structure learning and (2) parameter learning. The structure learning step involves finding a set of candidate directed acyclic graphs which best explain the data. Then these graphs are scored against some bayesian metric (Bayesian Dirichlet equivalence or Bayesian Information Criteria). The graph with highest scoring metric is selected. In the parameter learning step, the conditional probability distribution of each node is learned. The PKN is taken into account by assuming a particular distribution over a directed acyclic graph.

Hill et al. 2012

Hill et al. [HLM + 12] proposed a framework to learn cell specific protein signaling networks from phosphoproteomic data given a PKN. They started by building a structure of the network which best explains the data. The PKN is taken into account using an informative prior distribution on network structure. Then the posterior distribution over the network is calculated using the marginal likelihood and prior distribution. The marginal likelihood employed here penalizes models with complex parameters, hence reducing the model complexity. work consist of probabilities associated to each edge in all these 32 networks. For validation, the changes in these context specific networks under inhibition were compared with the true descendants existing in the canonical pathway. Based on this, a true and false positive rate was calculated to draw the area under the operating curve (AUROC). It was observed that on average 8 proteins showed changes in one context which were not observed in other context.

Experimental validation was performed for 78 % of these observations (104 out of 134 observations). To verify robustness of the approach, some part of the data (between 1 to 6 time points) for all proteins was deleted and replaced with adjacent time points. Initial and final time points were kept outside from random removal of the data. The robustness of the learning approach to the prior network was verified by adding and deleting an edge in the PKN. They generated a set of testable hypotheses by discovering some novel interactions. They discovered 235 novel interactions and some (six) of these interactions were validated using western blot analysis. Their key findings suggest that the signaling networks varies according to biological background. Heterogeneity of inferred networks strongly supports the fact that existing computational methods should be improved to take into account the context specificity.

Conclusion

The afore-mentioned approaches provide a useful way to model signaling networks in specific contexts. However, the dynamic bayesian network approach (like any other data driven approach) suffers from missing variables (intermediate nodes) due to limitation of measurements available in the experimental dataset. This hinders the mechanistic interpretation of the learned networks. Mostly, approaches based on dynamic bayesian network modeling assume the homogeneity of network structure and parameters through time. The softening of this assumption can be computationally expensive due to rapid increase in graph space and number of parameters [HLM + 12]. Moreover, predictions discovered through these networks are relatively less accurate. Also, these networks do not provide the activation and inhibition information [SHM + 15]. Generally, bayesian approaches are computationally expensive. They also cannot perform well with small datasets produced by protein microarrays[NSS + 08].

Integer linear programming

Integer linear programming (ILP) [START_REF] George | Constraint Classification for Mixed Integer Programming Formulations[END_REF] is used to solve optimization problems. An ILP solver takes a description of such a problem consisting of a set of constraints and an objective function as input, calculating a solution that satisfies all constraints and minimizes the objective function. The ILP approach uses a set of variables, and a set of equalities and inequalities to represent the constraints. The variables are required to have integer values. A typical ILP program is expressed as follows:

minimize c T x (2.8) subject to Ax ≤ b (2.9)
where c and b are vectors, A is a matrix, x ≥ 0, and x ∈ Z n where Z is the set of integers.

Here, we give an example of the coloring problem for a graph G = {V, E} with nodes V = {1, . . . , n} and edges E ⊆ V × V . This problem states that adjacent nodes of graph G cannot have the same color. We assume that there are n colors in total. The binary variable y k where k ∈ {1, . . . , n} is used to denote whether the color k is used for coloring. Another variable x ik is used to denote the color k for node i. This problem can be expressed in ILP as follows:

minimize n k=1 y k (2.10) subject to n k=1 x ik = 1 i ∈ V (2.11) x ik -y k ≤ 0 i, k ∈ V (2.
12)

x ik + x jk ≤ 1 (i, j) ∈ E, k ∈ V (2.13) 0 ≤ x ik , y k ≤ 1 i, k ∈ V (2.14) x ik , y k ∈ Z i, k ∈ V (2.15)
The first three constraints ensure that each node is colored and the last two constraints are used to restrict x ik and y k to binary values 1 . Their results revealed the main known alterations in the pathway and also several unknown drug effects. Their key finding suggests that lapatinib, erlotinib, and gefitinib have similar kind of effects, very few alterations differentiating them, while sorafenib drug induce totally different signaling alterations. + 11] extended the above mentioned ILP approach to build an extended pathway to link signaling behavior to a cellular response such as cytokine release or cell growth. They constructed two kinds of datasets for normal and cancer hepatocytes. One dataset (D 1 ) consisted of measurements of 16 phosphorylated proteins under approximately 50 perturbations (7 stimuli and 5 inhibitors). The other dataset (D 2 ) contained measurements of 33 cytokines under same set of perturbations. The learning process started with the construction of the canonical pathway via literature search, which covered the key phosphoproteins and stimuli. Then the ILP framework [MMS + 09] was modified to add edges (to link proteins with cytokines) not existing in the canonical pathways. This extended network consisting of canonical and non-canonical edges was optimized to fit both experimental datasets (D 1 , D 2 ). The optimization of the objective function was done in three main steps: 1) by removing the canonical edges which contradict signaling dataset (D 1 ), 2) removing non-canonical edges which contradict (D 2 ), and 3) removing the edges which have no effect on the network to decrease the overall network size. The ILP solver was used to optimize the objective functions to generate 100 different solutions for normal and cancer hepatocytes (Huh7). These solutions differed by 10% in the value of their objective function. The extended pathway for normal hepatocytes consisted of a total of 47 edges (19 canonical and 28 non-canonical edges). The extended pathway for cancer hepatocyte (Huh7) consisted of 43 edges (26 canonical and 17 non-canonical edges). These pathways are significantly different from each other, supporting the fact that Huh7 cells are not as responsive to stimuli (Toll Like Receptor) as normal hepatocytes are. The predictions performed through these models were also experimentally validated. The validation of the learning approach was performed by assessing the model's sensitivity to the canonical pathway, experimental design, and measured data. To verify sensitivity to the generic topology, 10 % of edges were replaced with random edges. The results showed that optimized pathways are sensitive to the canonical pathway. For verifying sensitivity to the experimental design, 5 to 50 % of experiments were left out. Results suggest that fitness error increases with the removal of number of experiments. The sensitivity to the measured data was checked by replacing part of the data-points with random numbers between 0 and 1. The fitness error also increased with the removal of number of data-points. + 09] provided an ILP based framework to study the effects of drugs on cell specific signaling networks. They successfully identified drug induced alterations in a network. Similar to our method (caspo-ts), the ILP based approach [MMS + 09] is topology dependent and cannot determine the alterations outside the constructed networks. Melas et al.

Mitsos et al.

The authors in [MMS

Melas et al.

Melas et al. [MMM

Conclusion

Mitsos et al. [MMS

[MMM + 11] proposed an extended ILP framework to study the cellular responses (growth, phenotypic, cytokines) under different signaling events. Their framework is based on the assumption that signaling and cellular dataset must be measured under same experimental conditions.

The proposed approach is quite sensitive to changes in topology and datasets as compared to the caspo-ts. Similar to caspo-ts, ILP solvers enumerate a family of solutions and can find global optimal solutions. However, open source ILP solvers are not as efficient as commercial ILP solvers are [MMS + 09]. All in all, ILP provides a useful paradigm to model signaling networks. However, ILP lacks a simple programing or modeling language which makes ILP coding a difficult task to write and maintain.

Conclusion

ODEs are suitable for detailed modeling of small systems where a complete description is available to model the system. Stochastic methods (more precisely genetic algorithms) cannot guarantee global optimal solutions as ILP and ASP offers (shown by [VGE + 12]). However, the lack of an accessible programming language makes ILP hard to model the problem (depicted in graph coloring problem). Moreover, open source ILP solvers are not adequate to model problem of inferring BNs from time series data [MMS + 09]. This thesis is focused on inferring BNs using an ASP based modeling approach, by combining traditional signaling networks with complex phosphoproteomic time-series data generated by RPPA in [HNJC + 17, HHC + 16]. This approach allows to build large scale networks without having to deal with large scale parameters as in case of ODEs, stochastic and bayesian approaches. Since ASP is an important part of this thesis. Therefore in the next chapter, we are introducing the modeling approach (ASP) used in this thesis, to formulate the problem of finding BNs from PKN and phosphoproteomic time series data. In this thesis, we use ASP to solve the combinatorial optimization problem of learning BNs by confronting phosphoproteomic data with PKNs. Available ASP solvers use high performance solving techniques, allowing us to enumerate all optimal solutions for a given problem. Indeed, there exist other solving techniques like satisfiability testing (SAT) [START_REF] Biere | Handbook of satisfiability[END_REF], Integer Linear Programming (ILP) [START_REF] Schrijver | Theory of linear and integer programming[END_REF] and Constraint Programming (CP) [START_REF] Dechter | Constraint processing[END_REF]. ASP solving techniques build upon SAT solving techniques. However, SAT-based systems lack ASP's rich modeling language, which forces the user to write complex programs for generating encodings, making it difficult to maintain and modify them. Another modeling difference is that the problem instances and the problem encodings are inseparable in SAT. Further, SAT encodings can be easily translated into the ASP programs in a modular way, but not vice versa [Nie99]. In general, the translation of an ASP program to a SAT encoding can lead to an exponential space explosion unless the language is extended [LR06]. Another difference between SAT and ASP solving is that the ASP semantics is non-monotonic while SAT is monotonic. This makes it easy to model concepts such as reachability and inertia in ASP. Constraint Programming and Integer Linear Programming also rely on monotonic foundations. Moreover, these formalisms focus on solving problems with constraints over real-valued variables or integers. However, there are extensions to ASP (for example clingcon), which allow for using similar constraints in ASP encodings. This allows for modeling a large class of problems with a natural and expressive modeling language using state-of-the-art solving techniques [SW18].

Answer Set Programming

Basic ASP syntax

ASP borrows basic terminologies from classical logic. Variables are denoted by strings starting with capital letters, such as U , V , Color , etc. Symbols for predicates are represented using strings starting with lowercase letters, like for example, p, q, node, etc. Function symbols are denoted using strings starting with a lowercase letter, too, for example, f , g, h, etc. 1 Each predicate and function symbol is associated with an arity written as p/n and f /n where n ≥ 0 gives the arity. A term is a variable, integer, or function f (t 1 , . . . , t n ) such that f is a function symbol of arity n and each t i is a term. An atom has form p(t 1 , . . . , t n ) such that p is a predicate of arity n and each t i is a term. Function symbols with arity zero are called constants. Furthermore, the parenthesis after a function or predicate are omitted if the arity is zero. Atoms or their negations are called literals. Terms and atoms without variables are called ground terms [START_REF] Gebser | Answer set solving in practice[END_REF]. Table 3.1, describes the basic connectives used to form logic programs.

1. In general the alphabets of predicate and function symbols can overlap because this information can be derived from context within a logic program. Both alphabets are kept disjoint in the logic programs presented in this thesis.

h ← b 1 ∧ • • • ∧ b m ∧ ¬b m+1 ∧ • • • ∧ ¬b n where h is an atom, 0 ≤ m ≤ n,

ASP extensions

Next, further extensions [START_REF] Simons | Extending the stable model semantics with more expressive rules[END_REF] to logic programs are described, which are frequently used in practice and ease modeling problems with ASP.

Choice rules

A choice rule has form

{h 1 ; . . . ; h o } ← b 1 ∧ • • • ∧ b m ∧ ¬b m+1 ∧ • • • ∧ ¬b n
where 0 ≤ o and each h i is an atom. Unlike with the normal rule above, a choice rule can be used to prove any subset of the atoms h 1 to h o whenever its body is satisfied.

Example 3.2.2. One nice aspect of choice rules is that they allow for writing more compact programs. Remember that in the previous example we needed two atoms to model that a switch is on or off. Choice rules allow for a more compact and easier to read representation.

Integrity constraints

An integrity constraint has form

← b 1 ∧ • • • ∧ b m ∧ ¬b m+1 ∧ • • • ∧ ¬b n
It cannot be used to prove atoms, but instead to remove candidate solutions satisfying its body. The first rule can prove arbitrary subsets of {switchOn, light}. Then the second integrity constraint discards all solutions where the light is on but the switch is not. Similarly, the last integrity constraint prunes all solutions where the switch is on but the light is not. Hence, we

obtain the same answer sets as in Example 3.2.2.

Cardinality aggregates

A cardinality aggregate can be used to express that a set of atoms has to have a certain cardinality. A cardinality aggregate has form

l {e 1 ; . . . ; e n } u
where n ≥ 0, the lower bound l and upper bound u are integers, and each element e i has form l 0 : l 1 , . . . , l m for m ≥ 0 and each l j is a literal. An aggregate element e i is satisfied whenever the literals l 0 to l m are satisfied. The head literals of satisfied elements are collected in a set whose cardinality has to be between l and u (inclusive) for the cardinality aggregate to be satisfied. One of the bounds can be omitted; if the lower bound is omitted, then the cardinality has to be smaller than u and greater than l, otherwise. Cardinality aggregates can occur in both rule heads or rule bodies. In a rule head they can be used to prove the literals l 0 of true aggregate elements. This leads to another way to model the light domain example. Each switch has an associated weight and the sum of weights of active switches has to be larger than 2 for the light to turn on. Let's consider the solution {on(s 1 ), on(s 3 ), light}.

This solution satisfies the first and third aggregate element, hence, we obtain the set of tuples {(1, s 1 ), (1, s 3 )}. Both tuples in this set have weight 1 and thus their sum is 2, which satisfies the lower bound. This allows us to derive atom light and we indeed have a solution.

So far we have not considered variables in aggregates. The above example can actually be written more compactly by introducing a variable in the aggregate.

{on(s 1 ); on(s 2 ); on(s 3 )} ← light ← 2 sum{1, S : on(S)}

Variable S in the aggregate is a local variable because it occurs only between the braces {. . .}.

In this case, this does not lead to multiple ground rules but in the three aggregate elements as in the preceding program by substituting s 1 , s 2 , and s 3 for S.

Optimization statements

An objective function to maximize or minimize a lexicographically ordered tuple of sums can be described in the ASP program using optimization statements. For example, a minimize 3. If t 0 is not an integer it is simply ignored.

statement can be expressed as:

minimize{e 1 ; . . . ; e n }
where n ≥ 0 and e i has form t 0 @p 0 , t 1 , . . . , t o : l 0 , . . . , l m for o ≥ 1, m ≥ 0, p 0 is a term, each t k is a term and each l i is a literal. The accumulation of integers works almost exactly as with sum aggregates. The only difference is the optional priority p 0 (which defaults to 0) in the elements e i . Sums are calculated for each priority individually and then ordered by priority and compared lexicographically when determining optimal solutions.

{max (a, 20); max (b, 30)} 1 ←

1 {min(c, 30); min(d, 10)} ← maximize{W @1, X : max (X, W )} minimize{W @2, X : min(X, W )}
In the above program, priority levels indicate that minimization (with priority level 2) is more important than maximization (with priority level 1). Solution {max (b), min(d)} induces the tuple of sums (10@2, 30@1), which is indeed an optimal solution for the above program because at most one atom over max /1 and at least one atom over min/1 can be true.

Graph coloring problem

Here, we give a logic program for the graph coloring problem (defined for ILP in Section 2.4.4). We first describe the input graph in the form of facts. Predicate node/1 determines the available nodes and predicate edge/2 the edges between the nodes. The example graph in Figure 4 . It ensures that colors are selected successively beginning with the first color (color (1)). Having a rule like this helps when enumerating all optimal solutions by drastically reducing the resulting set of solutions. Note that there are still more isomorphic solutions that could be pruned with additional rules. Finally, the optimization statement 3.8 is used to describe the objective function. This objective function minimizes the number of colors used and corresponds to Equation 2.10 in the ILP formulation.

{color (C)} ← node(C) (3.4) 1 {assign(U, C) : color (C)} 1 ← node(U ) (3.5) ← edge(U, V ) ∧ assign(U, C) ∧ assign(V, C) (3.6) ← color (C) ∧ ¬color (C -1) ∧ node(C -1) (3.7) minimize{1, C : color (C)} (3.8)
By specifying the above constraints, an ASP solver computes 6 optimal stable models. One of the computed stable models is shown below and depicted in Figure 3.3:

{assign(1, 1), assign(2, 3), assign(3, 3), assign(4, 2), assign(5, 1), assign(6, 2)}

4.

Here an important feature of ASP is used: terms can contain arithmetic operations. Such operations are evaluated during the grounding phase, e.g., C -1 evaluates to 2 if C is three.

Related work

Here, we discuss the related work of finding optimal Boolean Networks (BNs) from PKN and phosphoproteomic data using ASP. We also discuss the related work of finding diverse solutions of a logic program having a large solution space. Furthermore, they started with building two PKNs (P 1 and P 2 for mid and large scale case studies) representing growth and inflammatory signaling from literature. These PKNs were compressed to remove non-observable and non-controllable nodes, creating four compressed networks. Then, hyper-graphs compatible to these four compressed networks were created. A hyper-graph refers to the extension of the graph with AND and OR gates. For each compressed network, five hyper-graphs (of three different sizes) were selected randomly to generate 240 in silico datasets. The size of the hyper-graph is calculated by summing up the number of source nodes for all child nodes.

Learning Boolean networks

Results. Videla et al. inferred 8 global optimal models for the mid-scale case study and 2 global optimal models for the large-scale case study. They also executed the same problem with the CellNOpt system ([SRAE + 09]) and the comparison is drawn in Table 3.2. The biggest advantage of ASP over genetic algorithms is the calculation of an optimal score and the enumeration of all optimal solutions having this score. They also studied the impact of size on the solution space (number of models). They concluded that if they relax the size limit of the model then the size of the solution space increases exponentially. They also discovered that both approaches are able to find all global optimal solutions in 240 datasets. However, results suggest that the models of minimal sizes are best identified with an ASP solver. Computation time was also compared with CellNOpt and results showed that ASP solvers are faster than CellNOpt by 5 orders of magnitude. Case study. The authors illustrated the working of this approach on three case studies. For each case study they used three different versions of PKNs. Hyper-graphs compatible with these PKNs were used to generate in-silico datasets with 10 perturbations. The first case study (C 1 ) represents EGF-TNFα signaling ([MTH + 12]). Case study C 1 consists of 13 nodes and 16 edges. The second case study (C 2 ) represents TCR signaling ([KSRL + 06]). From this case study, three different case studies (C Results. Their findings suggest that caspo-ts is able to find the first optimal solution within seconds while the full enumeration can take longer. They also observed that the full enumeration of optimal solutions in the case of C They also compared the results with the previously proposed approach ([VGE + 12, GVE + 13]) using case study C 1 . They used the approach of Videla et al. to learn BNs with one time point and compared the fitness error with caspo-ts. Results show that caspo-ts outperforms Videla et al. approach to find an optimal model with the data in most of the experiments.

Conclusion

Regarding inferring BNs, the methods [SRAE + 09, VGE + 12] (of Videla and Saez) restrict themselves to learn BNs from two time points (start, end), assuming the system has reached an early steady-state when measurements are performed. This assumption prevents one from capturing interesting characteristics like loops as shown in [MTH + 12]. To overcome this issue, Ostrowski et al. proposed the caspo-ts method. The first part of this thesis is in the context of applying and ameliorating the caspo-ts method on a real case study of four breast cancer cell lines.

Learning diverse solutions

Another important part of my work is focused on the improvement of the caspo-ts search method in case of a large case study. The caspo-ts system uses a backtracking algorithm (described in Chapter 5, Section 5.2.3) to exhaustively generate solutions. It can lead to a situation where successive solutions share very similar properties. This can be problematic, especially in the case of a large solution space where discovering or analyzing all solutions becomes computationally hard. To resolve this issue, a diverse enumeration scheme has been introduced. In the following, we focus on discussing the work on finding diverse or similar solutions using ASP; similar work for constraint programming and SAT solving can be found here [START_REF] Hebrard | Finding diverse and similar solutions in constraint programming[END_REF][START_REF] Nadel | Generating diverse solutions in sat[END_REF].

Eiter et al.

Method description. Eiter et al. [START_REF] Eiter | Finding similar or diverse solutions in answer set programming[END_REF] proposed an approach to study similar or diverse solutions using ASP. They have introduced online/offline methods to enumerate diverse solutions.

The offline methods are based on computing all solutions in advance, while the online methods follow an incremental solving approach. For finding diverse or similar solutions in the case of the offline methods, they start with a complete set of solutions. These solutions form the nodes in a graph and distance-labeled edges between nodes represent the distances between solutions. The idea is to find a clique in this graph using edges with distances greater than k to find diverse solutions and edges with distances less than k to find similar solutions. This clique is computed using ASP.

They introduced three variations of online methods. The first variation encodes everything in one program consisting of three parts: (1) compute n solutions in parallel, (2) compute distances between solutions, and (3) remove answers where the distance between any two solutions is not equal to k. The second online method does not calculate solutions in advance but rather updates the solution set with a new solution in each run with a solution that has exactly distance k from all solutions found so far. The third variation of the online methods is similar to the second method but modifies the ASP solver clasp to incorporate a distance function in it.

Results. Eiter et al. illustrated this approach on phylogeny reconstruction explaining the evolutionary development of biological species or other entities (taxonomic units). They built diverse or similar phylogenies (solutions) for Indo-European languages. The methods are used to find diverse or similar solutions among eight phylogenies.

The first offline and online methods guarantee to find an optimal solution since they consider the whole solution space of a given problem, while online methods two and three cannot guarantee an optimal solution given the fact that the solving process is based on the first solution found. In their experiments, the offline method three performed better than all other methods.

Zhu et al.

Method description. Zhu et al. [ZT13] considered Eiter's approach in a more general setting where optimal solutions are calculated based on user preferences. Their approach first identifies optimal solutions and then computes a set of diverse or similar solutions. They studied three methods to solve this problem.

The first method is an iterative algorithm using an additional program to test optimality. This program decides if a computed answer set is optimal according to the preference rules.

The current answer set is updated if the next solution is better than the current one. At the end, an optimal answer set is returned. Next, to find an alternative optimal answer set, their algorithm finds another answer set which is also optimal and different from the solutions found so far. Similar or diverse solutions are found by adding distance constraints. This process continues iteratively to find further solutions.

The second method is based on introducing a modification to ASP solver clasp so that it checks the optimality of the answer sets using a tester program during the search. Otherwise, it is equivalent to the previous method.

Both of the previous methods are based on multiple calls to the ASP solver while the third method models all requirements in one disjunctive logic program. In general, such programs are more complex to solve than normal logic programs.

Results. Zhu et al. illustrated the working of their approach on four computational problems: finding one optimal solution, an alternative optimal solution, three similar optimal solutions, and eight diverse optimal solutions. They compared their results on the basis of ranked (where one preference is prioritized over the other) and unranked preferences.

The results show that the iterative method performs best with unranked preferences when finding optimal solutions while the third method (based on modeling) performs well for finding similar or diverse solutions. However, for ranked preferences the third method is slower than the iterative method.

Romero et al.

Method description. Romero et al. [START_REF] Romero | Computing diverse optimal stable models[END_REF] introduced a framework (asprin 2) to compute diverse or similar optimal solutions in case of computational problems with preferences. They introduced three methods to find diverse optimal solutions: (1) enumeration, (2) replication, and (3) approximation.

The enumeration method first finds all solutions for a given problem and then computes diverse solutions.

The replication method converts a logic program into an equivalent logic program with a maxmin preference; diverse optimal solution are then computed for the rewritten program.

Finally, the approximation method iteratively identifies diverse optimal solutions w.r.t. previously found solutions. They implemented three different variations of this method. The first variation is based on maximizing the minimum distance to previously found solutions. The second variation calculates a partial interpretation by maximizing the minimum distance to previously found solutions and then calculates a model close to this partial interpretation. And the third variation uses a heuristic to derive diverse optimal solutions.

Results. They illustrate the working of their approach on datasets from six different domains: biological network repair, metabolic network expansion, crossing minimization, circuit diagnosis, design space exploration, and timetabling.

The enumeration and replication methods are not as effective as the approximation methods since enumeration and replication methods traverse through an exponential number of optimal solutions to find an exact solution. The approximation techniques do not promise to provide an exact result but performs well in case of the large scale problems. Results suggest that the third variation of the approximation method performs better than the others in terms of computation time. They also compare the quality of the diversification where the third method based on heuristics performed worse than the other methods. A clear trade-off emerged between quality and performance.

Conclusion

The best performing approach of Eiter et al., [START_REF] Eiter | Finding similar or diverse solutions in answer set programming[END_REF] modifying the clasp solver, cannot guarantee to produce optimal solutions as the solving process is biased toward the first enumerated solution. Moreover, they did not consider optimization problems. Hence, their methods are not directly applicable to caspo-ts, since it enumerates optimal (subset minimal) solutions in order to produce simpler and more relevant solutions. Zhu et al. [START_REF] Zhu | On optimal solutions of answer set optimization problems[END_REF] can handle programs with preferences and also proposes modifications to the clasp solver but their approach becomes impractical as the search space increases. Therefore, we decided to extend the approach of [START_REF] Romero | Computing diverse optimal stable models[END_REF] for computing optimal diverse solutions in ASP. The novelty of this extension is that our approach uses heuristics for both the computation of optimal (subset minimal) solutions and the diversification. In the thesis, we refer to the modified caspo-ts as caspo-ts D . modifications. Furthermore, we also validate the cell line specific regulatory mechanisms using the inferred Boolean models.

Caspo-ts framework

We chose the caspo-ts method [OPS + 15, OPS + 16] for the inference of BNs. This method was tailored to handle phosphoproteomic time series data. The input of the method consists of a prior knowledge network (PKN) and normalized phosphoproteomic time series data under different perturbations. It generates a family of BNs whose structure is compatible with the PKN and that can also reproduce the patterns observed in the experimental data. The caspo-ts workflow is shown in Figure 4.1. It consists of two main steps, ASP solving and model checking. Here, we formally describe the inputs, outputs, solving and model checking processes of the caspo-ts framework.

Prior knowledge network

A prior knowledge network (PKN) is modeled as a labeled (or colored) directed graph are used to bound the system and also serve as interaction points of the system, these nodes can be experimentally stimulated, e.g. cellular receptors. Inhibitors are those nodes which remain inactive or blocked over all time points of the experiment by small molecule inhibitors.

(V, E, σ) with V = {v 1 , v 2 , . . . , v n } the set of nodes, E ⊆ V × V
Stimuli and inhibitor nodes take Boolean values {0, 1}, representing the fact that the node was stimulated (1) or inhibited (0). Readouts are experimentally measured given a combination of stimuli and inhibitors. They usually take continuous values in [0;1] after normalization. Unobserved nodes are neither measured nor experimentally manipulated. In this study, we use the term perturbation to refer to the combination of stimuli and inhibitors, similarly to other studies such as [Hei16, HHC + 16, HNJC + 17]. Thus, the set of all possible perturbations is the set of sets P = (S ∪ I). The PKN is one of the main inputs of the caspo-ts framework.

Phosphoproteomic time series data

Phosphoproteomic time series data are obtained by measuring the temporal changes in phosphorylated proteins under a perturbation (Figure 4.1). Without loss of generality, we as-Figure 4.1 -Caspo-ts workflow. Caspo-ts receives as input data a prior knowledge network (PKN) and a discretized phosphoproteomic dataset. In this example the phosphoproteomic data consists of two perturbations involving akt (inhibitor) and hgf (stimulus): 1) akt = 0, hgf = 1, and 2) akt = 1, hgf = 0. A black colored perturbation means the inhibitor or stimulus was perturbed (1) while white represents the opposite (0). Readouts are specified in blue and describe the time series under given perturbations. Using this input data, caspo-ts, performs two steps: ASP solving and model checking. In the ASP solving step: (i) a set of BNs compatible with the PKN is generated, (ii) afterwards an over-approximation constraint is imposed upon each candidate BN to filter out invalid BNs that do not result in an over-approximation of the reachability between the Boolean states given by the phosphoproteomic dataset, and finally (iii) BNs are optimized using an objective function minimizing the distance to the experimental measures. The ASP step also introduces repairs in some data points of the time series that added penalties to the objective function. These corrected traces will be given to the model checker. In the model checking step, the exact reachability of all the (binarized and corrected) time series traces in the family of BNs is verified.

Binarized & CorrectedTraces Time Series Corrections

Caspo-ts sume that the time series data are related to the observation of m ≤ n nodes for the nodes {v 1 , . . . , v m } (so the nodes {v m+1 , . . . , v n } are not observed). The observations consist of normalized continuous values: a time series of k data points is denoted by T P = (t 1 P , . . . , t k P ), where P ⊆ S ∪I is a perturbation and t j ∈ [0; 1] m for 1 ≤ j ≤ k. These data will be discretized in order to link it with further BNs discovery (see the ASP solving and model checking steps). Here the values between zero and one of three proteins are shown in different colors. In this figure, we see the time series of one readout protein (blue) under a perturbation of one stimulus (green) and one inhibitor (red). Stimuli have value 1 and inhibitors have value 0 across all time points of an experimental perturbation. Readouts take continuous values in [0;1] after normalization. In some phosphoproteomic datasets an inhibitor can also act as a readout protein, which means that there are perturbations where it will be measured.

Boolean network

A Boolean Network (BN) [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF][START_REF] Inoue | Logic programming for boolean networks[END_REF] is defined as a pair B = (N, F ), where -N = {v 1 , . . . , v n } is a finite set of nodes (or variables/proteins/genes), -F = {f 1 , . . . , f n } is a set of Boolean functions (regulatory functions) f i : B k → B, with B = {0, 1}, describing the evolution of variable v i .

A vector (or state) x = (x 1 , . . . , x n ) captures the values of all nodes N at a time step, where x i represents the value of the node v i , and is either 1 or 0. There are up to 2 n possible distinct states for each time step. Next, we define the transition x → x between two states of a BN. If there is no update for node v i then x i = x i . If there is an update for node v i then the state of a node v i at the next time step is determined by x i = f i (x 1 , . . . , x n ). Note that usually only a subset of the nodes influence the evolution of node v i . These nodes are called the regulatory nodes of v i . The state of each node can be updated in a synchronous (parallel) or asynchronous fashion. In the synchronous update schedule, the states of all nodes are updated, while in asynchronous update schedule, the state of one node can be updated at a time. The work presented in this article is independent of the update schedule routine, hence any number of nodes can be updated at a time. We use x → * y to say that state y can be reached from state x with an arbitrary number of steps. Notice that the BN is the main output of the caspo-ts method.

ASP solving PKN expansion

Given a PKN and a phosphoproteomic dataset, a family of candidate BNs, compatible with this PKN, is exhaustively enumerated. The BNs are represented by Boolean formulas in Disjunctive Normal From (DNF), i.e., as a disjunction of conjunctive clauses. A clause can be seen as a reaction, where the proteins represented positively are available, and the proteins represented negatively are absent. A Boolean formula in DNF encompasses all possible reactions to update the value of a protein. We refer to these BNs as subset minimal BNs.

Here, we use two properties to describe the BNs compatible with the PKN [VGE + 12].

Intuitively, the first property states that if there is an edge in a BN then this edge must exist in the PKN as well. The second property states that the BNs use the smallest DNF formulas possible, in the sense that no conjunctive clause can be removed from a DNF formula without changing the Boolean function it represents. It means a Boolean formula should have as few nodes as possible.

Evidence property:

Given a PKN (V, E, σ) and v ∈ V , a logical formula φ in DNF has an evidence in (V, E, σ) with respect to v if and only if for every propositional variable w that occurs positively (resp. negatively) in φ, there exists an edge (w, v) ∈ E and ((w, v), 1) ∈ σ (resp. ((w, v), -1) ∈ σ).

Redundancy property:

Given a logical formula φ in DNF, with φ = j≥1 c j where each c j is a conjunctive clause, φ is a redundant formula if and only if for some k, l ≥ 1 with k = l and some logical conjunction r it holds that c k = c l ∧ r.

Over-approximation

After generating a set of BNs (compatible with a PKN), a constraint is imposed to check the reachability of the states from another within state graph of a BN. Since this reachability is a computationally hard problem (PSPACE-complete) [START_REF] Cheng | Complexity results for 1-safe nets[END_REF], we use an over-approximation criteria. This constraint is imposed upon each candidate BN to filter out invalid BNs, that do not result in an over-approximation of the reachability between the Boolean states given by the phosphoproteomic dataset.

A meta-state u = (u 1 , u 2 , . . . , u n ) is a vector of dimension n over non-empty subsets of B, noted M = {{0}, {1}, {0, 1}}; the set of meta-states is M n . Meta-states characterize a set of Boolean states: a state x ∈ B n belongs to a meta-state u, written x ∈ u, iff each Boolean component x i belongs to the set u i . Given a state x, we use x for the corresponding meta-state ({x 1 }, . . . , {x n }). We define the transition relation u ⇒ v between the meta-states u and v as follows: u = v and v = (u 1 , . . . , u i ∪ {f i (x) | x ∈ u}, . . . , u n ) for some 1 ≤ i ≤ n.

In [OPS + 16], it has been shown that if y is reachable from x (x → * y) then there exists a meta-state u such that y ∈ u and x ⇒ * u. This definition is further refined to describe the necessary condition for reachability called support consistency. A state x is support consistent with state y denoted by x * y, if and only if there exists a meta-state u with x ⇒ * u such that y ∈ u and for all 1 ≤ i ≤ n either y i = x i , or y i = x i and u i = {0, 1}, or y i = x i , u i = {0, 1}, and there exists z ∈ u such that f i (z) = y i .

If state y is reachable from state x (x → * y) then x * y. Since we are using the overapproximation criteria, it is possible that some of BNs may fail to reproduce some trajectories of the time series data. These BNs are called false positive (FP). To filter out the false positive BNs, exact model checking is applied.

Optimization

The goal is to generate a BN that can reproduce the experimental data as well as possible.

For this purpose, BNs are optimized by defining an objective function minimizing the distance between the actual time series, T P , and the predicted time series, Y P :

RMSE = 1 m * k * |P| m i=1 k j=1 P ∈P ((t j P ) i -(y j P ) i ) 2 , (4.1)
where m is the number of observed nodes, k is the number of time points, and P ∈ S ∪ I is the set of perturbations. If a BN can verify the time series traces of experimental data then it will have minimal RMSE . In addition, the optimization step highlights the data points in the time series which added penalties to the RMSE. Such data points are automatically corrected before the model checking step.

All the analyses described in this step are performed using ASP, using the clingo 4.5.4

solver [START_REF] Gebser | Clingo= asp+ control: Preliminary report[END_REF]. This solver guarantees finding optimal solutions, and all BNs outputted by the ASP solver step will be identically optimal.

Model checking and true positive BNs

From the previous step, a set of optimal BNs that over-approximate the phosphoproteomic time series data is produced. This set of BNs is verified with exact model checking to detect true positive (TP) BNs. TP BNs are guaranteed to reproduce all the (binarized) trajectories under all perturbations by verifying exact reachability in the BN state graph. For this, we have used computational tree logic (CTL) implemented in the NuSMV 2.6.0 [CCG + 02], which is a symbolic model checker. Properties are written using the existential (E) and finally (F ) operators. The existential quantifier asserts that the property holds for at least one path. The finally operator asserts that the property will hold eventually. They are specified in the following manner:

Exp t 0 e → EF (Exp t 1 e ∧ EF (Exp t 2 e ∧ EF (. . . EF (Exp tn e ) . . . ))), (4.2)
where Exp t i e refers to the state under perturbation e ∈ P at time point t i . The above property (4.2) states that there exists a way to reach time point t n from time point t 0 in the perturbation e. We specify properties for all perturbations (P). If a BN verifies all these properties then it is a true positive (TP) BN. If any of the properties is dissatisfied then the BN is called a false positive (FP).

Computation of root mean square error

The caspo-ts method can also generate the RMSE for the data absent (testing data) from the training data. It generates over-approximated traces of the BN and compares them with the traces existing in the testing data (see Algorithm 1). For this computation it uses the same equation as given in Equation 4.1. There are two types of RMSEs returned by caspo-ts -discrete and model. The discrete RMSE is imposed by the discretization of the method. Since we use a discrete learning approach, our recovered traces will be in {0,1} and this introduces an error with respect to continuous measurements in [0;1]. The model RMSE refers to the learned BN error with respect to the normalized time series data; that is, the model RMSE is at least as large as the discrete RMSE. When the difference between these two is zero then the Note that, this Algorithm 1 is given to illustrate a general idea about RMSE calculation.

The ASP program does more than generating the over-approximated values of proteins of BNs.

It also optimizes the minimum distance between over-approximated Boolean traces and actual time series traces. It guarantees to capture the minimum distance that is possible of Boolean traces of BN to the actual time series traces.

Caspo-ts logic program

Here, we provide the pseudo logic program to describe the input data given in the caspo-ts modeling framework Section. The logic program in written in the ASP language. ASP is a powerful declarative logic programming language for knowledge representation and reasoning

[Lif08]. The basic idea is to encode the problem using a non-monotonic logic program and then feed it into the ASP solver, which computes the solution of the problem in the form of models (also known as answer sets). Note that we only provide pseudo encodings here, please refer to [OPS + 16] for details.

Modeling PKN and data

Facts, rules and constraints are the building blocks of ASP programs. Here we use facts to describe the inputs. The PKN (V, E, σ) is described by the following facts:

node(v) ← for v ∈ V edge(u, v, s) ← for (u, v) ∈ E and ((u, v), s) ∈ σ
Facts over the predicate node are used to model the nodes existing in the PKN. A fact edge(u, v, s)

is used to describe an edge between node u and v with sign s.

For each perturbation P ∈ P and phosphoproteomic time series T P , we have the following facts:

clamped(P, v, 0) ← for v ∈ P ∩ I clamped(P, v, 1) ← for v ∈ P ∩ S obs(P, j, v i , s) ← for s = (t P j ) i , 1 ≤ j ≤ k and 1 ≤ i ≤ m
Facts over predicate clamped set the value of node v in perturbation P to either 0 or 1 depending on whether the node belongs to inhibitor set (I) or stimuli set (S) respectively. Facts over predicate obs set the value of the readout node v i in perturbation P at time point j to value s.

Over-approximation

To filter the compatible BNs with the over-approximation criteria, we start with generating Boolean values for each node of the BN with the following rule:

1{state(P, T, U, (0; 1))}1 ← P ∈ P, node(U ), T = 1..k.

Here, state predicate is implemented as a choice rule. It has four arguments. The first argument represents the perturbation P , the second represents the time-point T , the third is the node U , and the fourth represents the guessed value for the the node U . States are subject to further constraints which are omitted here for brevity.

Further, we calculate the meta-states from the given states of the BNs with the following rules:

meta(P, T, U, V ) ← state(P, T, U, V ) update(P, T, U, V ) ← "node U takes value V if its formula can evaluate to V" meta(P, T, U, V ) ← update(P, T, U, V )

The first rule initializes the meta-states from the state predicate. The second rule selects possible updates by evaluating the Boolean functions associated with the nodes. The third rule updates the current meta-state with the update predicate. We omit the exact encodings here and refer to the literature [OPS + 16] for details.

The following integrity constraints encode the support consistency condition described in the over-approximation Section:

← state(P, T + 1, U, V ), ¬meta(P, T, S, V )

← state(P, T + 1, U, V ), ¬update(P, T, U, V ), update(P, T, U, W ), V = W

The first rule ensures that the state at time-point T + 1 is in the meta-state reachable from timepoint T . The second rule implements the condition in the third item in the definition of support consistency described in the over-approximation (see Section 4.2.4). It states that if a state at time-point T + 1 has value V and in the update predicate it has value W then there must exist an update in-between to recover the final value V .

Objective function

Here, we describe the objective function to minimize the difference between the phosphoproteomic time series data and the guessed time series for the BN. The following rules declare the objective function by using minimize directive:

minimize{M -49, P, T, U : obs(P, T, U, M ), state(P, T, U, 0), M ≥ 50} minimize{50 -M, P, T, U : obs(P, T, U, M ), state(P, T, U, 1), M < 50}

If the conditions after the " : " do not evaluate to true, then the guessed time series does not match with the corresponding measurements (phosphoproteomic data) and a penalty is accumulated by the minimize statement. Note that since ASP does not support continuous values, we converted them to integers in the range from 0 to 100 by using the formula: x * 100 .

Graph similarity measure

Here, we present a graph similarity measure in order to check the variability among the family of BNs generated by caspo-ts. This measure serves to compare the reactions existing in a gold standard network (A) with a family of BNs (B) and is based on the Jaccard similarity coefficient which measures the similarity of these models (see Algorithm 2). The Jaccard index between A and M is defined as length of the intersection divided by the union:

J(A, M ) = | A ∩ M | | A ∪ M | = | A ∩ M | | A | + | M | -| A ∩ M | (4.3)
4.5 HPN-DREAM challenge case study

Data acquisition

The DREAM portal provides unrestricted access to complex, pre-tested data to encourage the development of computational methods. In this study, we are focused on the HPN-DREAM challenge, which was motivated by the fact that the same perturbation may lead to different signaling behaviors in different backgrounds, making it necessary to build a model which can perform unseen predictions (absent from the learning data). The main goal of the HPN-DREAM challenge is to learn signaling networks efficiently and effectively to predict the dynamics of breast cancer [HHC + 16, HNJC + 17].

Learning data

Reverse Phase Protein Array (RPPA) quantitative proteomics technology was used for generating the dataset of this challenge. The measurements focus on short term changes on up to 45 proteins and their phosphorylation over 0 to 4 hours. The HPN-DREAM dataset includes temporal changes in phosphorylated proteins at seven different time points (t 1 = 0min, t 2 = 5min, t 3 = 15min, t 4 = 30min, t 5 = 60min, t 6 = 120min, t 7 = 240 min). The learning data consists of four cancer cell lines (BT20, BT549, MCF7 and UACC812) under different perturbations (8 stimuli and 2 or 3 inhibitors). The number of perturbations varies from 24 to 32 depending on the cell line. In each cancer cell line up to 45 phosphorylated proteins are measured against different sets of perturbations over multiple time scales.

Testing data

Test data is available for assessing the performance of networks learned from the learning data. The HPN-DREAM portal provides testing data for four cancer cell lines (BT20, BT549, MCF7 and UACC812) under different perturbations (8 stimuli and 1 inhibitor). They contain gold standard datasets of time series predictions of up to 45 proteins having the same time scale as learning data [Hei16, HHC + 16, HNJC + 17]. The number of perturbations varies from 7 to 8 depending on the cell line. These data will be used to test the quality of the BNs inferred by caspo-ts.

Normalization

The protein measurements of the HPN-DREAM challenge have variable ranges as depicted in Maximum value based normalization was used to set the measurements between a common scale, i.e., 0 and 1, in order to assign activation or inactivation values to variables or species of the BN. Equation (4.4) describes the formula used for the normalization. Given time series T P , we obtain time series T P :

(t P j ) i = (t P j ) i max{(t Q l ) i | Q ∈ P, 1 ≤ l ≤ k} , (4.4)
where i ∈ {1, . . . , m} are the measurements, j ∈ {1, . . . , k} are time-points, and P ∈ P are the perturbations. Here (t P j ) i represents the value of protein i under perturbation P at timepoint j and the denominator denotes the highest value of protein i under all perturbations and time-points.

Data preprocessing

The learning and testing datasets used in this study were extracted from the HPN-DREAM challenge and correspond to time series protein measurements upon different perturbations of four breast cancer cell lines (UACC812, BT20, BT549, and MCF7) [HNJC + 17, HHC + 16].

Since readout signals are measured on variable ranges depending on the protein, a normalization step was necessary. The learning dataset had a few noisy, inconsistent and incomplete time series data points. The caspo-ts system identified these inconsistencies existing in the time series data. A recurrent experimental inconsistency observed was an oscillation in the protein signal upon experimental inhibition of the same protein.

To resolve the above mentioned issues, we performed the following data preprocessing steps on the learning dataset:

1. Set the protein values between a common scale, i.e., 0 and 1, using a maximum-valuebased normalization scheme (4.4).

2. For time point 0 the expression of some readout proteins under some perturbations was not available. Thus, control experimental readings were used as the time point 0 for such proteins.

3. In some cases readout measurements were duplicated for the same time point. To solve this noise issue we chose one time point arbitrarily.

4. We removed inconsistent perturbations where the protein AKT was inhibited and was having a dynamic behavior as a readout protein.

5. We considered only perturbations with complete time series data, since guessing the missing time points automatically with caspo-ts for this case study will be computationally expensive.

The experimental errors pointed in steps 2-5 were raised as warning or exceptions by caspots. Steps 1 to 5 were applied to the learning dataset. Only step 1 was applied to the testing dataset. After removing perturbations with inconsistent behaviors or incomplete time series from the learning dataset, we had 15, 13, 13 and 18 perturbations for MCF7, BT20, BT549 and UACC812 cell lines respectively measuring 23 readouts. Figure 4.5 shows the preprocessed dataset for cell line BT20.

Prior knowledge network

Our first result was to build the consensus protein signaling network associated with the HPN-DREAM challenge dataset. The structure of the prior knowledge network (PKN) was generated by mapping the experimentally measured phosphorylated proteins (HPN-DREAM dataset) to their equivalents from literature-curated databases and connecting them together

Cell line specific Boolean networks

In this section, we show the generated BNs for each cell line. For this, we used caspo-ts to learn the BNs from the PKN (Figure 4.6) and the phosphoproteomic data of four breast cancer cell lines -BT20, BT549, MCF7, and UACC812. We inferred a family of cell line specific BNs for each cancer cell line. As explained in the caspo-ts modeling framework (section 4.2), the caspo-ts method produces BNs fulfilling two criteria, (i) satisfaction of the over-approximation criteria and (ii) optimality with respect to the RMSE objective function. ASP-optimal solutions were fast to collect, their computation time ranged from 36 seconds to 3 minutes depending on the cell line, as shown in the Table 4.1.

Afterwards, these ASP-optimal BNs were given to the model-checker for further verification. This second step is more complex and we put a restriction on the computation time of 7 days for each cell line. The number of verified BNs varies from one cell line to another, depending on a number of factors such as the number of perturbations, the order of answer sets in the solutions space, and the perturbation order. The total number of verified ASP-optimal BNs within the 7 day time-frame were 188, 231, 52, and 150 for the BT20, BT549, MCF7 and Figure 4.6 -Breast cancer signaling pathway. This figure shows the reconstructed signaling network from a combination of databases. An arrow shows the positive regulatory relationship between two proteins, while a T shaped arrow indicates inhibition. Green nodes are stimuli, blue nodes are readouts, white nodes are unmeasured or unobserved, and blue nodes with a red border represent inhibitors and readouts at the same time. Note that in the node labels, we have added the phosphorylation sites to the protein names in order to connect them to the experimental measurements. UACC812 cell lines respectively. We obtained 72, 191, and 21 true positive BNs for BT20, BT549, and MCF7 cell lines respectively with an optimal fit to the data. For the UACC812 cell line, we were unable to obtain true positive BNs within the 7 day time limit for verification.

Hence, we kept the first 20 BNs from the 150 ASP-optimal BNs for the UACC812 cell line.

Union of these BNs are shown in Figures (4. can be seen that 23% of the Boolean functions are shared among BT549 and MCF7, and also between BT20 and UACC812. BT20 shares the least number of Boolean functions (15%) with BT549. This table revealed pronounced differences among different cell lines of breast cancer.

We also analyzed the diversity of Boolean functions among the family of BNs within the same cell line. The similarity among Boolean functions from BT20 (0.73) and MCF7 (0.63) is higher than the ones from BT549 (0.43) and UACC812 (0.46) cell lines.

Heterogeneity among cell lines

There are a total of 69 distinct Boolean functions shown in Figure 4.12 along with their respective frequencies. It is interesting to note that the B549 and UACC812 cell lines have The aggregated graph for all cell lines. Blue, red, green and brown colors have been used for each cell line BT20, BT549, MCF7 and UACC812, respectively. The nodes are connected by logic gates (AND or OR) to their direct predecessors. Edges are used to show influences (→ for positive and for negative).

An AND gate is depicted by a small black circle where the incoming edges correspond to the inputs of the gate. An OR gate is depicted by multiple incoming edges to the node. A different color scheme is used to represent different types of nodes. The green color is for stimuli, the red for inhibitors, the blue for readouts, and the white for unobserved nodes. Black edges denote common hyper-edges across cell lines; the thickness of the black hyper-edge denotes the number of cell lines sharing this hyper-edge. Interestingly, only 4% of the Boolean functions are shared in all cell lines and 99% of these shared functions have the same frequency. In this figure, there is only one Boolean function which is frequent in 3 cell lines and has a lower frequency in BT20.

Biological literature related to the Boolean functions discovered by caspo-ts

The caspo-ts method revealed that cell line specific reactions are clustered around the AKT, MAPK3, and PIK3R1 proteins. PI3K is an important factor for cancer development in HER2 amplified cancers (UACC812) as compared to non-HER2 amplified (BT20, BT549 and MCF7) cancer cell lines. The HER2 amplified tumor is more aggressive than other types of tumor (see Chapter 2 Section 2.1.1). We can see from Figs 4.7, 4.8, 4.9 and 4.10 that PIK3R1 exists in all cell lines but is rather more connected in the UACC812 cell line with 10 incoming edges, while in others with only 1 incoming edge. The PIK3R1 node in UACC812 (Figure 4.10) has a centrality measure of 0.37 while in the other three cell lines the centrality measure is less than 0.11. The centrality measure is used to quantify the most important node within a network i.e., the number of times a node has been used as a bridge (along the shortest path) to connect to other nodes in the network [START_REF] Abboud | Subcubic equivalences between graph centrality problems, apsp and diameter[END_REF].

It has been established that P1K3R1 (the regulatory unit of PI3K) plays an important role in suppressing tumors (breast, lung, ovarian, bladder, and prostate) [SWF + 05, TWK + 10, SHGAL + 08]. Recently, it has been found that PIK3R1 is mutated in 3% of breast cancer cell lines [N + 12]. Therefore, it is worth studying the impact of the PIK3R1 regulatory unit in breast cancer.

Evaluation

The performance of the caspo-ts method is evaluated using three criteria: 1) RMSE calculation using a typical learning and testing data approach, 2) random data comparison, 3) AUROC (Area Under the Operating Curve) score. The BNs are learned using the learning dataset only.

The prediction accuracy is evaluated by comparing the RMSE of trajectories in the testing dataset with those recovered by the learned networks (see Equation 4.1 in Section 4.2.4). To check how our method performs in case of random time series, we calculated the RMSE score for random data and compared it with the one obtained using the learning and testing data.

Next, the validity of these networks was verified by comparing them with the canonical MTOR Figure 4.12 -Heterogeneous Boolean functions. The Boolean functions are represented on the y-axis and the frequency of each Boolean function is shown on the x-axis. A Boolean function, or hyper-edge, is of the form node ← expr, where node is the receiver of the Boolean clause expr in the BN. In the Boolean clause, the not operator is represented by a "!" symbol and the AND operator by a "+" symbol. The disjunction of clauses is represented by multiple reactions upon the same receiver node.

AKT2<-AKT1 P90<-GSK3A HER3<-AR+HER2 p70<-AKT1 BCL2L1<-BAD p70<-mTOR HER3<-HER2 GSK3B<-AKT1+GSK3A+P90 PGR<-ER BAD<-!AKT1+!AKT2+BCL2L1+MAPK3 GSK3B<-AKT1+P90 GSK3B<-GSK3A BAD<-!AKT1+!AKT2+BCL2L1+MAPK3+!P90 ER<-PGR+PIK3R1 ER<-AR+PGR+PIK3R1 AKT1<-IRS1+PIK3R1+STAT3+mTOR PIK3R1<-FGF1 Src<-PIK3R1+SERUM STAT3<-MAPK3 JNK<-IRS1+PIK3R1 MAPK3<-AKT1+NRG1 MAPK3<-AKT1+IGF1+NRG1 MAPK3<-AKT1+EGFR MAPK3<-EGFR+IGF1+NRG1 MAPK3<-AKT1+EGFR+NRG1 MAPK3<-AKT1+EGFR+IGF1 MAPK3<-EGFR+NRG1 AKT2<-AKT1+EGFR AKT2<-AKT1+STAT3 AKT2<-AKT1+EGFR+mTOR AKT2<-AKT1+EGFR+STAT3 AKT2<-AKT1+STAT3+mTOR AKT2<-AKT1+PIK3R1+STAT3 AKT2<-AKT1+EGFR+PIK3R1 EGFR<-PIK3R1 AKT1<-EGFR+GSK3A+STAT3+mTOR PIK3R1<-EGF+FGF1 IGF1<-PIK3R1 MAPK<-IRS1+MEK1 STAT3<-EGF+EGFR MAPK3<-NRG1 NRG1<-EGF+EGFR STAT3<-EGF AKT2<-EGFR+PIK3R1 AKT2<-EGFR+PIK3R1+STAT3+mTOR P90<-GSK3A+MAPK MAPK<-MEK1 HER3<-AR AKT1<-PIK3R1+mTOR BAD<-!AKT1+!AKT2+MAPK3+!P90 IGF1<-IRS1+PIK3R1 PIK3R1<-EGF+HER3+IGF1+NRG1 PIK3R1<-HER3+IGF1+NRG1 AKT2<-PIK3R1 PIK3R1<-HER3+NRG1 AKT1<-PIK3R1+STAT3 SMAD3<-EIF4E+GSK3B PIK3R1<-HGF PIK3R1<-HER3+IGF1 PIK3R1<-EGF+HER3 PIK3R1<-EGF+IGF1+NRG1 PIK3R1<-EGF+HER3+IGF1 PIK3R1<-EGF+HER3+NRG1 Src<-SERUM STAT3<-EGF+EGFR+MAPK3 NRG1<-EGF AR<-SMAD3 EIF4E<-!INSULIN SERUM<-IRS1+PIK3R1 Boolean Functions Frequency BT20 BT549 MCF7 UACC812 0 1 0.8 0.6 0.4 0.2
GSK3A<-!AKT1 MEK1<-c-Raf Rb<-!CASP3 CASP3<-CASP9 mTOR<-IRS1 IRS1<-IGF1 CASP9<-!AKT1 GSK3B<-AKT1+GSK3B Boolean Functions 0 0.2 0.4 0.6 0.8 1 Frequency BT20 BT549 MCF7 UACC812
signaling pathway using two parameters, i.e., true positive rate (TPR) and false positive rate (FPR).

Validation using root mean square error criteria

The goal is to verify that the BNs learned using caspo-ts applied to the HPN-DREAM dataset are able to recover trajectories that do not exist in the learning data. For this purpose, we used experimental testing data to check the specificity of the trajectories of the proposed networks. These testing data were provided by the HPN-DREAM challenge organizers. Table 4.3 shows the corresponding RMSE in case of learning and testing data. It can be seen that the inferred BNs are able to produce the discrete trajectories without any error in the learning dataset for all cell lines. These BNs can reproduce trajectories of the learning data with an error of 0.3464, 0.3498, 0.3207, and 0.3464 in BT20, BT549, MCF7, and UACC812 cell line respectively. It is encouraging to see that the inferred BNs are able to recover the discrete testing trajectories without any error in MCF7, and with a minimal error of 0.0009, 0.0106, and 0.0094 in BT20, BT549, and UACC812, respectively. These BNs can reproduce trajectories of the testing data with an error of 0.3302, 0.3113, 0.2772, and 0.3178 in BT20, BT549, MCF7, and UACC812 cell line, respectively.

We also compared the RMSE score with the top two best performers of the HPN-DREAM challenge. We got the top position with an RMSE score of 0.31 as compared to their RMSE scores of 0.47 and 0.50. Notice that in comparison to other HPN-DREAM challenge methods based on Bayesian inference, Regression, and Granger Causality among others, caspo-ts does not make new predictions but it checks the recoverability of the testing trajectories with the inferred BNs (see Section 4.2.6). was different and therefore the ASP-optimal BNs did not recover the full discretized time series data. Overall these results demonstrate that our learning process is robust and is able to handle noisy time series data.

Validation using MTOR canonical pathway

To perform the validation of the structure of the BNs, we calculated a set of standard nodes from our PKN which are downstream nodes of MTOR and belong to the canonical MTOR pathway. We then evaluated how many of these standard nodes are also downstream nodes of MTOR in the learned BNs. In the following, the set of downstream nodes of MTOR in the learned BNs is referred to as inferred set. The inferred set is specific to each cell line. True positive rate (TPR) and false positive rate (FPR) are defined by Equation (4.6) and Equation (4.7) respectively:

T P R = T P T P + F N , (4.6)
F P R = F P F P + T N . (4.7)
Here, TP is the number of nodes in the intersection between standard and inferred sets, FP is the number of nodes in the inferred set but not in the standard set, FN is the number of nodes in the standard set but not in the inferred set and TN is the number of nodes which are not in the standard set nor the inferred set.

Discussion

In this chapter, we built cell line specific signaling networks for the DREAM time series dataset of 4 breast cancer cell lines (BT20, BT549, MCF7, and UACC812) using caspo-ts. suggest that the behavior of cell line specific signaling networks is highly variable even under the same perturbations, agreeing with the heterogeneity of breast cancer and specifically with previous analysis on this data [HNJC + 17]. The inferred Boolean models of each cell line were analyzed to identify commonalities as well as discrepancies. Moreover, these inferred models can be executed computationally to identify potential drug targets or to see the effect of unseen perturbations. The predictive power of these models can be increased with improvements in protein interaction databases and comprehensive experimental data.

We have discovered 38% of the cell line dependent behaviors as compared to the 33% of the HPN-DREAM challenge winner [START_REF] Daniel E Carlin | Computational evaluation and derivation of biological networks in cancer and stem cells[END_REF]. We have implemented an algorithm to analyze the variability among cell lines and observed pairwise similarities among these cell lines. The similarity index varies from 15% (BT20 & BT549) to 23% (MCF7 & BT549, BT20 & MCF7).

We have analyzed the similarity among the family of BNs of the same cell line as well, which varies from 43% to 73%. We have evaluated the accuracy of our method with RMSE and AUROC scores. The average RMSE of the inferred BNs was 0.31 placing caspo-ts in first place in comparison with the top scores reported in the HPN-DREAM challenge. Various choices made during this study may have an impact on the final score. The caspo-ts method allowed us to remove noisy and faulty experiments, leaving us with the reliable experimental settings only.

Here, we made the choice to use only the reliable experiments of the learning dataset instead of using all experimental settings. Also, we did not observe all 45 proteins as we could not find connections in our PKN for all the studied proteins, leaving us with approximately 23 proteins for each cell line.

Nonetheless, the obtained results are quite promising, making caspo-ts a good candidate computational method for learning models given time series datasets and a prior knowledge network. In addition, caspo-ts can be used to pinpoint the errors in the experimental data. In particular, we discovered four experiments where the protein AKT was inhibited and had a dynamic behavior as a readout protein. Our work therefore provides a novel approach to show erroneous experiments, which is crucial and complementary to current approaches. Finally, the HPN-DREAM dataset contained some noisy readings of experiments. Noisy experimental data reduces the efficiency of computational methods by increasing the variability among constructed Boolean models. To overcome this, we suggest to build automated methods to filter out the noisy experiments. This approach provides a step forward in building cell line dependent networks in the case of phosphoproteomic data.

Conclusion

We have applied the caspo time series (caspo-ts) approach, which is a combination of logic programming and model checking, over the time series phosphoproteomic dataset of the HPN-DREAM challenge to learn cell specific BNs. The learned BNs can be used to identify the cell specific topology. Our analysis suggests that caspo-ts scales to real datasets, outputting networks that are not random with a lower fitness error than the models used by the 178 methods which participated in the HPN-DREAM challenge. On the biological side, we identified the cell specific and common mechanisms (logical gates) of the cell lines.

Since caspo-ts uses an ASP solver to enumerate BNs, in the resulting sequence of solutions similar BNs are typically clustered together. This can be problematic for large scale problems where we cannot explore the whole solution space in reasonable time. In the next chapter, we discuss how we resolved this issue by sampling to randomly select BNs from the solution space, which allows for shuffling the order in which solutions are enumerated [START_REF] Romero | Computing diverse optimal stable models[END_REF]. We implemented this by dynamically modifying the heuristic of the ASP solver at execution time.

We show an application of caspo-ts with diversification on the UACC812 breast cancer cell line and TCR signaling network. We obtained TP BNs for UACC812 cell line using extended caspo-ts.

Computing Diverse Boolean Networks from Phosphoproteomic Time Series Data

Introduction

In this chapter, we introduce the new enumeration scheme for the caspo time series (caspots) method. The article related to this work is published in the International Conference on Computational Methods in Systems Biology [RKR + 18].

The caspo-ts system uses an over-approximation criteria to learn candidate BNs, which can lead to some false positive (FP) BNs. These BNs are not guaranteed to reproduce all traces of the time series data. To resolve this issue, it uses model checking to filter out FP BNs, by checking the exact reachability of these BNs w.r.t. time series data. The caspo-ts method uses the clingo ASP solver [START_REF] Gebser | Clingo= asp+ control: Preliminary report[END_REF], which is able to exhaustively enumerate all solutions. The clingo solver by default uses an enumeration scheme in which, once a solution is found, it backtracks to the first point from where the next solution can be found. This typically leads to the situation where successive solutions (BNs) only change in a small part. As a result, caspo-ts may enter a solution space where FP BNs are clustered together. Given the size of the PKN and the small number of perturbations in the experimental data, the solution space can be very large containing billions of BNs making it difficult to enumerate true positive (TP)

BNs in reasonable time if it gets stuck in a cluster of FP BNs. We extend caspo-ts with a new enumeration scheme for breaking up clusters of similar solutions. In the following, we refer to the modified caspo-ts as caspo-ts D , where "D" is used to represent diversification. We apply both systems to two datasets: (1) an artificial dataset for modeling TCR signaling networks, and (2) the HPN-DREAM challenge dataset. Our results show substantial improvements of caspo-ts D in solution quality by discovering more signaling behaviors than caspo-ts. Moreover, caspo-ts D is able to find solutions in cases where caspo-ts is unable to find any in a reasonable time period. We also present an algorithm to improve the computational time of the model checking process. The results suggest huge computational time improvements in the model checking process.

Materials and methods

In this section, we describe the datasets, the caspo-ts system, the new algorithm to enumerate diverse BNs implemented in caspo-ts D , and the new algorithm to improve the model checking process. Note that we are using the phosphoproteomics time series dataset which has been described in details in Chapter 4, Section 4.2.2. We are using two types of datasets: an artificial (TCR Signaling) and a real dataset (HPN-DREAM).

Artificial dataset

The PKN was derived from the TCR signaling model of [KSRL + 06] and consists of 16 nodes and 25 edges. The artificial dataset for TCR signaling was generated in [OPS + 16] by simulating the PKN using logic based ODEs. This dataset consists of 4 readouts, 3 stimuli and 2 inhibitors. The readout proteins were measured at 16 time points under 10 perturbations.

HPN-DREAM

The HPN-DREAM dataset is described in the Chapter 4, Section 4.5. The PKN (associated to HPN-DREAM) consists of 64 nodes (7 stimuli, 3 inhibitors, and 23 readouts) and 178 edges (see Figure 4.6 described in detail in Chapter 4, Section 4.5.3).

The caspo-ts system

The caspo-ts system is based on a combination of ASP and model checking. The ASP part of the caspo-ts system is used to solve the combinatorial optimization problem of finding BNs compatible with a PKN and time series data. All learned BNs are optimized using an objective function, minimizing the distance between the original and the time series data determined by the BN learned with caspo-ts. The ASP solver guarantees finding all optimal solutions w.r.

Optimization

The caspo-ts method optimizes the distance of over-approximated traces of BNs w.r.t. the time series traces of the phosphoproteomic dataset. For this purpose, it uses Equation 4.1 described in detail in the Chapter 4, Section 4.2.4. From this optimization step, caspo-ts finds a bound, which is an optimal distance that can be achieved from traces of BNs to the traces of experimental data. Then, a constraint is added that only admits solutions that satisfy this bound. We can refer to this solution space as optimal solution space.

Then, different parameters of the ASP solver are configured to generate subset-optimal solutions. This means that enumerating subset-optimal solutions is a projection operation over the optimal solution space. We refer to this solution space as the subset-optimal solution space.

Notice that the optimal solutions space is generated by writing a logic program, while the subset-optimal solution space is the projection of optimal solution space by setting different parameters of the ASP solver.

Subset minimal BNs. The BNs learned by caspo-ts are represented by Boolean formulas in Disjunctive Normal Form (DNF), i.e., as a disjunction of conjunctive clauses 1 . The BNs inferred by caspo-ts use the smallest DNF formulas possible, in the sense that no conjunctive clause can be removed from a DNF formula without changing the Boolean function it represents. We refer to these BNs as subset minimal BNs.

Sampling of solution space of BNs. Sampling techniques are used to enumerate diverse solutions from the solution space [START_REF] Eiter | Finding similar or diverse solutions in answer set programming[END_REF]. Existing sampling techniques [START_REF] Zhu | Identification of potential drug targets in cancer signaling pathways using stochastic logical models[END_REF][START_REF] Romero | Computing diverse optimal stable models[END_REF] can be applied to the optimal solution space. In our case, we have to sample the subset-optimal solution space of the caspo-ts system. We discovered that existing sampling techniques can only be used with standard ASP solving, which means for the optimal solution space but not for the subset-optimal solution space. This is the main technical difficulty, which hindered our ability to enumerate diverse subset minimal solution for the caspo-ts system, as sampling 1. A clause can be seen as a reaction, where the proteins represented positively are active, and the proteins represented negatively are absent. A Boolean formula in DNF encompasses all possible reactions to update the value of a protein.

cannot be implemented on top of subset minimization in the current implementation of the ASP solver.

To resolve this issue, we enumerate subset minimal solutions using heuristics instead of configuring the ASP solver's parameters. Finally, a heuristic based approach is used to enumerate diverse optimal subset minimal solutions. This process of enumeration is explained in Algorithm 4.

In the following, we give a brief description of ASP solving algorithms used by caspo-ts.

But first let us have a look at an example Boolean formula.

ASP solving

Backtracking algorithms are used to solve computational problems by finding (all or some) candidate solutions to a problem. Davis Putnam Logemann Loveland (DPLL) [START_REF] Davis | A computing procedure for quantification theory[END_REF]DLL62] is a backtracking algorithm which is used to solve the Boolean satisfiability problem (SAT) [START_REF] Biere | Handbook of satisfiability[END_REF]. Given a Boolean formula, the goal is to find an interpretation assigning truth values to literals satisfying the formula. DPLL starts with arbitrarily choosing a literal and then assigning a truth value to it. If this assignment does not lead to a conflict then chose another literal and assign a truth value to it. In case of conflict, it backtracks to the last choice. ASP solvers use the conflict driven clause learning algorithm (CDCL), which is a variant of the DPLL backtracking algorithm (Figure 5.1). The main difference difference between CDCL and DPLL lies in the backtracking step. In the CDCL algorithm, when a conflict arises then instead of going just one step back it tries to analyze it and go back to the origin of the conflict (as shown in Figure 5.1).

Next, we explain in detail how the ASP solver clingo used by caspo-ts discovers solutions (BNs) using the CDCL algorithm [START_REF] Gebser | Clingo= asp+ control: Preliminary report[END_REF], shown in Algorithm 3.

Example 5.2.2. Before continuing any further, we build up our running example. Our running The algorithm works by extending a Boolean assignment over the atoms occurring in the given logic program P until a solution is found. The assignment is initialized in line 1. Then it is extended by the decision heuristic and propagation in the loop in lines 2-10. The call to Decide() in line 3 at the beginning of the loop uses a heuristic to select an atom, makes it either true or false, and adds it to the assignment. The consequences of this decision are then propagated in the following line extending the assignment accordingly. Then it is checked if propagation leads to a conflict. If this is the case, the conflict is analyzed in line 6 and the assignment adjusted in line 8 accordingly. Note that a call to Backjump() takes back one or more decisions together with their consequences, and adds an additional consequence to the assignment. This property ensures that the algorithm always terminates. It can also happen that a conflict cannot be recovered from. In this case, the problem is found unsatisfiable and the algorithm returns in line 10. Once the assignment is complete, the corresponding solution (set of true atoms) is returned in line 11.

Example 5.2.3. We can now apply Algorithm (3) to our running example (c.f. Example 5.2.2).

Starting with an empty assignment, we set a to false as the first decision. There are no immediate consequences and, hence, no conflict can arise. Then we decide to make b false. The consequences of this decision are that c is false via rule (5.3), and c is true via rule (5.4). Hence, we get a conflict, which is resolved and followed by a backjump. Since the conflict was caused by deciding a truth value for b, but is independent of the decision for a, the algorithm takes back all decisions and adds b as a consequence (we now know it must be true in all solutions).

We can then decide to make a false again, which sets c to true via rule (5.2). This decision does not cause a conflict and the assignment is no longer partial. Hence, the algorithm terminates with solution {b, c}.

Improvements in caspo-ts

Modifications in the ASP solver: caspo-ts D

Here, we describe the algorithm used for enumerating diverse subset minimal solutions.

In caspo-ts, the algorithm is implemented in the Python programming language using clingo's multi-shot solving API [START_REF] Kaminski | A tutorial on hybrid answer set solving with clingo[END_REF]. The API allows us to customize the solving process, in particular, it allows us to customize the decision heuristic of the solving component, which is

the key feature to find subset minimal answer sets. Note that we implemented the algorithm using the multi-shot solving API of clingo version 5. For that, we have upgraded the solver clingo of caspo-ts from 4.5.4 to 5.

Input: program P and atoms T to subset minimize 1 Prepare(P); Algorithm 4 is used to enumerate diverse optimal subset minimal answer sets. The idea is to configure the decision heuristic in lines 8 and 10, so that it first makes all atoms subject to subset minimization false before deciding truth values for other atoms. This decision heuristic is used in function Decide() in Algorithm 3 line 3. This modification ensures that the solution obtained from Algorithm 3 by calling Solve() from Algorithm 4 is a subset minimal solution (see [START_REF] Brewka | Implementing preferences with asprin[END_REF] for more details). Such a solution is output in line 12 of the algorithm. Furthermore, the algorithm calls Solve() in line 5 multiple times to find all subset minimal solutions. To not enumerate solutions twice, a constraint preventing to find the same solution again is added to the logic program P in the following line. This constraint is violated whenever a superset of the atoms in the previously found solution is true. This process is repeated in the loop in lines 4-12 until the program is no longer satisfiable and, hence, all solutions have been enumerated. So far we only discussed how to enumerate subset minimal solutions. Now we explain how to extend the method in order to compute diverse solutions. The key idea is to make the next solution different from the previous one. For this purpose, we assign atoms appearing in the last solution to false before assigning any other atoms. To modify the heuristic, we use function SetSign(a, t, l), which instructs the decision heuristic to assign atom a to truth value t on level l. The decision heuristic assigns free atoms with the highest level to the designated truth values before assigning atoms on lower levels. By default all atoms have level 0 and the decision heuristic is free to make them either true or false. The loop in lines 7-8 instructs the decision heuristic to assign atoms that appeared in the last solution to false on level 2. Any other atoms subject to subset minimization are assigned to false on level 1 in lines 9-10. We see in the experiments in the next section that this strategy breaks up clusters of similar solutions in the solution sequence.

Example 5.2.4. We continue with Example 5.2.2. Let us assume that a, b, and c are the atoms subject to subset minimization. Note that in Example 5.2.3 all decisions assigned atoms to false, so the first solution {b, c} obtained is in fact a subset minimal solution. Let us further assume that Algorithm 4 produced this solution in the first iteration (line 5). First, the constraint ← b∧c preventing any superset of {b, c} as solution is added in line 6. Then, the decision heuristic is configured to set atoms b and c to false on level 2 (line 8) and atom a to false on level 1 (line 10). In the next iteration, the only possible decision is to set c to false because b is already irrevocably assigned. The consequence of this decision is to set a to true via rule (5.2). Hence, we obtain solution {a, b}, which is a subset minimal solution. This is followed by adding the constraint ← a ∧ b in line 6, which in turn makes the program unsatisfiable and causes the loop in lines 4-12 to terminate. We correctly obtain the subset minimal solutions {b, c} and {a, b}.

Solution {a, b, c} is not subset minimal and not enumerated.

Proofs for Subset Minimal Solution Enumeration

Lemma 1. In line 4, Algorithm 4 computes a subset minimal w.r.t. a set of atoms S and the currently active program.

Proof. We only change the order and truth value with which the decision heuristic of the solving algorithm from [START_REF] Gebser | Conflict-driven answer set solving: From theory to practice[END_REF] assigns atoms. With this modification, the algorithm is still guaranteed to find a solution if there is any. All that remains to show is that the solution found is subset minimal.

When clingo finds a stable model of a program P , it has successively assigned truth values from {true, false} to all atoms. This sequence of assignments has the following property: (P)

Each assignment is either a decision or a result of propagation. Let a be any atom that is the result of propagation. Then the assignment of a is caused by the set of decisions D assigned before a. Every stable model of P that satisfies the literals in D also assigns a to the same truth value. Now consider that we want to subset minimize the set of atoms S and we use the heuristic method as in Algorithm 4 to obtain a stable model M . As long as there are unassigned atoms in S, the decision heuristic assigns them to false, so every atom a ∈ S that is true in M is a result of propagation, and all decisions before a are atoms from S assigned to false. Hence, using property (P), we have that every stable model that satisfies the atoms from S assigned to false in M also satisfies the atoms from S assigned to true. Then, there can be no stable model that is a subset of M w.r.t. S.

Lemma 2. Let P be a logic program and M be a subset minimal stable model of P w.r.t. a set of atoms S. Let the intersection of M and S be the set {a 0 , ..., a n }, and let C be the constraint

← a 0 ∧ • • • ∧ a n .
Then, the subset minimal stable models of P are exactly M and the subset minimal stable models of P ∪ {C}.

Proof. The constraint C eliminates all stable models M with {a 0 , . . . , a n } ⊆ M , so C eliminates M if (M ∩ S) ⊆ (M ∩ S). This means, the constraint eliminates all stable models that are a superset of M w.r.t. S, and the result follows.

Theorem 1. Algorithm 4 calculates all subset minimal stable models w.r.t. a set of atoms S.

Proof. By Lemma 1, the first stable model calculated by the algorithm is a subset minimal model. By Lemma 2, also the follow up stable models are subset minimal. To see that all subset minimal stable models are computed, note that by Lemma 2 in each iteration exactly one subset solution is eliminated and that the solving algorithm guarantees to find a solution if there is any.

Modifications in the model checking

Here, we describe the modifications to reduce the runtime of the model checking process.

Recall that a BN is given to a model checker to verify a set of traces under a set of perturbations (P) existing in the experimental data. These traces are specified using a combination of existential (E) and eventually (F ) operators (as explained in Chapter 4, Section 4.2.5).

The model checking component of the caspo-ts system takes a BN and verifies these traces.

Model checking itself is a PSPACE-hard problem. Giving traces as one big specification to the model checker can result in large overhead in computational time. Because in practice, a model checker performs better when solving independent problem separately. In order to resolve this issue, we propose to split the big specification in smaller parts and verify the traces in a parallel way where as soon as the traces under one perturbation are not reachable, we stop the model checking process and reject the BN as a false positive. The following Algorithm 5 explains the proposed model checking process. associated with the same node are grouped together by shading the background alternatingly in light gray and white. The frequencies of the clauses are depicted on the y axis. The red line depicts the frequency considering all 68338 solutions, while the green and blue lines depict the frequencies of the first 100 solutions computed by caspo-ts and caspo-ts D , respectively. In total, there are 49 clauses appearing in the family of BNs. While caspo-ts D (blue line) discovered 48 clauses, caspo-ts (green line) learned only 29 clauses by enumerating the same number of BNs (100). We also observe that the blue line is often much closer to the red line (with an average distance of 0.06) than the green line is (with an average distance of 0.20). This shows that the diverse enumeration scheme is able to produce solutions that are less similar to each other and better represent the solution space. The underlying ASP solver of caspo-ts by default uses an enumeration scheme that backtracks to the first point from which the next solution can be found. This approach typically leads to the situation where successively enumerated solutions only change in a small part. We can observe this behavior in Figure 5.3, where some clauses are overrepresented. axis stands for a block of 1000 solutions. The y axis depicts the percentage of true positives in a block of solutions. The red line depicts the overall true positive rate (78%), while the green and blue lines depict the true positive rates of caspo-ts and caspo-ts D , respectively 2 . We observe that for the caspo-ts system there are a lot of blocks with either a lot of true positives or very few. This suggests that true positives are clustered in the sequence of enumerated solutions. We observe that the diverse enumeration scheme does not show this behavior. This is especially 2. For example, when x has value 3000, the y value in blue gives the true positive rate among the solutions 2001 to 3000 computed by caspo-ts D . important for enumerating true positive solutions of real world instances where only a limited number of solutions can be checked because of time constraints. With the original caspots system, it can happen that the first cluster does not contain any true positives, making it impossible to find any true positive solution within a given time budget. The graph also shows that the diverse enumeration scheme does not sample over the full solution space. We see that before around 23000 solutions, the true positive rate is below the ideal 78% and then jumps up afterward. Thus, we conjecture that our enumeration scheme mainly breaks up local clusters of solutions with the artificial benchmark. Still, it is able to discover almost all clauses compared to the whole solution set (the frequency is 0 only once), while caspo-ts does not discover 20 clauses at all.

HPN-DREAM challenge dataset

Next, we show the results of applying diverse solution enumeration to the HPN-DREAM challenge dataset [HHC + 16]. We discuss the results according to three aspects: (1) time to compute the first true positive BN, (2) similarity among the family of solutions, and (3) Boolean functions computed by the original caspo-ts and the extended caspo-ts D system. We start the analysis with four cell lines, and then we provide a detailed analysis of the Boolean functions of one cell line discovered by caspo-ts and caspo-ts D .

Given that model checking is a computationally hard problem, we stop an experiment after a system verifies (using the model checker) 46 BNs per cell line3 . The model checking task was performed on a server with 1.5 Tb of RAM. Table 5.1 shows the number of TP BNs obtained for each cell line. We see that the number of TP BNs differs comparing caspo-ts and caspo-ts D .

For MCF7 we obtain 0 TP BNs with caspo-ts and 4 TP BNs with caspo-ts D , while for BT549 we obtain 2 and 14 TP BNs, respectively. For the other two cell lines BT20 and UACC812, we obtain a comparable number of BNs. We observe that we can get more TP solutions by checking the same number of BNs with caspo-ts D . This is an important improvement given the fact that model checking BNs is a computationally hard problem. Next, we consider the time column showing the time to compute the first TP BN for each cell line. We see that we are unable to get TP BN with caspo-ts in case of the MCF7 cell line, which shows that the caspo-ts system is stuck in a part of the search space where there are only FPs. Otherwise, for the other cell lines the time to get the first TP BN is comparable. We conclude that the difficulty to model check a BN depends on the cell line and not on the order in which solutions are found. Finally, the similarity column shows the similarity score among the set of TP BNs for each cell line. This score is calculated by comparing the clauses of one cell line with each other. We observe that the similarity among the solutions is much higher for caspo-ts than for caspo-ts D . From this we conclude that, studying the same number of networks, caspo-ts D can discover more clauses representing diverse signaling behaviors. Boolean functions with clauses of size one. We also noticed that the learned BNs are very similar to each other, as we see in Table 5.1 with the similarity score of 0.81. This relates to the fact that the ASP solver used by caspo-ts uses a backtracking algorithm to enumerate solutions and, hence, the solutions only change in small parts.

Next, we analyze the UACC812 cell line using caspo-ts D . Figure 5.6 shows the union of 10 TP BNs obtained by caspo-ts D . Nodes, edges and colors have the same meaning as in are represented by black rectangles where the nodes of incoming edges are their elements.

Additionally, we use dashed edges to represent clauses that were also discovered by caspo-ts.

In total, caspo-ts D discovered 66 clauses, 41 more than caspo-ts. It identified 23 out of the 25 clauses discovered by the caspo-ts system, and 43 additional clauses, studying the same number of BNs. It is important to note that even though for the UACC812 cell line caspo-ts D learned 5 TP BNs less than caspo-ts, the number of clauses learned by caspo-ts D is 3 times higher. It is interesting to see that 4 new stimuli appear in the solution set using caspo-ts D while there is only one stimulus in the previous solution set discovered by caspo-ts. We also discover one additional readout "JNK_pT183_pT185" using caspo-ts D as compare to caspo-ts. Since we find much more clauses with caspo-ts D , we can get an impression of the whole solution space by just inspecting a limited number of solutions. This analysis shows the efficacy of the extended caspo-ts D system in a real case scenario, where it is difficult to study the complete solution space because of time constraints.

Computation of root mean square error

Recall that the UACC812 cell line was the most difficult one to obtain true positive BNs (see Section 4.5.4). We discussed previously that it was because of the enumeration scheme used by the ASP solver. In this chapter, we have resolved this issue and obtained true positive BNs for this cell line. In Chapter 4 and Section 4. The previous model checking process could not finish verifying all 231 BNs for 7 out of 16 samples while parallel model checking finished verifying all BNs for all 16 samples. The total time spent on verifying 16 samples using the previous model checking approach is more than 1016 days while the total time spent on verifying reachability using parallel model checking is 148 days. This analysis shows the efficiency of the presented algorithm, where it is extremely time consuming (in some cases more than 3 months) to verify the BNs using the previous implementation.

Discussion

We have presented an algorithm to enumerate diverse optimal solutions with the caspo-ts system. The new algorithm extends the approach of [START_REF] Romero | Computing diverse optimal stable models[END_REF] for computing diverse optimal solutions. The novelty of this extension is that by modifying the heuristic of the solver we manage to enumerate solutions that are both optimal and diverse. There are other approaches for computing diverse solutions [EEEF09, HHOW05b, Nad11b] but they do not consider optimization problems.

Our key findings suggest that we retrieved a more complete set of mechanisms explaining the experimental data and better approximate biological reality, by sampling the large solution space of BNs. We were able to discover almost all behaviors (48 out of 49) existing in the and therefore can be applied to any problem modeled in ASP.

We have also improved the issue of the huge amount of time required for categorizing the BN as TP or FP using model checking. An algorithm is presented to check reachability in a parallel fashion. The substantial improvements have been shown in computation time by demonstrating the huge time difference (868 days) between the previous and the new implementation of the model checking process.

Conclusion

In this chapter, we have presented an extended version of the caspo-ts system to enumerate diverse solutions. We applied caspo-ts D on an artificial dataset (TCR signaling) as well as a real case study (HPN-DREAM) to learn diverse BNs. We compared the results with the caspots system, showing a substantial improvement in solution quality. For one, we discovered more signaling behaviors (clauses) comparing solutions enumerated with both systems. For another,

Conclusions and Future Work

In this chapter, we summarize the work and outline future perspectives.

Summary

Protein signaling networks are static views of dynamic processes where proteins go through many biochemical modifications such as ubiquitination and phosphorylation to propagate signals that regulate cells and can act as feedback systems. Understanding the precise mechanisms underlying protein interactions can elucidate how signaling and cell cycle progression occur within cells in different diseases such as cancer. In this thesis, we have advanced our understanding in breast cancer by using the caspo-ts method to learn BNs from multiple perturbation phosphoproteomic time series data given a Prior Knowledge Network. We have improved and adapted caspo-ts to deal with a midscale PKN with 64 nodes and 178 edges in order to learn the BNs of four breast cancer cell lines (BT20, BT549, MCF7, UACC812) from their time series phosphoproteomic datasets. Importantly, the PKN did not contain any information about the temporal changes or dynamic properties of the proteins. This information was learned from a dataset describing the dynamics of signaling processes for those breast cancer cell lines as part of the HPN-DREAM challenge. The initial results highlighted several characteristics of this method when applied on HPN-DREAM dataset (a real case study):

1. high false positive rate, 2. large solution space (billions of BNs), 3. high computation time of caspo-ts-model checking component.

Future perspectives 6.2.1 Simulation of Boolean models

It has been shown that the HER2 amplified breast cancer demonstrates significant increase in BCL2 and ER expression when treated with anti HER2 therapy. The co-regulation of ER or BCL2 with anti HER2 drug is shown to be useful for decreasing the BCL2 and ER overexpression in HER2 positive UACC812 cell line. These molecular entities (BCL2 and HER2/3) exist in our Boolean model of UACC812 cell line. In the near future, we plan to execute the Boolean model of UACC812 cell line derived using caspo-ts D presented in Chapter 4 (Figure 5.6), to study the impact of inhibiting these molecules in our Boolean models and compare the results with literature. We believe this study will be useful for final validation of our models.

We studied some software to simulate boolean models but we faced a lot of trouble configuring different parameters for the toy example. Different software produces different simulation results. We also had problem reproducing the results published with these software. Therefore, we plan to develop our own tool to simulate the Boolean models and embed it in the caspo-ts D system .

Extension of diversity algorithm

For the first time, with our study, we observe that solutions are not randomly distributed along the search tree of caspo-ts, as shown in Chapter 5. In the near future, we plan to extend the diversity algorithm in two directions. First, we are planning to experiment with solver parameters in order to introduce some randomness into the search. Second, we intend to extend the algorithm to call the model-checker only on answer sets which are diverse (according to some measure to be defined). This is possible because the time to enumerate over-approximated solutions using the ASP solver is much lower than the time needed to check solutions using the model-checker. We expect both enhancements to further improve the diversity of the discovered solutions.

Extension of model checking algorithm

In the current implementation, we parallelize the properties on the perturbation level. It can be improved further by parallelizing the traces within a perturbation as well. This is possible because we can prove in a parallel fashion that one time point is reachable from the preceding time point, too. Typically, it can be quickly determined if a property does not hold but it might take a long time to verify that property holds. A BN is a false positive if one of the reachability property does not hold, and a lot of computation time might be spent before detecting a false positive BN. However by parallelizing verification at time point level, we can stop the model checker from verifying reachability of other time points and other perturbations as soon as one transition cannot be verified. We believe that this improvement will help in gaining computation time particularly in case of many time points.

Experimental design

Currently, our research is in the context of identifying BNs representing a part of a signaling pathway from real time series of phosphoproteomics data using caspo time series method which is based on answer set programming and model checking. Phosphoproteomics data contain a lot of noise, and also our approach is exhaustive, so it results in learning a family of Boolean Networks. This highlights the need for designing efficient experimental strategy. To deal with this issue, we proposed to: (1) identify a way to characterize different dynamical behaviors among a family of Boolean networks and (2) develop an algorithm based on the maximum entropy measure to select the efficient experimental perturbation to make this family of dynamic behaviors less variable.

7.5 Supervision We have used the Epidermal Growth Factor (EGF) and Tumor Necrosis Factor Alpha (TNFα)

as benchmark [MTH + 12]. Our results suggest that the initial condition affect the evolution of the dynamics of the Boolean model. Abstract: Traditional canonical signaling pathways help to understand overall signaling processes inside the cell. Large scale phosphoproteomic data provide insight into alterations among different proteins under different experimental settings. Our goal is to combine the traditional signaling networks with complex phosphoproteomic time-series data in order to unravel cell specific signaling networks. On the application side, we apply and improve a caspo time series method conceived to integrate time series phosphoproteomic data into protein signaling networks. We use a large-scale real case study from the HPN-DREAM Breast Cancer challenge.

We infer a family of Boolean models from multiple perturbation time series data of four breast cancer cell lines given a prior protein signaling network. The obtained results are comparable to the top performing teams of the HPN-DREAM challenge. We also discovered that the similar models are clustered together in the solutions space. On the computational side, we improved the method to discover diverse solutions and improve the computational time.

Figure 1 . 1 -

 11 Figure 1.1 -Central dogma of biology (Image source: Genome Research Limited).

Figure 1 . 2 -

 12 Figure 1.2 -Caspo-ts workflow. Prior Knowledge Networks (PKNs) are extracted from literature curated databases containing information about interactions between different proteins or genes. PKNs are available in different databases such as Reactome, PID, etc. Phosphoproteomic time-series data show the measurement of different proteins at different time points under multiple perturbations. A BN consists of a set of nodes where a Boolean function is assigned to each node. The state of each node is updated by evaluating the Boolean function. Given phosphoproteomic time series data we construct a PKN by querying pathway databases. After normalizing the time series data, we use it together with the PKN as input of caspo-ts (ASP component) for learning BNs. Finally, caspo-ts, uses a model checking step to filter false positive BNs. In this figure, the two main components of caspo-ts are shown in orange.

Figure 2 . 1 -

 21 Figure 2.1 -Breast structure (Image source: National Cancer Institute www.cancer.gov).

  10, FAC13]. They became a powerful tool to investigate how apoptosis, migration and proliferation are deregulated during breast cancer progression. Breast cancer cell lines help to generate quantifiable and reproducible results [VGR07]. They have provided valuable insights into understanding different aspects of breast cancer [LL04]. They have created a large amount of knowledge about breast cancer biology [Eth96, SLB + 10, YXW + 10, TCJ + 10]. They have suggested novel cancer therapies over the past decades [WBB + 99, WBM + 98, Mas00]. They have been used in preclinical studies and predicted accurate clinical outcomes [CGL10]. In this thesis, we are using four breast cancer cell lines (BT20, BT549, MCF7, and UACC812) to study signaling behaviors.

  mutants. BRCA1 germ mutations put female mutation carriers to high risk of ovarian and breast cancer. Four new cell lines were identified with BRCA1 mutations [EHN + 06]. This cell line belongs to the triple negative breast cancer, which means hormone therapy cannot be applied to this kind of cancer. Triple negative breast cancer patients have a poor treatment outcome because of the lack of validated molecular targets. There is a clear need to enhance knowledge about this type of cancer to develop better therapies. The BT20 cell line has been used to study the triple negative breast cancer type [CGL10]. Results suggests that the BT20 cell line has very high expression of the EGFR protein and a genetic amplification of the EGFR gene [LG87].

  The single cell cloning of cells without classification of ER positive or negative clones, highlighted the heterogeneity existing in breast cancer [OZG + 01]. MCF7 cell lines do not have an amplified level of HER2. However, Osborne et al. have developed MCF7 cells with amplified HER2 [BSS + 92]. They showed how anti HER2 inhibitors can block growth and can be used for treatment. The studies have been validated clinically. Moreover, many laboratories have developed anti HER2 resistant cell lines. Despite the limitation imposed by cell culture, this cell line has tremendously advanced the knowledge about breast cancer and gave new directions to breast cancer research [LOD15]. BT549 In 1978, the BT549 cell line was derived from a 72 year old women. This cell line was established by Coutinho and Lasfargues. It is a breast ductal carcinoma cell line. This cell line represents the triple negative breast cancer type [LC81]. Lehmann et al. analyzed gene expression profiles from 21 breast cancer datasets. They identified 587 triple negative breast cancer cases from these breast cancer datasets. Cluster analysis was performed to identify 6 subtypes of triple negative breast cancer. They derived gene expression signatures from these subtypes to select the representative cell line of these subtypes.

  cell line was established by Meltzer and colleagues [MLD + 91]. It represents the HER2 positive breast cancer type. Wang et al. investigated mechanisms of resistance to HER2 targeted drugs using a panel of HER2 positive cell lines. They used two drugs (trastuzumab and lapatinib) alone and in combination to study drug resistance. Their results revealed that the resistance to trastuzumab is due to the reactivation of the HER2 pathway. Furthermore, resistance to lapatinib or lapatinib combined with trastuzumab is associated with an alternative signaling through the ER pathway.

  .2). Hence, RPPA is ideally suited for measuring the level of a single protein for multiple samples simultaneously [NBK12, MLE10]. For example, it can be used to study the regulation of proteins in healthy and cancerous cells [HRS + 09]. It is mainly used to: (1) define new and efficient targeted drug therapies for individual patients, (2) identify and validate biomarkers, and (3) verify drug effects (on or off targets, downstream signaling) [SGTP17].

Figure 2 . 2 -

 22 Figure 2.2 -Forward and reverse phase protein array format.

Figure 2 . 4 -

 24 Figure 2.4 -Work-flow of the Computational Modeling (This work-flow was inspired by [KHJ + 06]).

Figure 2 .

 2 5 shows general steps involved in solving problems using genetic algorithm. It starts with defining an initial population. Each individual in this population represents a solution to the problem. A fitness function is evaluated for each individual according to some objective function. Then a stopping criteria is evaluated to create a new population or to terminate the search process. Genetic algorithms have been widely used for solving optimization problems. They are also used for parameter estimation of mathematical models based on ODEs.

Figure 2 . 5 -

 25 Figure 2.5 -Genetic algorithm.

2 )

 2 was inferred by fitting the simulated data to the cancer dataset using the genetic algorithm. This dataset contained six phenotypic responses of 35 proteins at five time points , for three triple negative breast cancer cell lines BT20, MCF7 and MDA-MB-453. The genetic algorithm was iterated 150 times to discover 50 cancer networks. The rewired events discovered by the genetic algorithm were different in different cancer networks. The rewired events with high frequency were kept to construct N 2 network. The N 2 network is a modified version of the N 1 network, containing two newly added edges while 3 edges were deleted. Various inferred rewired events were verified and confirmed through existing literature. From the N 2 network, the cancer network (N 3 ) is inferred by fitting the simulated data to the drug treatment dataset using the genetic algorithm. Th drug treatment dataset consisted of measurements of the signaling and cell fate data under six treatments. The genetic algorithm was repeated 150 times to discover the drug sensitive network. The N 3 network is a modified version of the N 2 network, containing two newly added edges while 2 edges were deleted. Most of the rewiring events inferred by this method are consistent with the literature. In this work, authors modeled only 18 out of 35 proteins existing in the experimental dataset.ConclusionHere, I have described two studies based on non-deterministic stochastic search method (genetic algorithms). Saez et al. [SRAE+ 09] used genetic search algorithm to train Boolean networks against experimental data. However, it was not guaranteed that genetic search would yield the lowest value of the objective function, so all interactions were exhaustively checked to decrease the model size. Mishra et al. [MBC+ 15] also used the genetic algorithm to infer models which fit experimental data. They discovered different rewiring events over different iterations of the algorithm. Hence a decent number of iterations were required to converge the solutions. Overall, stochastic search methods are useful for modeling signaling networks,

  The prior distribution contained the expected number of edges (interactions existing in the PKN) and penalized the unusual edges (not part of expected edges). Finally, the posterior probabilities of edges are calculated by model averaging. Model averaging is used to highlight edges which are common in many network structures. This network building approach was applied to learn signaling network for the MDA-MB-468 breast cancer cell line. The learning data contained measurements of 20 proteins at 8 time points under 4 growth conditions. The PKN consisted of 20 protein and 74 edges. The robustness of the cell line specific network to the PKN was investigated by adding and deleting an edge in the PKN. The 25 perturbed PKNs were generated with changes in one third (25) of the total edges[START_REF]OSM +[END_REF]. The results showed the robustness to modification in the PKN. The robustness of the results to the data perturbation was also investigated, by deleting part of the data and replacing it with the average of adjacent time points. The analysis showed that results are robust to the perturbation of the data. They generated a set of testable hypothesis, discovered various known and unknown signaling behaviors. The known inferred edges were validated through existing literature knowledge. Some of the novel edges were validated through experimental validation.Hill et al. 2017Hill et al. [HNJC+ 17] learned cell line specific networks given PKN and phosphoproteomic data for four breast cancer cell lines (BT20, BT549, MCF7, UACC812) using a variant of dynamic bayesian networks. The learning data contained measurements of 35 proteins at 7 time points in all cell lines under a combination of perturbations (8 stimuli and 5 inhibitors). The PKN contains information about the connection of these 35 proteins with 65 edges curated from literature. A dynamic bayesian network approach was used to learn 32 context specific networks (4 cell line * 8 stimuli) using learning dataset given PKN. Each network contained 49 edges on average, and 40% of these edges were not part of prior network. The learned net-

  + 09], presented a framework based on integer linear programming (ILP) to learn cell type specific pathways and to study the effects of different drugs on these pathways. Two types of datasets were constructed; one to learn the cell specific signaling network, and the other to learn signaling alterations when four drugs were induced. The first experimental dataset (D 1 ) consisted of 13 phosphoproteins measured under 55 perturbations (stimuli and inhibitors). The second experimental dataset (D 2 ) contained measurements of 13 phosphoproteins under 55 conditions (stimuli and drugs). The experimental data were normalized to convert values to a range of [0,1]. The prior network consisted of 74 proteins and 105 reactions. They built a cell type specific network (Boolean model) given the prior network and the experimental dataset (D 1 ). They described their problem in the form of equations and stated two objective functions. The first objective function is used to infer the network which best fits the experimental data. The second objective function is used to minimize the size of the network. Then the ILP solver is used to infer networks satisfying the constraints and optimizing the objective functions. The inferred cell line specific networks consisted of 49 proteins and 44 reactions. This network contained 5 stimuli and 13 phosphorylated proteins. They evaluated the effects of four drugs on this cell type specific network: Raf kinase inhibitor (Sorafenib), two potent EGFR kinase inhibitors (Erlotinib), the dual EGFR/ErbB-2 inhibitor (Lapatinib), and Gefitinib. They used similar learning process described above to learn 4 drug induced networks; using the cell specific network as the prior network and experimental dataset (D 2 ).

3. 1

 1 IntroductionAnswer Set Programming (ASP) is an emerging simple, yet powerful, framework for reasoning and knowledge representation. The simple modeling language and flexible reasoning modes of ASP have led to the use of ASP in numerous domains such as systems biology [BCT + 04, GSTV11, VGE + 12, OPS+ 15], a decision support system for NASA shuttle controllers [NBG + 01, BGN06], and product configuration [SN99]. ASP is a declarative problem solving approach, in which a problem is modeled rather than telling the computer how to solve it[START_REF] Gebser | Answer set solving in practice[END_REF]. Figure3.1 shows the work-flow of the ASP paradigm. Modeling aims at describing the problem in the form of a logic program using first order variables. Afterwards, solving generates stable models using searching algorithms to find a solution for the stated problem.These resulting solutions are called stable models or answer sets.

Figure 3 .

 3 Figure 3.1 -ASP workflow

  The program {switchOn} ← light ← switchOn has the same solutions as the program in Example 3.2.1 without using the extra atom switchOff to capture that the switch is off. Since the first choice rule can be used to prove any subset of {switchOn} there are two solutions: ∅ and {switchOn, light}.

Example 3 . 2 . 3 .

 323 Coming back to our light switch domain, we model our toy problem just using choice rules and integrity constraints. {switchOn; light} ← ← light ∧ ¬switchOn ← ¬light ∧ switchOn

Figure 3 . 3 -

 33 Figure 3.3 -Directed graph after coloring

  Videla et al. Method description. Videla et al. [VGE + 12] proposed an approach based on ASP to infer BNs from PKNs and experimental data. It takes three inputs a PKN, perturbations, and experimental observations. A PKN specifies the canonical information about connections among proteins existing in the experimental data. Their method can only handle PKNs without any feed-back loops. A perturbation consists of a combination of stimuli and inhibitors. Experimental observations consist of measurements of readouts proteins. They model the combinatorial optimization problem of finding BNs using facts, rules, and the objective function. The objective function takes two parameters into account: the fitness of the model to the data, and the model size. They derive the structure of the model (BN) from the PKN. Then, the values of readouts are calculated for each model by fixing the value of stimuli and inhibitors. The optimization is applied on these models to minimize the fitness score and the size of the model. Case study. Videla et al. illustrated their approach on mid and large scale datasets, compared results with the genetic algorithm optimization framework [SRAE + 09] (Section 2.4.2).

caspo-ts) 5

 5 to learn Boolean models from phosphoproteomic time series data given a PKN. A family of candidate BNs compatible with the PKN is exhaustively enumerated. Afterwards an over-approximation constraint is used to filter out BNs which cannot reproduce the phosphoproteomic dataset. An objective function is used to minimize the distance between the actual time series and the time series produced by the BNs. Because of the over-approximation criteria, some of the returned BNs may not reproduce the time-series traces. Such false positives 5. This approach is based on ASP and model checking as described in the Figure1.2 (formal description is given in Chapter 3). can be ruled out using model checking. The remaining BNs are true positive. True positive BNs are guaranteed to reproduce all traces of the phosphoproteomic time-series data.

  3 and C 2.3 can take hours because the search space is huge. Next, they used model checking to determine the true and false positive rate of inferred models. The true positive rate was high (≥ 78%) in all case studies except C 2.3 and C 3 . They also verified the sensitivity the caspo-ts to the modified and incomplete PKN for case study C 1 . For this, they generated different versions of PKNs by adding, removing and replacing edges randomly. Results depicted sufficient robustness to the changes in the PKN.

  the set of directed edges and σ ⊆ E × {+1, -1} the signs of edges. The set of nodes is denoted by V = S ∪ I ∪ R ∪ U where S are stimuli, I are inhibitors, R are readouts, and U are unobserved nodes. Stimuli, inhibitors, readouts, and unobserved nodes are encoded by different colors in the graphs. Stimuli are shown in green, inhibitors in red, readouts in blue, and unobserved nodes in white (Figure 4.1). Moreover, the subsets S, I, R, U are all pairwise disjoint except for I and R. Stimuli

Figure 4 .

 4 Figure 4.2 depicts an example of phosphoproteomic time series data. The phosphoproteomic dataset is the second input of the caspo-ts system.

Figure 4 . 2 -

 42 Figure 4.2 -Phosphoproteomic time series data.Here the values between zero and one of three proteins are shown in different colors. In this figure, we see the time series of one readout protein (blue) under a perturbation of one stimulus (green) and one inhibitor (red). Stimuli have value 1 and inhibitors have value 0 across all time points of an experimental perturbation. Readouts take continuous values in [0;1] after normalization. In some phosphoproteomic datasets an inhibitor can also act as a readout protein, which means that there are perturbations where it will be measured.

Figure 4 .

 4 3 describes the model checking process of the caspo-ts system. It takes the Boolean Model and properties as input and generates two types of output:(1) true, which means the property is satisfied, or (2) false, with a counter example which means the property is not satisfied.

Figure 4 . 3 -

 43 Figure 4.3 -Model checking process.

inferred 3 4

 34 BNs are able to recover the discrete trajectories without any error. If the model RMSE is greater than the discrete RMSE then the inferred BNs have some errors in the recoverability of the discrete time series data. The delta denotes the difference between discrete and model RMSE. The lower value of delta means that the BN can successfully reproduce the traces of the experimental data with minimum error possible. Data: BNs and Time Series Traces Result: RMSE scores 1 for each model i in B do 2 Generate over-approximated traces for i ;Calculate the distance between the over-approximated traces of i and the actual traces using Equation4.1;Minimize the distance between the over-approximated and the actual traces; 5 end Algorithm 1: Computing Root Mean Square Error.

Data: 3 Calculate the intersection of A with M ; 4 Apply 5 end Algorithm 2 :

 3452 Gold Standard Network (A) and Family of Boolean Networks (B) Result: Similarity index for each model in the family of Boolean networks 1 Initialize two structures A and B; 2 for each model M in B do Jaccard similarity coefficient (J(A, M )) ; Graph Similarity Algorithm

Figure 4 . 4 .

 44 The values of 4EBP1 are oscillating between 2 and 5 while the values of BAD are oscillating between 0.6 and 1.8.

Figure 4 . 4 -

 44 Figure 4.4 -Protein values over time. Here, the x-axis represents the time points and the y-axis represents the protein measurements.

Figure 4 . 5 -

 45 Figure 4.5 -BT20 cell line dataset. Cues are the perturbations representing combination of stimuli and inhibitors. The x-axis shows the time point while the y-axis shows the normalized protein measurements.

  7, 4.8, 4.9 and 4.10). It is worth noting that we generated true positive BNs for UACC812 cell line by allowing the model checker to run without bounding it to the 7 day time limit. This result pointed us to the fact that the similar BNs are clustered together in the solution space generated by caspo-ts. We analyze this aspect and present the detailed results of the UACC812 cell line in the chapter 5.An aggregated network was built (Figure 4.11) by combining the BN families (with 191, 21, 72, and 20 BNs for BT549, MCF7, BT20, and UACC812 cell lines respectively) obtained for the four cell lines by keeping the hyper-edges (Boolean functions) having a frequency higher than 0.3 within each BN family. The frequency is calculated by counting the number of common Boolean functions and dividing it by the total number of Boolean functions within the BN family of each cell line. This aggregated network contains 34 nodes and 74 Boolean functions involving 36 AND gates. As compared to the PKN (Figure 4.6), the inferred networks are highly specific to each cell line. In Figure 4.11, all cell lines share only 4% of Boolean functions which are shown in thick black colored edges. This shows that the inferred BNs of these four breast cancer cell lines are very diverse and different from each other. To measure cell line similarity, we calculated the similarity score by applying the Graph Similarity Measure on the family of BNs (with 191, 21, 72, and 20 BNs for BT549, MCF7, BT20, and UACC812 cell lines, respectively). This algorithm receives two parameters as input: (1) one gold standard BN and (2) a family of BNs. It outputs a score in [0; 1], measuring the average of the similarity scores between each BN in the family and the gold standard BN. In our case, the gold standard BN is the aggregation of one family of BNs. The similarity scores between all pairs of breast cancer cell lines are shown in

Figure 4 . 7 -Figure 4 . 8 -Figure 4 . 9 -

 474849 Figure 4.7 -Union of BNs of BT20. Here, we show the union of BNs for the cell line BT20. This network is generated by combining 72 true positive BNs. It contains 31 nodes and 41 boolean functions with 12 AND gates. There are 2 stimuli, 2 inhibitors and 21 readouts.

Figure 4 . 10 -Figure 4 . 11 -

 410411 Figure 4.10 -Union of BNs of UACC812. Here, we show the union of BNs for the cell line UACC812. This network is generated by combining 20 BNs. It contains 33 nodes and 54 boolean functions with 29 AND gates. There are 6 stimuli, 2 inhibitors and 18 readouts.

  among their family of BNs with a variable frequency range. This shows that these cell lines have different mechanisms agreeing with the results obtained through graph similarity measure given in Table4.11.

Figure 4 .

 4 Figure 4.13 shows the common Boolean functions along with their frequency in all BNs.

Figure 4 . 13 -

 413 Figure 4.13 -Common Boolean functions across all four cell lines. The Boolean functions are represented on the x-axis and the frequency of each Boolean function is shown on the y-axis.

  We calculated the RMSE score for each Boolean model in a family of BNs of each cell line. We show these results in Figure4.14. We observe that the all models of MCF7 and BT20 cell lines have the same RMSE w.r.t testing data. The family of BT549 cell line has an average RMSE score of 0.0116 and UACC812 cell line has an average RMSE score of 0.0119.

Figure 4 .

 4 Figure 4.14 -RMSE for a family of BNs w.r.t testing data, where x-axis shows the true positive Boolean models and y-axis shows the delta. The delta captures the difference between discrete and model RMSE (see Section 4.2.6). Recall that we obtained different number of true positive BNs in each cell line (see Section 4.5.4).

Figure 4 . 15 -

 415 Figure 4.15 -Performance assessment with learning, testing and random datasets. The x-axis shows the cell line and the y-axis shows the RMSE ratio (see Equation (4.5)) of the inferred BNs from the HPN-DREAM data for each cell line with respect to the three datasets. The three datasets are encoded by different color codes. The RMSE ratio with respect to the HPN-DREAM learning and testing datasets is shown in blue and green colors, respectively. The random dataset RMSE ratio distribution is shown as red boxplots.

Figure 4 .

 4 Figure 4.16 shows the Receiver Operating Characteristic (ROC) curve of each cell line. ForBNs of each cell line, TPR and FPR was calculated using Equation (4.6) and (4.7). BT549 cell line models are the most accurate, followed by MCF7 and BT20. We can observe the clear distinction between true positive and false positive BNs. The BNs inferred by caspo-ts have an average AUROC score of 0.77 which is comparable to the AUROC score of 0.78 of the top performing method of HPN-DREAM challenge. A number of assumptions made during the modeling phase may have influenced our ranking. First, since our method can pinpoint the noisy, incomplete and erroneous experiment, it allows us to use only the reliable experimental settings. Second, our method constrains the solution space to the proteins existing in the PKN, and so anything outside the prior knowledge cannot be found. From Figure4.16, we can see that the caspo-ts method shows promising results for the inferred true positive BNs.

Figure 4 .

 4 Figure 4.16 -ROC curve across all cell lines. The x-axis shows the false positive rate and the y-axis denotes the true positive rate. The average AUROC score is 0.77.

Example 5 . 2 . 1 .

 521 To use Boolean formulas we discretize the phosphoproteomic data: values greater or equal to 0.5 are set to 1, and others are set to 0. Let protein A be associated to Boolean variable A have the formula (B) ∨ (¬B ∧ C) containing the two conjunctive clauses (B) and (¬B ∧ C), where B and C are Boolean variables representing proteins B and C. This formula can be used to update the value of A. If the update is applied, A is set to 1 if either the value of B is 1, or the value of B is 0 and the value of C is 1. Otherwise, the value of A is set to 0.

Figure 5 . 1 - 4 )

 514 Figure 5.1 -Difference between CDCL and DPLL algorithms. The dashed arrows represent back-jumps. The red and green arrows denote false and true assignment.

Figure 5 . 2 -

 52 Figure 5.2 -Conflict driven clause learning for example 5.2.3. The red arrows denote false assignment. The dashed arrow represents the back-jump. The gray box shows the assignment of truth values to literals including assignments from decisions and consequences. The green box represents the solution obtained.

5 S 11 if S is a true positive then 12 OutputAlgorithm 4 :

 511124 2 foreach x ∈ T do 3 SetSign(x, false, 1) 4 while satisfiable do ← Solve(); 6 AddConstraint(← a 0 , . . . , a n for {a 0 , . . . , a n } = T ∩ S); 7 foreach x ∈ T ∩ S do 8 SetSign(x, false, 2); 9 foreach x ∈ T \ S do 10 SetSign(x, false, 1); Diverse subset minimal solution enumeration.

Figure 5 . 3 -

 53 Figure 5.3 -Frequency of clauses per node in all 68338 BNs (red line), and in the first 100 BNs enumerated by caspo-ts (green line) and caspo-ts D (blue line).

Figure 5 .

 5 Figure 5.4 depicts the true positive rate of blocks of successive solutions. Each tick on the x

Figure 5 . 4 -

 54 Figure 5.4 -True positive rate of BNs grouped in blocks of 1000 networks.

Figure 5 .5 shows the union of 10 TP

 510 BNs obtained by caspo-ts. There are four different kinds of nodes in the graph: (1) stimuli shown in green, (2) inhibitors shown in red, (3) readouts shown in blue, and (4) unobserved nodes shown in white. Note that blue nodes with red borders are readouts, which are also inhibitors. There are two different kinds of edges shown in red and green color. Green edges are used to show a positive influence (←), and red edges are used to show a negative influence ( ). We discovered 25 clauses with caspo-ts, and observed that the learned BNs only contain

Figure 5 . 5 .Figure 5 . 5 -

 5555 Figure 5.5. Unlike with caspo-ts, here we identified clauses with more than one element. They

Figure 5

 5 Figure 5.6 -caspo-ts D : 10 optimal TP BNs concatenated for cell line UACC812. The dashed edges are used to represent clauses that were also discovered by caspo-ts. AND gates are represented by black boxes.

5. 3 . 4 Figure 5 . 7 -

 3457 Figure 5.7 -RMSE score for true and false positive networks of UACC812 cell line.

7. 5 . 1

 51 Simulator for Boolean models I have supervised three students to develop a tool in python which can simulate and visualize boolean models. We have managed to develop a first version of the tool. This tool takes a Boolean model as input, simulates it under given experimental settings and visualize the results.

7. 6

 6 Collaborative Research 7.6.1 Interaction graph for dream 11 challenge I also participated in the Dream 11 challenge focusing on discovering molecular signatures to virus exposures [Sie16]. More precisely, I worked on generating PKNs for the Dream 11 challenge for merging knowledge driven methods with data driven methods. An influence graph is extracted from different literature databases based on the selection of core variables. I inferred six influence graphs from literature curated databases i.e BioCarta, Go Ontology, Kegg, NCI, Reactome, and Wiki Pathways using the enrichnet web tool. In order to cover all core variables existing in dataset, I merged influence graphs to generate a new enriched influence graph. According to the results obtained, there is still a lot of room for improvement in the data driven algorithm to improve overall performance. Titre : Intégration de données de séries temporelles phosphoprotéomiques dans des réseaux de connaissances antérieurs Mots clés : réseaux de signalisation, programmation logique, vérification de modèles, lignées cellulaires Résumé: Les voies de signalisation canoniques traditionnelles aident à comprendre l'ensemble des processus de signalisation à l'intérieur de la cellule. Les données phosphoprotéomiques à grande échelle donnent un aperçu des altérations entre différentes protéines dans différents contextes expérimentaux. Notre objectif est de combiner les réseaux de signalisation traditionnels avec des données de séries temporelles phosphoprotéomiques complexes afin de démêler les réseaux de signalisation spécifiques aux cellules. Côté application, nous appliquons et améliorons une méthode de séries temporelles caspo conçue pour intégrer des données phosphoprotéomiques de séries temporelles dans des réseaux de signalisation de protéines. Nous utilisons une étude de cas réel à grande échelle tirée du défi HPN-DREAM Breast Cancer. Nous déduisons une famille de modèles booléens à partir de données de séries temporelles de perturbations multiples de quatre lignées cellulaires de cancer du sein, compte tenu d'un réseau de signalisation protéique antérieur. Les résultats obtenus sont comparables aux équipes les plus performantes du challenge HPN-DREAM. Nous avons découvert que les modèles similaires sont regroupés dans l'espace de solutions. Du côté informatique, nous avons amélioré la méthode pour découvrir diverses solutions et améliorer le temps de calcul. Title : Integrating Phosphoproteomic Time Series Data into Prior Knowledge Networks Keywords : signaling networks, logic programming, model checking, cell lines

Table 2 .

 2 1 -Molecular classification of breast cancer cell lines.

	Cell Lines ER	PR HER2	Subtype
	BT20	-	-	-	Triple Negative A
	BT549	-	-	-	Triple Negative B
	MCF7	+	+	-	Luminal A
	UACC812 + + / -	+	Luminal B

  Giuliano et al. analyzed the effects of HER2 targeted therapies on ER and Bcl2 expression in clinical tumor samples and preclinical models. Results suggested that Bcl2 and ER expression increased significantly in breast tumor xenografts treated with anti HER2 therapies. They revealed that the co-regulation of Bcl2 or ER with an anti HER2 therapy can prevent the increased Bcl2 and ER expression in HER2 positive breast cancer patients. They also reported that tumor progression slowed down with endocrine therapy in the presence of restored ER expression in xenograft tumors treated with anti-HER2 therapy [GHW + 15]. Zhang et al. investigated the role of EPOR (erythropoietin receptor) in the inhibition and the resistance to the trastuzumab drug, in HER2 positive breast cancer. They identified EPOR

mRNA and protein expression in HER2 positive breast cancer cell lines (UACC812, MDA-MB-453, and SKBR3). They suggested that EPOR expression may influence tumor proliferation and progression in HER2 positive breast cancer [ZDX

+ 12]

.

Brennan et al. studied the connection between the JAM-A (junctional adhesion molecule

A) and aggressive tumor, using breast cancer cell lines and clinical datasets. They suggested that over-expression of JAM-A may increase breast cancer progression in HER2 positive breast cancer. They also suggested that JAM-A can serve as a potential drug target and biomarker of HER2 positive breast cancer [BMH

+ 13]

.

Overview of analysis of signaling networks using RPPA technology. RPPAs offer excellent data for studying deregulated signaling networks in cancer. It allows

		CHAPTER 2. BACKGROUND AND RELATED WORK
	Health in 2018 and has been given access to by some insurance companies. This assay can be
	used to guide personalized treatment by linking drug targets with available therapies [The18].
	In contrast to RPPA, there exist other proteomic technologies such as mass spectrometry,
	western blot, and enzyme linked immunosorbent assay (ELISA) (See Table 2.2). Mass spec-
	trometry based methods can analyze thousands of proteins at a time but they are unable to
	resolve all proteins within a sample. In contrast, RPPAs are highly sensitive and can detect
	even low abundance proteins. Mass spectrometry based approaches can be used to discover
	Cell Culture new candidate biomarkers by comparing healthy tissues with cancerous tissues. RPPA can Protein Lysates Arrays probed Arrays with Antibodies be used to validate these markers in small samples (taken from patients), thereby deriving an Slide 1 Slide 2 Antibody 1 Antibody 2 individual precise cancer therapy. There are 380 validated antibodies available for RPPA. As
	Slide 3 compared to ELISA which requires two antibodies for a protein, RPPA requires only one anti-Antibody 3
	Slide 4 ….. body against the same protein. Western blot and ELISA require high amount of protein lysates Antibody 4 ….. as compared to the RPPA technology [BB15, CH15].
	Tissue
		Signal
		Antibody
	Signaling Network Analysis	Data Analysis
	Figure 2.3 -for investigating and comparing samples treated at different time points under different condi-
	tions with different doses. It is a useful tool to quantify multiple phosphorylated proteins even
	in a very small sample, which makes them suitable for individual patient therapy. Several stud-
	ies have investigated heterogeneity among protein levels within a primary tumor, tumor and
	metastases of the same patient using data generated by the RPPA technology. Results showed
	significant heterogeneity in a subset of proteins within a tumor, and between metastases and the
	primary tumor. This suggests to analyze the samples from multiple locations instead of a single
	sample. Another study used RPPA to distinguish colon and ovarian cancer by finding molec-
	ular markers. Other studies found key proteins which play important role in the regulation of
	signaling pathways in different cancers such as breast cancer, leukemia, and glioma. A first
	commercial RPPA assay (the TheraLink HER Family Assay) was developed by Theranostics

Table 2 .

 2 2 -Characteristics of different protein analysis platforms. OPS + 15, MTH + 12, GVE + 13, SRAE + 09, VGE+ 12]. In this thesis, we use a PKN (built using ReactomeFIViz) combined with the phosphoproteomic time series dataset of four breast cancer cell lines to generate cell line specific Boolean Networks.

		Pros	Cons
	Western Blot	High specificity	High sample, Low-medium throughput
	ELISA	Quantitative, High sensitive	High sample, Costly setup for high-throughput
	Mass Spectrometry High multiplex, Discovery of new biomarkers Low throughput, Complex sample preparation
	RPPA	High sensitive, High-throughput	Costly setup, Specific antibodies
	In this thesis, we are using the HPN-DREAM dataset, which was generated by RPPA quan-
	titative proteomics technology [HNJC + 17, HHC + 16]. These data contain phosphorylated mea-
	surements of multiple proteins under sets of perturbations. Perturbation refers to the combina-
	tion of stimuli and inhibitors.
	2.3 Prior knowledge networks
	Prior Knowledge Networks (PKNs) are available in different databases such as Reactome,
	PID, and kegg among others [DW15, WDD + 14a, KG00, C + 04, KvIH + 12, Nis01, SBR + 06,
	SMC + 17, XRS + 00, PNK + 04, HMPL + 04, ZMPQ + 02, RMD08]. PKNs are graphs where
	molecules are represented by nodes and interactions are represented by edges. We can con-
	struct a PKN through different tools or softwares such as ReactomeFIViz [WDD + 14b] which
	is available as a Cytoscape [SMO + 03] plugin. A PKN alone cannot be used to build reliable
	dynamical models or to explain underlying biological behaviors [RCAdS15], because signaling
	behaviors are rewired in specific contexts. The signaling behavior may differ in cancerous and
	normal cells due to many genetic modifications [HLM + 12]. Therefore, it is extremely impor-
	tant to manifest how these networks are regulated in different diseases. In order to overcome

Table 3 .

 3 1 -Basic connectives to form logic programs

		Logic Program
	if	←
	and	∧
	default negation	¬
	A logic program consists of rules of form

  and each b i is an atom. If a rule contains variables, we This program consist of three atoms switchOn, switchOFF and light. The first rule states that if the switch is not off then the switch is on. The second rule states that if the switch is not and the second solution consists of two atoms {switchOn, light}. Let's inspect the first solution. Clearly, all rules are satisfied by this model but only the body of the second rule is satisfied. Hence, the reduct consists of the second rule only. Since the empty set does not satisfy the reduct, {switchOff } is a stable model. The second solution can be verified similarly but this time the reduct consists of the first and the third rule, which means that both switchOn and light have to be included in the model of the reduct.

	switchOn ← ¬switchOff	(3.1)
	switchOff ← ¬switchOn	(3.2)
	light ← switchOn	(3.3)

obtain a set of ground rules by substituting all possible ground terms for the variables. Given a logic program, the corresponding ground logic program is the set of all possible ground rules obtained from its rules. Such a (ground) logic program induces a set of stable models determined by the stable model semantics

[START_REF] Gelfond | The stable model semantics for logic programming[END_REF]

. Each stable model is a subset of the atoms occurring in the logic program. Atoms appearing in this set are said to be true, and false otherwise. A rule is satisfied if its body (the part after the ←) is not satisfied, or its head atom h is true. A rule body is satisfied if all the atoms b 1 to b m are true and all the atoms b m+1 to b n are false. A rule with an empty body is called a fact. A model satisfies all (ground) rules of a logic program and a stable model also satisfies a minimality criterion. This criterion requires that each atom in a stable model is proved by some rule. For this, a true atom has to appear in at least one rule head with a satisfied body. Or more formally, a model X is a stable model if there is no smaller model Y ⊂ X satisfying the reduct of the program; the reduct comprises all rules whose bodies are satisfied by X. In the following, we simply refer to the stable models of a logic program as its solutions. Next, we give a simple example of a logic program 2 .

Example 3.2.1. Here, we model the light switch scenario. on then the switch is off. The last rule states that if the switch is on then there is a light. This problem consists of two stable models as solutions. The first solution consists of one atom {switchOff }

  ≥ 0, the lower bound l and upper bound u are integers, and each element e i has form t 0 , . . . , t o : l 1 , . . . , l m for o ≥ 1, m ≥ 0, each t k is a term and each l j is a literal. Unlike cardinality aggregates, sum aggregates allow for expressing that the sum of the integers given by a set of tuples has to be between a lower and upper bound. An aggregate element is satisfied if the literals l 0 to l m are satisfied. In this case, the term tuple t 0 , . . . , t o contributes the integer t 0 to sum up3 . Unlike cardinality aggregates, this form of aggregate can only be used in a rule body. Like with cardinality aggregates, at least one bound has to be specified.

Example 3.2.4. The program 1 {switchOn; switchOff } 1 ← light ← switchOn has the same solutions as Example 3.2.1. Here, the first rules states that exactly one of the atoms switchOn and switchOff has to be true. Sum aggregates A sum aggregate has form l sum{e 1 ; . . . ; e n } u where n Example 3.2.5. The following program extends the light switch scenario:

{on(s 1 ); on(s 2 ); on(s 3 )} ← light ← 2 sum{1, s 1 : on(s 1 ); 1, s 2 : on(s 2 ); 1, s 3 : on(s 3 )}

  3.2 -Directed graphFigure 3.2 is defined by the following facts: We use rule 3.4 to define that each node can have a different color. So predicate color /1 corresponds to variable y in the ILP formulation of the graph coloring problem. Another rule (3.5) is added to assign exactly one color to each node. Here, predicate assign/2 corresponds to variable x in the ILP formulation and the rule captures Equations 2.11 and 2.12. Next, rule 3.6 is added to ensure that two adjacent nodes must not have the same color. This rule corresponds

	node(1..6) ←
	edge(1, 2) ←	edge(1, 3) ←	edge(1, 4) ←
	edge(2, 4) ←	edge(2, 5) ←	edge(2, 6) ←
	edge(3, 1) ←	edge(3, 4) ←	edge(3, 5) ←
	edge(4, 1) ←	edge(4, 2) ←
	edge(5, 3) ←	edge(5, 4) ←	edge(5, 6) ←
	edge(6, 2) ←	edge(6, 3) ←	edge(6, 5) ←

to Equation

2

.13 in the ILP formulation. Rule 3.7 is used to eliminate isomorphic solutions w.r.t. the selection of colors

Table 3 .

 3 2 -Comparison of ASP solver and CellNOpt as presented by Videla et al.

		Mid-scale case study	Large-scale case study
		ASP solver	CellNOpt	ASP solver	CellNOpt
	Execution Time	0.03 s	9.2 h	0.07 s	27.8 h
	Best Solution	0.06 s (optimal)	7.2 h	0.04 s (optimal)	24.5 h
	Number of models	8 (optimal)	66	2 (optimal)	206
	Size of model	16	16 to 24	26	27 to 36
	Ostrowski et al.

Method description. Ostrowski et al. [OPS

+ 15]

proposed a new method (caspo time series:

  2.1 , C 2.2 , C 2.3 ) were created by modifying the experimental design. Case study C 2.1 consists of 14 nodes and 22 edges, C 2.2 consists of 16 nodes and 25 edges, and C 2.3 consists of 40 nodes and 58 edges. Third case study C 3 represents the ERBB receptor regulated G1/S transition model ([SFL + 09]. Case study C 3 consists of 19 nodes and 50 edges. Finally, they utilized two optimization schemes (cardinality and subset) of the ASP solver to derive minimal models for these case studies.

Table 4 .

 4 1 -Computation summary. Here, we show the number of verified solutions, true positive and false positive BNs, and their computation (ASP solving and Model Checking steps) time for each cell line. The ASP solving was performed on a standard laptop machine. The model checking task was performed on a cluster with 560 cores and 1.9 Tb of RAM.

	Cell Line Number of Solutions True Positives False Positives	Time
					ASP solving Model Checker
	BT 20	188	72	116	210 seconds	< 7 days
	BT 549	231	191	40	93 seconds	< 7 days
	M CF 7	52	21	21	36 seconds	< 7 days
	U ACC812	150	0	150	197 seconds	7 days

Table 4

 4

.2. Figure 4.11 agrees with the results presented in Table 4.2 as we can see the clear discrepancies among the four cell lines. It

Table 4 .

 4 2 -Similarity scores among breast cancer cell lines.

	Cell Lines Size of BNs' family		Similarity Score
			BT20 BT549 MCF7 UACC812
	BT 20	72	0.73	0.15	0.17	0.23
	BT 549	191	* *	0.43	0.23	0.20
	M CF 7	21	* *	* *	0.63	0.21
	U ACC812	20	* *	* *	* *	0.46

Table 4 . 3 -

 43 Root mean square error. This table summarizes the RMSE results for each cell line. We have calculated the discrete RMSE (error related to the discretization of the data) and the model RMSE (caspo-ts error). The Delta column shows the difference between model and discrete RMSE.

	Cell Line	Learning		Testing
		Discrete Model Delta Discrete Model Delta
	BT 20	0.3464 0.3464	0	0.3293 0.3302 0.0009
	BT 549	0.3498 0.3498	0	0.3007 0.3113 0.0106
	M CF 7	0.3207 0.3207	0	0.2772 0.2772	0
	U ACC812 0.3464 0.3464	0	0.3084 0.3178 0.0094
	Analyzing family of BNs using root mean square error

Table 4 .

 4 4 -Computation summary. This table summarizes the analysis performed on the BT549 cell line. The # of solutions column shows the number of ASP optimal BNs. The true and false positive columns show the number of true positive and false positive BNs returned by the model checker. The Delta RMSE column shows the difference between the discrete and model RMSE. The ASP solving column represents the time spent on learning ASP optimal BNs for each sample. The model checking column shows the time spent on verifying the BNs for each sample. The last column RMSE represents the corresponding RMSE of each sample with respect to the testing data. The ASP solving was performed on a standard laptop machine. The model checking task was performed on a server with 1.5 Tb of RAM.

	BT549	# of Solutions		Learning Data		Time		Testing Data
	Sample #		True Positives False Positives Delta RMSE ASP solving Model Checking	RMSE
	1	231	231	0	0	124 seconds	18 days	0.3113
	2	231	231	0	0.0009	125 seconds	15 days	0.3113
	3	231	0	231	0.0002	122 seconds	3 days	0.3113
	4	231	0	231	0.0004	112 seconds	10 days	0.3113
	5	231	0	231	0.0007	101 seconds	7 days	0.3113
	6	231	0	231	0	110 seconds	9 days	0.3113
	7	231	0	231	0	133 seconds	5 days	0.3113
	8	231	0	231	0	133 seconds	6 days	0.3113
	9	231	0	231	0.0003	134 seconds	6 days	0.3113
	10	231	0	231	0	136 seconds	18 days	0.3113
	11	231	0	231	0.0001	90 seconds	22 days	0.3113
	12	231	0	231	0.0029	148 seconds	30 days	0.3113
	13	231	0	231	0	104 seconds	30 days	0.3113
	14	231	0	231	0.0003	91 seconds	23 days	0.3115
	15	231	0	231	0.0001	329 seconds	22 days	0.3113
	16	231	120	111	0.0018	100 seconds	25 days	0.3113
	17	231	0	231	0	150 seconds	30 days	0.3113
	18	231	0	231	0.0004	131 seconds	30 days	0.3115

  This method combines Answer Set Programming and Model Checking techniques to infer true positive BNs verifying the experimental data. Caspo-ts allowed us to handle a midscale PKN (64 nodes and 178 edges) and a real dataset subject to experimental error. Caspo-ts enabled us to learn key dynamic mechanisms within the BNs explaining the time series data. Our results

  t.an objective function. The model checking part of the system detects TP BNs by checking the reachability of time series traces given compatible BNs generated by the ASP part of the system. TP BNs are guaranteed to reproduce all the (binarized) traces under all perturbations by verifying exact reachability in the BN state graph. Since checking this reachability is a PSPACE-hard problem, the second step can be very time consuming for large BNs. To resolve this issue, the ASP part over-approximates solutions [OPS + 16]. This over-approximation removes a large set of BNs that have no reachable traces, reducing the number of calls to the model checker.

Table 5 .

 5 1 -Number of TP BNs out of 46 BNs, time to get the first TP solution, and similarity among TP solutions per cell line.

			caspo-ts			caspo-ts D
	Cell Line TPs Time Similarity TPs Time Similarity
	MCF7	0	-	-	4	6.7h	0.51
	BT549	2	8.4m	0.92	14 7.9m	0.44
	UACC812 20	26s	0.81	15	27s	0.45
	BT20	13	20h	0.86	7+	20h	0.32

Now, we analyze in more detail the UACC812 cell line using caspo-ts. This cell line was the most difficult one as explained in Chapter 4, Section 4.5.4.

Table 5 .

 5 2 -Comparison of previous and parallel model checking. BNs of TCR signaling dataset by exploring only 100 BNs using caspo-ts D . We also identified TP BNs in a case (MCF7) where caspo-ts could not find any TP BN among 46 solutions. Our method is applicable to gene or protein expression time series datasets measured upon different perturbations. Moreover, the proposed Algorithm (4) is not specific to our biological application. It computes diverse subset minimal solutions in ASP,

	BT549	Previous Model Checking	Parallel Model Checking	Difference
		Verified BNs TPs FPs Time Verified BNs TPs FPs Time
	1	231	0	231 64d	231	0	231 7.2d	56.8d
	2	221	0	221 >90d	231	0	231 15.7d	>74d
	3	173	0	173 >90d	231	0	231 11.9d	>78d
	4	231	0	231 7.2d	231	0	231 17h	6.5d
	5	37	0	37 >90d	231	0	231 20.4d	>70d
	6	231	0	231 78.5d	231	0	231 14d	64.5d
	7	22	0	22 >90d	231	0	231	6h	>84d
	8	172	0	172 >90d	231	0	231 8.5d	>82d
	9	231	0	231 76.7d	231	0	231 9.5d	67.2d
	10	121	0	121 44d	231	0	231 2.5d	41.5d
	11	135	0	135 >90d	231	0	231 7.5d	>82d
	12	231	231	0	16.2d	231	231	0	4.1d	12.1d
	13	231	0	231 9.3d	231	0	231 10h	8.9d
	14	231	231	0	50d	231	231	0	11.8d	38.2d
	15	231	0	231 40d	231	0	231 3.8d	36.2d
	16	105	0	105 >90d	231	0	231 29.7d	>60d
	complete solution space of 68338

This example was inspired by [BL13].

Note that the model checker could only verify 32 out of 46 solutions within one month for cell line BT20 in case of caspo-ts D . There may exist more TPs for this cell line.

List of Figures Computational Discovery of Dynamic Cell

Line Specific Boolean Networks from Multiplex Time-Course Data

Introduction

In this chapter, we improve and apply the caspo time series (caspo-ts) method conceived to identify optimal BNs from time-series phosphoproteomic measurements and a prior knowledge network. The article related to this work is published in the Plos Computational Biology journal [RPS + 18].

We use a large-scale real-case study from the HPN-DREAM challenge. This dataset contains phosphoproteomic time series data of four Breast Cancer cell lines. Here, we show that the discrete approach, based on Answer Set Programming and model-checking techniques, provides accurate results. We identify optimal networks explaining the dynamic behavior of the experimental data for a network of 64 components and 178 interactions. Our models have a Root Mean Square Error (RMSE) of 0.31 with respect to the testing data. It is well known that the methods addressing dynamic properties on large networks do not scale well because of the considerable size of the respective state transition graph. Using caspo-ts, we identify the BNs given the inherent combinatorial explosion that appears when considering real case data. Our key findings suggest that this method is capable of constructing cell line specific Boolean models, which is extremely valuable given the heterogeneity of breast cancer due to many genetic

Validation using random data samples

The objective of the following analysis is to show that the BNs obtained with caspo-ts using the HPN-DREAM datasets for the four cell lines have a worse RMSE score with respect to random trajectories, and therefore are very specific to the HPN-DREAM datasets. For this purpose, we generated 100 random data samples per cell line. In each sample, we generated a random value in [0; 1] for each readout protein in each time point without changing the perturbations. Then, we calculated the RMSE of these samples with respect to the inferred BNs of each cell line, and finally compared it with the learning and testing RMSE of these BNs. In Figure 4.15, the RMSE ratio for random datasets is displayed by red boxplots for each cell line, and the RMSE ratio for testing and learning datasets is shown as clear outliers in green and blue colors respectively. It is worth noting that the caspo-ts method has failed to recover random data time series, hence proving the specificity of the learned networks with respect to the HPN-DREAM challenge dataset.

Validation by introducing noise in the learning data

Additionally, we computed the RMSE of the testing data by using a leave one out approach.

The analysis is performed for the cell line BT549. Recall that checking the exact reachability of the BNs with respect to a binarized time series trace is a computationally expensive task. The time to model check BNs differs from one network to another. For this case study, the model checking time varied from minutes to days. Because of this reason, we have done this analysis for only one cell line.

We generated 18 modified learning datasets by randomly introducing 5% of noisy data in each dataset; that is, the 5% of the total time point values were generated randomly. We used these datasets to learn 231 ASP-optimal BNs. Then we passed these BNs to the model checker to detect true positive BNs.

Results

We discuss the results of applying Algorithm 4 on two different datasets. We start with the artificial benchmark, where the solution space is small enough (68338) to compute all solutions. This allows us to study this benchmark in more detail. We can analyze how well a limited number of solutions enumerated with caspo-ts and caspo-ts D represents the solution space.

Then we move to the real dataset, where we cannot enumerate all solutions and can only consider a limited number of solutions because the solution space is too large. Nevertheless, we can show improved results with caspo-ts D over caspo-ts by being able to enumerate more TP BNs and also more diverse BNs.

We also discuss the results of applying Algorithm 5 to the BT549 cell line. For this, we compare the computation time of the caspo-ts system based on previous and new implementation of the model checking process. We can show that we spend less time in verifying reachability, especially in case of FP BNs, as we might be able to stop the verifying process as soon as one of the perturbation fails to verify.

Artificial dataset

Here, we use the TCR signaling dataset [KSRL + 06] we were able to find solutions for cell line MCF7, where caspo-ts could not find solutions before. Furthermore, we have presented an algorithm to parallelize the model checking process.

The results show substantial improvements in computation time of the model checking process.

CHAPTER 6. CONCLUSIONS AND FUTURE WORK

The high false positive rate of BNs highlighted issues existing in a real dataset, i.e., missing time points and inconsistent or incompatible measurements of proteins (see Chapter 4 Section 4.5.2). It also pointed us to the fact that a limited number of perturbations is not enough to derive a unique model representing one specific cell line, but rather leads billions of solutions (BNs). We also discovered that the shuffling of perturbations leads to the generation of solutions in a different order. By resolving these problems, we managed to generate true positive BNs for the case study of the HPN-DREAM challenge. Our results point to measurements in the time series phosphoproteomic dataset that contradict the experimental setting and to perturbations that show contradictory dynamics. This analysis highlighted the erroneous experiments, so this method can be used as complementary to the existing approaches to discover errors (see Chapter 4).

A detailed analysis of the large solution space generated by the caspo-ts method led to the discovery of clustered solutions. The caspo-ts method uses the clingo ASP solver, which is able to exhaustively enumerate all solutions. The clingo solver by default uses an enumeration scheme, in which, once a solution is found, it backtracks to the first point from where the next solution can be found. This typically leads to the situation where successive solutions only change in a small part. As a result, caspo-ts may enter a solution space where BNs are clustered together. We have observed that given the size of the PKN and the small number of perturbations in the experimental data, the solution space of the caspo-ts can be very large containing billions of BNs making it difficult to enumerate true positive BNs in reasonable time if it gets stuck in a cluster of false positive BNs. We have improved the ASP-based caspo-ts system to resolve the above mentioned issue. We have extended the caspo-ts system using heuristics to enumerate diverse solutions. Our results showed that this extension proved to be very useful to breakup cluster of solutions and to get the better impression of the solution space (see Chapter 5).

The high computation time in case of a real case study (HPN-DREAM challenge) directed us to analyze the model checking component of the caspo-ts system. Since the caspo-ts system uses an over-approximation criteria, it can lead to false positive BNs. These false positive BNs are ruled out using a model checker by verifying a reachability property. We observed that the properties are verified in a manner, which can cost the expensive computation time spent on verification of properties. To resolve this issue, we proposed a new parallel verification of the properties. Our results showed that the parallel verification of properties saves us months of computation time (see Chapter 5).

FUTURE PERSPECTIVES

Conference paper