BioSystems 150 (2016) 52–60

Contents lists available at ScienceDirect

BioSystems journal homepage: www.elsevier.com/locate/biosystems

Networks and games for precision medicine Célia Biane, Franck Delaplace ∗ , Hanna Klaudel IBISC Laboratory, Evry Val d’Essonne University, Evry, France

a r t i c l e

i n f o

Article history: Received 17 February 2016 Received in revised form 20 July 2016 Accepted 11 August 2016 Available online 16 August 2016 Keywords: Precision medicine Network medicine Game theory Boolean network Therapy prediction Systems biology Cancer therapy

a b s t r a c t Recent advances in omics technologies provide the leverage for the emergence of precision medicine that aims at personalizing therapy to patient. In this undertaking, computational methods play a central role for assisting physicians in their clinical decision-making by combining data analysis and systems biology modelling. Complex diseases such as cancer or diabetes arise from the intricate interplay of various biological molecules. Therefore, assessing drug efficiency requires to study the effects of elementary perturbations caused by diseases on relevant biological networks. In this paper, we propose a computational framework called Network-Action Game applied to best drug selection problem combining Game Theory and discrete models of dynamics (Boolean networks). Decision-making is modelled using Game Theory that defines the process of drug selection among alternative possibilities, while Boolean networks are used to model the effects of the interplay between disease and drugs actions on the patient’s molecular system. The actions/strategies of disease and drugs are focused on arc alterations of the interactome. The efficiency of this framework has been evaluated for drug prediction on a model of breast cancer signalling. © 2016 Elsevier Ireland Ltd. All rights reserved.

1. Introduction The analysis of the patients’ omics profiles (genome, metabolome, proteome, etc.) will soon become a standard for customized molecular-based diagnosis and treatment by contrast to a “one-size-fits-all” strategy based on a one-to-one correspondence between diseases and drugs (Ginsburg, 2009; Mirnezami et al., 2012). Precision medicine is an emerging branch of medicine that uses omics data to improve clinical decision-making by designing new tools for the personalization of therapies and their risk/benefit assessment. Addressing this challenge puts the focus on the causality study of the pathogenesis at molecular level. However, the causal relationship between omics information and disease phenotypes remains elusive. Indeed, a disease phenotype is rarely a consequence of an abnormality in a single gene product but involves complex interplays of various biological molecules (Barabási te al., 2011). For instance, patients with sickle cell anaemia, which is caused by a unique well-defined genetic event in a single gene (classic Mendelian disease) can exhibit highly variable phenotypes in the clinic (Ballas, 2011; Schadt, 2009). This variability is due to the interaction of the mutated

∗ Corresponding author. E-mail addresses: [email protected] (C. Biane), [email protected] (F. Delaplace), [email protected] (H. Klaudel). http://dx.doi.org/10.1016/j.biosystems.2016.08.006 0303-2647/© 2016 Elsevier Ireland Ltd. All rights reserved.

gene with other individual-dependant genetic variants (Sebastiani et al., 2010; Schadt, 2009). Therefore, understanding the pathogenesis at molecular level requires to conceive frameworks facilitating the discovery of causes altering the molecular system of a living organism. This challenge logically focuses on biological networks modelling the causal interplays of molecules (Delaplace et al., 2010). The main approaches in this field study the location of dysfunctional molecules in networks and the nature of network alterations leading to diseases. The works (Barabási, 2007; Barabási te al., 2011; Gustafsson et al., 2014) study the formation of specific molecular subnetworks, called modules, supporting diseases. This approach is motivated by the hypothesis that network modules support key molecular functions disrupted in disease (Davidson, 2010; Milo et al., 2002). In Oti et al. (2006), authors confirm the fact that proteins involved in the same disease have a higher propensity to interact with each other, forming a tightly interconnected entity in the interactome. Thereby, disease should likely alter a functional module or being themselves modules supporting a dysfunctionality (disease modules). In Zhong et al. (2009), a network-perturbation approach is used to explain molecular dysfunctions underlying human diseases. The genetic events causing diseases are expressed as perturbation of both edges and nodes of the interactome. Schematically, a genetic event leading to the expression of an inoperative protein is modelled by a node deletion while genetic events inducing loss or gain of interaction are respectively modelled by an edge deletion or

C. Biane et al. / BioSystems 150 (2016) 52–60

53

The paper is structured as follow: after recalling the main features of Boolean networks and Ordinal Game Theory, we detail the theoretical framework called Network-Action Game in Section 2 and show its application to drug prediction in breast cancer in Section 3. 2. Network-Action Game We first introduce the two theories composing the NetworkAction Game framework: Boolean networks and Ordinal Game Theory and we then show their coupling. 2.1. Boolean networks

Fig. 1. Overview of the Network-Action Game method.

A Boolean network is a discrete dynamical system of a population of agents defined by a family of propositional formulas determining the evolution of the agents. The dynamics is defined by a transition system on a set S representing all the possible states of the agents and transitions represent their evolution. Let A be a set of agents, a (Boolean) state is defined as a mapping s : A → B associating to an agent in A a value from B and S denotes the set of mappings representing the states. Let F = (fa )a ∈ A be a family of propositional formulas, each fa is a function defining the state of a from the states of all the agents (seen as propositional variables). An asynchronous1 Boolean network is defined as a pair A, F and its model of dynamics is a labelled transition system (S, −→, A) where the transition relation labelled a

addition (edgetic perturbation). They uncovered experimental and computational evidences that these network alterations occur in human Mendelian diseases. It is worth noting that the perturbation on networks induced by diseases are formalized by elementary topological modifications of molecular interaction networks: nodes or edges are added or deleted. Hence, the complexity of disease should rely on the impact of these topological perturbations on the physiological processes controlled by these networks. For instance, cancer cells acquire the ability to maintain proliferative signalling notably by defecting feedback loops regulating cell division (Hanahan and Weinberg, 2011). Therefore, a deeper understanding of disease/therapy mechanisms relies on the prediction of the consequences of these elementary actions on the dynamics of molecular networks. In this paper, we combine two theoretical frameworks: Game Theory and discrete models of dynamics (Boolean networks) to determine the best drug to administrate to a patient. The clinical decision-making is modelled using Game Theory, that defines the process of selection by the players of an action among alternative possibilities (Chettaoui et al., 2007; Osborne and Rubinstein, 1994), while Boolean networks are used for modelling the effects of the interplay between disease and drugs on the patient molecular system. Boolean networks are used in biology to study the dynamics of molecular networks (modelled as interaction graphs), which represent functional interactions between molecules (Abou-Jaoudé et al., 2009; Ciliberto et al., 2005; Thomas and Thieffry, 1995). Such a dynamics evolves towards equilibria interpreted at the molecular level as the patient health condition or illness. The physician and the disease are considered as players of a game, each of them having strategies of action that correspond to a drug administration and to a genetic event, respectively. In a game, combinations of strategies, called strategy profiles, modify the patient interaction graph, therefore modifying the associated Boolean dynamics and its equilibria. From the assessment of biomarkers at these equilibria, players’ preferences are determined, and then, the interpretation of the Ordinal Nash equilibrium leads to the discovery of the best physician action (drug). Fig. 1 recalls the main steps of the framework described above.

by agent a and denoted −→ ⊆ S × S, updates the state of agent a only. Formally, it is defined as: a

def

s−→s = (s [a] = fa (s) ∧ s[a] = / fa (s) ∧ (∀x ∈ A\{a} : s[x] = s [x])). Hence the global dynamics is the union of the transitions labelled



a

by agents (i.e., −→= a ∈ A −→). The signed interaction graph associated to F, G = A, →, ı represents all the signed interactions defined by F between agents in A. The sign of the arc is given by a labelling function ı and may be + for increasing relation, − for decreasing one and ± otherwise. Such a graph can be inferred from the syntax of the propositional formulas −

in minimal2 disjunctive normal form, where ai →aj stands for the +

occurrence of negative literal ¬ai in faj , ai →aj for the occurrence of ±

positive literal ai in faj , and ai →aj for both. Fig. 2 shows a Boolean network, its model of dynamics and signed interaction graph. A state s is an equilibrium for −→, if it may be reached infinitely often, i.e., ∀s ∈ Bn : s−→∗ s ⇒ s −→∗ s, where −→* denotes the reflexive and transitive closure of −→. We denote by E−→ the set of all equilibria of −→. An attractor is a set of equilibria that are mutually reachable and a steady state is an attractor of cardinality 1. In Fig. 2, the states (1, 0, 0) and (0, 1, 1) are steady states. 2.2. Ordinal game An Ordinal Game models strategic decision-making based on the definition of a preference relation amongst combination of players’ strategies. Each player has a set of possible strategies and a strategy profile represents a particular combination of strategies. A preference relation defines the preference, for a player, between each pair of strategy profiles and a preference graph represents the union of the preference relations of all players. Formally, an Ordinal Game is a triple P, (Cp )p ∈ P , ( p )p ∈ P  where:

1 Asynchronous means that the state of at most one agent is updated at each transition. 2 In a minimal disjunctive normal form the conjunctions are the prime implicants.

54

C. Biane et al. / BioSystems 150 (2016) 52–60

Fig. 2. A Boolean network, its model of dynamics and associated interaction graph.

• P = {p1 , . . ., p|P| } is a set of players; • for each player pi ∈ P with i ∈ [1, . . ., |P|], Cp is a non-empty set of i strategies Cpi of player pi . The set CP = Cp1 × . . . × Cp|P| represents the set of all strategy profiles; • for each player pi ∈ P with i ∈ [1, . . ., |P|], the relation p ⊆ CP × i CP is a transitive and reflexive relation (pre-order) on strategy profiles called preference relation of player pi . In the sequel, we denote by c−p the strategy profile c restricted to strategies of players in P \ {p}. For each player, the preference relation is restricted to strategy profiles where only his/her own strategy may differ:

∀p ∈ P : (c p c / c−p (c−p =

) ∧ ⇒ c−p = c−p

⇒ (c / p c ∧ c / p c)).

(1)

The solution concept of an Ordinal Game is an Ordinal Nash equilibrium, namely a strategy profile c such that, by considering the opponents’ choices, a player may only deviate into another Nash equilibrium. Formally, a strategy profile c ∈ CP is an Ordinal Nash equilibrium if:

∀c ∈ CP , ∀p ∈ P : c p c  ⇔ c  p c. In Le Roux et al. (2008), the authors propose a computational definition of the Ordinal Nash equilibrium as the strategy profiles belonging to a terminal strongly connected component of the preference graph. A preference graph represents  the union of the  preference relations of all players =

where the preferp p∈P ences are represented as arcs between strategy profiles and labelled by the player. As an illustration, let us take the example of the prisonner’s dilemma. In this famous example, two suspects, P = {p1 , p2 }, are arrested without the possibility to interact and are offered to choose between two strategies: betray (B) by testifying that the other committed the crime or remain silent (S), i.e., both players have the same sets of strategies Cp1 = Cp2 = {B, S}. A prison sentence is associated to all the possible strategy profiles, i.e., CP = {(B, B), (B, S), (S, S), (S, B)}: if both players betray (B, B), they serve 2 years in prison, if they both remain silent (S, S), they serve 1 year, if p1 betrays but p2 remains silent (B, S), p1 is set free and p2 serve 3 years and vice versa for (S, B). The determination of preference relations between strategy profiles is governed by the expectation that players would minimize their prison sentence. The preference graph of the prisoner’s dilemma is shown in Fig. 3, where the Nash equilibrium is highlighted in grey and corresponds to strategy profile (B, B) where both players betray each other. This strategy is not optimal for each player individually but none of them has an incentive in changing its strategy while considering at the other’s strategy. 2.3. Network-Action Game Our framework, called Network-Action Game, defines an Ordinal Game complying to two founding principles.

Fig. 3. Preference graph of the prisoner’s dilemma. The nodes correspond to strategy profiles and the arcs to the preferences labelled by players. The Nash equilibrium is highlighted in grey.

In the first principle, each strategy profile c ∈ CP is associated to a Boolean network A, Fc  giving rise to the unique dynamics (S, −→ c ) and its set of equilibria E−→c . The preference relation p between any strategy profiles c and c , for each player p, is deduced from an assessment of these equilibria. More precisely, p is defined from a pre-order p on sets of equilibria E−→c and E−→c , as follows: def

∀c, c ∈ CP : c p c = E−→c p E−→c ∧ c−p = c−p .

(2)

The second principle is directly related to the biological network perturbation scheme mentioned in the Introduction. Actually, the effects of a strategy profile c ∈ CP are instantiated as basic actions of arc deletion or maintenance on the interaction graph GF of a saturated Boolean network A, F. By definition, GF must gather all the arcs that will be deleted or maintained during the game. The application of each strategy profile c gives rise to a Boolean network A, Fc  computed from the saturated Boolean network by replacing each ai by 0 in faj if arc ai → aj is deleted by c, and by maintaining ai in faj otherwise. Thus, an arc deletion is interpreted as the absence of the source ai for the target aj . This rule is motivated by the biological application detailed in Section 3. 2.4. Example The following example illustrates a practical methodological framework for Network-Action Game application (cf., Section 3) decomposed into four stages detailed below. Let us consider a two-players Network-Action Game with a given saturated Boolean network A, F. For the example, we choose the Boolean network A, F = {a1 , a2 , a3 }, (a1 = ¬a2 ∧ ¬a3 , a2 = ¬a1 , a3 = a2 ). presented in Fig. 2. Each player has two strategies: a particular one, ˛ for p1 , ˇ for p2 , and an identical one  for both. Stage 1 We define first a reference Boolean network, obtained by application of strategy profile cref = (, ) on the saturated Boolean network A, Fcref . For the example, the interaction graph of A, Fcref , Gref is depicted in the top left square of Fig. 4. The notion of

C. Biane et al. / BioSystems 150 (2016) 52–60

reference is a matter of modelling convention, introduced here to clarify its application in Section 3. A reference indicates a particular Boolean network representing a basal condition and resulting from the application of a strategy profile with a strategy  common to all players. Stage 2. The actions of strategies are defined by arc action function ϕ assigning to each arc of GF a value defining a strength (cf., Table 1). The positive value is interpreted as arc maintenance, whereas the negative or null one as arc deletion from GF . It is worth noticing that an arc absent in Gref and for which the actions of strategies induce its maintenance from GF is added into Gref , by convention. Stage 3. In order to determine the actions of strategy profiles c ∈ CP in terms of arc deletion or maintenance in GF , we use function ϕ to compute the strength from the combination of these strategies,

55

Table 1 Strength ϕ of arcs per strategies.

 ˛ ˇ

a1 → a2

a2 → a1

a3 → a1

a2 → a3

1 1 1

−1 2 2

1 −2 1

1 1 −2

from which the resulting arc action is deduced. In our example, CP = {(, ), (˛, ), (, ˇ), (˛, ˇ)} and for each c ∈ CP we compute, for each arc ai → aj , the sum c,ai →aj =



ϕ(cpi , ai → aj )

pi ∈ P

Fig. 4. Four strategy profiles c ∈ {(, );(˛, );(, ˇ);(˛, ˇ)} with associated families of formulas Fc , interaction graphs, models of dynamics (S, −→ c ) and equilibria E−→c (in grey).

56

C. Biane et al. / BioSystems 150 (2016) 52–60

Table 2 Function pi . The table describes the scores per player defined from the computed equilibria for each strategy profile. Players/equilibria

p1 p2

E−→(,)

E−→(˛,)

E−→(,ˇ)

E−→(˛,ˇ)

011

100

011

100

100

010

100

010

2 1

1 2

2 1

1 2

1 2

1 2

1 2

1 2

Fig. 5. Preference graph of the example. The nodes are strategy profiles and the arcs are the preferences labelled by players. The Ordinal Nash equilibria are highlighted in grey.

of ϕ for all Cpi composing c. Thus, if c,ai →aj > 0, then arc ai → aj is maintained, otherwise it is deleted. Stage 4. Finally, the players’ preferences on strategy profiles are computed. They are determined from a scoring function on states, pi : S → R for each player and the score is computed from the markers’ states at steady states. To this end, the steady states E−→c are computed from each Fc and the set of markers is fixed. In the example, the scoring function p1 corresponds to the number of 1s in a state while p2 is the number of 0s. The set of markers equals the set of all agents, M = A. The four strategy profiles give rise to the families of formulas, models of dynamics and equilibria depicted in Fig. 4 and the scores at steady state are given in Table 2. Therefore, for a strategy profile c, a set of scores is computed from the application of the scoring function on the steady states of the Boolean network A, Fc . From these sets, a preference relations

pi is determined for each player as specified in Equation (2). In the example, we compare the scores according to the following formula, where c and c are strategy profiles: def

c pi c =

max({pi (s) | s ∈ E−→c }) ≤ min({pi (s) | s ∈ E−→c })

(3)

. ∧ c−pi = c−p i

This preference selects the best profiles, meaning that player p prefers c to c whenever the scores obtained by c are all less or equal to the scores obtained by c . For the example, the union of the preference relations of p1 and p2 is represented as a preference graph depicted in Fig. 5. The interpretation of the preference graph is application dependent and will be discussed in Section 3. 2.4.1. Complexity of the method The evaluation of the Network-Action Game complexity depends on the method used to determine the preference relation. To some extent, the example can be considered as a methodological framework for the application of Network-Action Game. To evaluate its complexity we set: n as the number of agents, p as the number of players and m as the (maximal) number of strategies per player. Assume that the maximal number of substitution of variables by 0 in the formulas corresponding to deletion on the saturated Boolean network is ı. Moreover, checking whether a state s is stable (i.e.,

F(s) = s) by a brute force algorithm requires to enumerate all the 2n states. Thus, applied to the mp strategy profiles, the computation of all the scores is then in O(mp .ı.2n ). The computation of Ordinal Nash equilibria requires to build the preference graph by pairwise comparison of the mp scores of profiles (m2p comparisons). It corresponds to the computation of the terminal strongly connected components using Tarjan algorithm (Tarjan, 1972) which is linear in the number of elementary preferences (arcs in the preference graph) also bounded by m2p . Thus, the complexity of the method is in O(mp · ı · 2n + m2p ) representing the score computation for each profile followed by the preference graph construction. However, it is worth noting that the exponential complexity does not reflect its performance in practice. For instance, the computation time3 for the application in Section 3 is 0.3 s with 14 agents (molecules) and 3 strategies per player. Indeed, as players are assimilated to some main modelling notions, their number is small and often restricted to two antagonistic concepts like Physician and Disease (Section 3). Besides, by focusing on steady states for equilibrium scoring, we can compute steady states symbolically using SAT-solvers and modular network decomposition (Bollman et al., 2007; Corblin et al., 2009; Delaplace et al., 2012), which drastically reduce the computation time. Therefore for twoplayers using symbolic steady states computation, the complexity is in O(m2 · ı · ω(n) + m4 ), where ω stands for a function estimating the average complexity of the symbolic steady state computation (bounded by 2n ). 3. Network-Action Game applied to best drug selection in breast cancer In this section, we introduce the modelling principles for the application of Network-Action Game in Precision Medicine, we describe the construction of a Breast Cancer model and show the interpretation of the result of the application of Network-Action Game on this model. 3.1. Modelling principles 3.1.1. From formal definitions to biomedical interpretation Boolean networks have been used in Biology to study the dynamics of signalling and gene regulatory networks, in particular in the case of cancer (Calzone et al., 2010; Grieco et al., 2013; Lu et al., 2015; Von der Heyde et al., 2014). In this context, the population of agents A represents biomolecules (genes, RNA, proteins, metabolites, etc.), the family of formulas F defines the dynamical evolution of the activation/inhibition or presence/absence status of each molecule. Finally, the equilibria of Boolean networks are considered as representative of cell behaviour (Thomas, 1995). We consider a two-players Network-Action Game where the players are the Physician (p1 ) and the Disease (p2 ) whose respective strategies correspond to sets of drugs (Cp1 ) and genetic events (Cp2 ). Thus, a strategy profile c is a combination (Drug, Genetic Event). The action of strategies on the patient are translated as addition or deletion of arcs in the patient’s healthy interaction graph Gref representing its interactome network. We define cref = (, ) as the strategy profile where neither the Disease nor the Physician acts and leading to the healthy interaction graph Gref . Notice that the saturated Boolean network is only used for computational purposes, in order to determine each Fc resulting from the application of strategy profiles, and has no biological interpretation. Each strategy profile is thus associated to an interactome network whose dynamics is

3 The evaluation has been performed on a HP ZBook 15 Quad Core (4 2.8 GHz-Core i7 processors with 16 Go of memory).

C. Biane et al. / BioSystems 150 (2016) 52–60

modelled by a Boolean network. To assess the result of the strategy profiles actions, a scoring function pi : Bn → R is defined for each player on the state of some molecules (the markers) assimilated to biomarkers in Precision Medicine (Mirnezami et al., 2012; Frank and Hargreaves, 2003). By comparing the scores at the equilibria of the dynamics, we determine the preferences of the players on strategy profiles, and the interpretation of the obtained preference graph leads to the prediction of the efficiency of a drug, with respect to the genetic events. 3.1.2. Modelling choices and stages Breast Cancer (BC) has been extensively studied in molecular biomedicine, unravelling the signalling pathways and the molecular mechanism underlying the proliferation of cancerous cells and consequently enabling the development of targeted drugs impairing this proliferation (Hanahan and Weinberg, 2011; Campisi, 2003). In this section we describe the choices and stages of the construction of the BC model. Stage 1: construction of Gref . We constructed an interaction graph corresponding to a signalling network involved in the control of healthy breast cell proliferation and that is impaired in cancer cells. As proliferation of cells relies on the balance between apoptosis and division, the signalling network model must integrate: • The main proteins of signalling pathways responsible for the regulation of the apoptosis/division balance; • The effectors of the proliferation behaviour controlled by these pathways. The identification of these molecules was performed based on literature search about genetic events and targeted drugs in BC together with key molecular pathways involved in cell proliferation. After the identification of the main molecules (nodes that seems essential and at the crossroads of several pathways to enable extensions of the model), we added the activatory or inhibitory interactions (labelled arcs) in order to form linear pathways controling the effectors of apoptosis and division. The interactions were placed according to the KEGG database (Kanehisa and Goto, 2000) page describing the signalling pathways. When two molecules were interacting through a chain of intermediates, the sign of the interaction was chosen as follows: when the sum of inhibitory interaction is even, the interaction is positive, otherwise it is negative. We then confirmed each interaction by literature search. When the linear pathways were established, we checked whether these pathways were interconnected by looking on the KEGG page of each pathway if the molecules of one pathway were present in the KEGG page of some other pathways. Again, we searched for confirmation in the literature. Stage 2: definition of arc action function ϕ. Then, we modelled the structural perturbations of this signalling network induced by genetic events and targeted drugs as deletion or addition of arcs with respect to Gref . We seek in the literature the targets of drugs and genetic events and applied the following protocol: when a genetic event was described as leading to the deletion of the gene, its action was modelled as the deletion of all the output arcs of the node corresponding to the protein encoded in the target gene. Conversely, when a genetic event was described as leading to the over-expression of a protein, its action was modelled by the addition of an arc from a node the state of which is fixed to True to the protein encoded in the target gene. Finally, as the action of a drug corresponds to an inhibition of its molecular target (Higgins and Baselga, 2011; Sawyers, 2004), the action of the drug suppresses the action of its target and is modelled by deleting this target output arcs. Stage 3: strategy profiles CP : opposite actions and synthetic lethality. We observed that when the genetic event causing the

57

cancerous behaviour of the cell corresponds to over-expression of the gene (addition of an arc), the corresponding targeted drug has the opposite action: it inhibits the protein that is over-expressed (deletion of the added arc). Therefore, the action of a strategy profile (Drug, Over − expression) leads back to Gref . On the other hand, when the causal genetic event leads to the deletion of a gene and its protein (deletion of the output arcs), the action of the corresponding drug is more complex and relies on the synthetic lethality concept. Two genes are said synthetic lethal if a genetic event leading to the loss of one of them is viable but the mutations on both genes lead to cell death (William, 2005). The rationale behind the development of anticancer drug is to specifically trigger cell death in the cancer cell only. Therefore, the targeted drug is developed in order to inhibit the synthetic lethal partner of the mutated gene in the cancer cell (Gibbs, 2000). Hence the action of a strategy profile (Drug, Disease) leads to an interaction graph different from the healthy one Gref . Stage 4: Boolean network and score functions pi . For computational reasons, the model of dynamics is defined for the interaction graph GF of the saturated Boolean network representing the signalling network comprising all the possible interactions between molecules (healthy or added by strategies). We modelled the dynamics of this network by assigning logical formulas to each of the node of the saturated signalling network. When only one source influences the target, the establishment of the formula was straightforward. Otherwise, the formulas were derived from the literature. The nodes activated in case of DNA damage were set to True to model the fact that the cell is in stress, and the nodes that are over-expressed were also forced to True. The proliferation behaviour of cells depends on the balance between apoptosis and cell division and the efficiency of an anticancer drug depends on its ability to impair the proliferation of tumoral cells (Higgins and Baselga, 2011; Sawyers, 2004). Moreover, from a Darwinian point-of-view, we can assume that a tumoral cell has an “incentive” in proliferating (Campisi, 2003). Finally, cancer cells in quiescent and dormant states are responsible for relapses occurring many years after the treatment and healing the patient (Spiliotaki et al., 2012). We discriminate between this four proliferation behaviours (quiescence, dormancy, proliferation and death) by considering the states of two biomarkers: the effectors of apoptosis and division. When considering two markers, four states are possible: either they are both active or inactive, interpreted respectively as dormancy and quiescence (Spiliotaki et al., 2012), either one of them is active while the other is inactive (proliferation and death). We therefore defined opposite scores for both players: the worst for the disease is the best for the physician and vice versa and assigned to the player Disease, a maximal score when mitosis is active and apoptosis is inactive (the cancer cell is proliferating), a minimal score when mitosis is inactive and apoptosis is active (the cell is dying), and intermediate scores for both players when the cell is in a quiescence or dormancy. 3.2. Description of the breast cancer model In the following, we describe the molecules and interactions of the healthy signalling network, the actions of drugs and genetic events, the saturated Boolean network and the score functions for each player. 3.2.1. Healthy signalling network The genes affected in BC cells are involved in the cell division and the apoptosis control in response to stress through the PI3K/Akt, p53 and DNA repair pathways involving PARP and BRCA (Gasco et al., 2002; Sherr and McCormick, 2002; Engelman, 2009; Shaw and Cantley, 2006; Xu et al., 2001; Deng, 2006). These pathways collaborate to regulate the activation of Cyclin D1 and Bax, which

58

C. Biane et al. / BioSystems 150 (2016) 52–60 Table 3 Arc action function ϕ indicating the actions on arcs by strategies for each player. Player

Physician

+

HER2 → PI3K

ER ˛ → PI3K



0 −2 1 1

0 1 −2 1

0 2 1 1

0 1 2 1

Trastuzumab Tamoxifen Olaparib

 Disease

Fig. 6. Healthy signalling interaction graph Gref (arcs in black). The arcs in grey are present in the interaction graph of the saturated Boolean Networks but not in the healthy one.

are respectively the regulator of the G1/S transition during the cell division and a pro-apoptotic factor initiating apoptosis (Musgrove et al., 2011; Ward et al., 2008) and therefore considered as the markers of these two processes. The PI3K/Akt pathway is a phosphorylation cascade that promotes cell cycle progression through the inactivation of GSK3ˇ and prevents apoptosis through the activation of Bcl2, an inhibitor of Bax (Chang et al., 2003). The PI3K/Akt pathway interacts with p53 signalling through the activation of its inhibitor Mdm2 (Ute et al., 2003). In turn, p53 inhibits PI3K signalling through the activation of its inhibitor PTEN (Stambolic et al., 2001), therefore forming a circuit (Mayo and Donner, 2002; Carracedo and Pandolfi, 2008). p53 is also involved in the activation of the apoptosis through direct activation of Bax transcription (Toshiyuki and Reed, 1995; Basu and Haldar, 1998). BRCA increases Bax activation transcription through p53 activation (Xu et al., 2001; Mullan et al., 2006). BRCA is also involved in cell cycle arrest at the G1/S checkpoint (Mullan et al., 2006; Deng, 2006); this mechanism is modelled by an inhibition of CycD1 by BRCA. Finally, PARP inhibition induces cell cycle arrest and enhances cell death in a p53dependent manner (Nguyen et al., 2011). This is modelled as PARP activation of Cyclin D1 and PARP inhibition of p53. The healthy interaction graph Gref is shown in Fig. 6. 3.2.2. Actions of genetic events and drugs Three different subtypes of Breast cancer are defined based on the genetic events observed in the tumoral cells: ER positive, HER2 positive and BRCA-deficient breast cancer cells and are respectively treated with three different drugs: Tamoxifen, Trastuzumab and Olaparib (Higgins and Baselga, 2011). ER-positive BC cells over-express the gene coding the Oestrogen Receptor (ER) while HER2-positive BC cells over-express the gene coding for the Human Epidermal Receptor-2 (HER2) and BRCA-deficient BC cells are characterized by genetic events in both alleles of the BRCA gene leading to a deficiency of BRCA protein (Welcsh and King, 2001). We modelled these disease actions as additions of arcs between ER and HER2 nodes and their targets, and by the deletion of arcs connecting the BRCA node to its targets. The drugs inhibit specific molecular targets, respectively ER, HER2 and PARPs (Higgins and Baselga, 2011). Action arc function ϕ indicating addition or deletion of arcs by strategies for each players is shown in Table 3. The interaction graph of the saturated Boolean network includes arc in grey in Fig. 6. This graph represents all the interactions of healthy interaction graph and the interactions that can be added by strategies. Three interactions are added compared to the healthy interaction graph: the activation of PI3K by HER2 that occurs in HER2 positive breast cancer cells (Hynes and Lane, 2005; Yarden and Sliwkowski, 2001) and the activation of PI3K by ER˛ that occurs in presence of E2 in ER-positive breast cancer cells (Lee et al., 2005; Sun et al., 2001).

+

Strategy/arcs

HER2 ER BRCA +



Strategy/arcs

PARP → CycD1

PARP → p53



HER2 ER BRCA

1 1 1 −2 1 1 1 1

1 1 1 −2 1 1 1 1

Strategy/arcs

BRCA → p53

BRCA → CycD1



1 1 1 1 1 1 1 −2

1 1 1 1 1 1 1 −2

Trastuzumab Tamoxifen Olaparib



Trastuzumab Tamoxifen Olaparib

 HER2 ER BRCA

+



Table 4 Logical formulas underlying the activation of nodes. Protein

Logical formula

E2 HER2 BRCA PARP ER˛ PI3K Akt MDM2 GSK3ˇ p53 PTEN Bax Bcl2 CycD1

True True True True E2 ¬PTEN∧(ER˛ ∨ HER2) PI3K Akt ¬Akt ¬MDM2∧(BRCA∨ ¬PARP) p53 ¬Bcl2∧p53 Akt ¬GSK3ˇ ∨(¬BRCA∧PARP)

3.2.3. Boolean network and score function on cell proliferation biomarkers The facts driving the establishment of the logical formulas (given in Table 4) controling the activation of the nodes of the saturated Boolean network (Fig. 6) are the following: • Both ER˛ and HER2 can activate PI3K signalling (or) (Hynes and Lane, 2005; Lee et al., 2005; Sun et al., 2001; Yarden and Sliwkowski, 2001). • In a cell under DNA damage stress, both BRCA and PARP genes are activated (Livraghi and Garber, 2015), thus we set the nodes to True. • p53 activity is regulated in two ways by Mdm2: the binding of Mdm2 to p53 blocks p53 transcriptional activity, and Mdm2 exports p53 to the cytoplasm and targets it for proteasomal degradation (Mayo and Donner, 2002), we thus consider that the absence of Mdm2 is a necessary condition for p53 activation and. Moreover the similar role of PARP and BRCA in DNA repair leads to an or (William, 2005; Gibbs, 2000).

C. Biane et al. / BioSystems 150 (2016) 52–60 Table 5 Score function pi with respect to the biomarkers Boolean profile. CycD1

Bax

Cell behaviour

Physician score ( p1 )

Disease score (p2 )

0 1 0 1

1 0 0 1

Death Division Quiescence Dormancy

100% 0%

0% 100%

50%

50%

59

Ordinal Nash equilibria, indicating a situation where Disease has permanently disappeared in presence of these drugs. Therefore from the preference graph, we can conclude that Trastuzumab is efficient to heal HER cancer, Tamoxifen for ER cancer and Olaparib for BRCA cancer. These conclusions are confirmed by clinical practice since the associations described by these strategy profiles are currently used in the clinic (Higgins and Baselga, 2011). Hence, in this case, the network action game framework has inferred the best drug strategy selection for three types of different genetic events causing BC without explicit knowledge on these associations.

4. Conclusion

Fig. 7. Preference graph.

• Bax gene transcription is directly activated by p53 (Toshiyuki and Reed, 1995) and its translocation to its active site (mitochondria) is blocked by the PI3K/Akt pathway (Tsuruta et al., 2002). We thus consider that the absence of Akt is a necessary condition for Bax activation (and). • Cyclin D1 degradation is regulated both dependently and independently of GSK3 levels (or) (Alao, 2007). Moreover, as PARP and BRCA have the opposite action on DNA repair than they have on division, we modelled their action on Cyclin D1 by the negation of their action on p53. We define the score functions pi for each player pi on the states of Cyclin D1 and Bax, which are respectively the regulator of the G1/S transition of the cell cycle and a pro-apoptotic factor initiation apoptosis and as a consequence can be considered as biomarkers of division and apoptosis respectively (Musgrove et al., 2011; Ward et al., 2008). The scores for each player are given in Table 5. 3.3. Result and interpretation of the application of Network-Action Game on the BC model The comparison of the scores, for each player, at the equilibria of the dynamics (with Eq. (3)) determine the preference relations on strategy profiles ( pi ) for each player. Fig. 7 shows the preference graph representing the union of the preference relations obtained from the application of Network-Action Game on the BC model. On this graph, the nodes represent strategy profiles (Drug, Genetic Event) and the arcs represent Physician preferences (black) and Disease preferences (grey). On the preference graph, one can observe that Disease does not prefer a strategy profile corresponding to a cancer treated with a drug (Tra, HER; Ola, BRCA; Ta, ER; highlighted in grey) to a strategy profile corresponding to its absence of action (Tra, ; Ola, ; Ta,  respectively). Disease is unable to discriminate between a healthy and a treated state with a Genetic Event in these situations and we therefore interpret them as the efficiency of the drug. Morevover, Physician does not change his strategy while on (Tra, HER) and (Ta, ER) or (Ola, BRCA). Also notice that (Tra, ), (Ola, ) and (Ta, ) are

Network action game couples Boolean networks with Game Theory in order to predict the best association between drug and disease. Arc action on the interactome is considered here as a paradigm of the causal explanation of disease and therapy prediction. The efficiency of the framework has been assessed on a breast cancer model. As a result we show that the proposed associations between drugs and malignant genetic events exactly match to those found in literature. The application of Network-Action Game requires to first, interpret genetic events and drugs mode-of-action in terms of network actions on a Boolean network, second, to determine the markers discriminating the phenotypic states and finally, to assign the scores based on a qualitative assessment on the objectives of the physician and the disease. In the application of the method to breast cancer these elements are derived from literature and databases. Therefore, the efficiency of the method relies on the availability and quality of biological knowledge for the studied case. Future prospects mainly concern the application of the network action game to drug re-purposing. The former investigates the repositioning of drugs to new indications leading to substantially reduce the duration and the cost of their development cycle since the necessary analysis of the molecules was already performed. This may be addressed by generalizing the method applied for breast cancer. Indeed, screening in-silico drug actions on arcs against disease ones may be realized by assessing the consequence on the dynamics from the evaluation of marker states. To be feasible, this approach requires to automatically characterize actions from data and knowledge on drugs and to get reliable and complete description of network encompassing all the actions on arcs. The computational challenge for drug action discovery is the inverse problem of the therapy prediction, where the effects described by states of markers are known but the causes defined as actions must be discovered. Hence the issue is to infer the necessary actions on a diseased network in order to bring the dynamics back to the healthy state.

References Abou-Jaoudé, W., Ouattara, D.A., Kaufman, M., 2009. From structure to dynamics: frequency tuning in the p53-Mdm2 network. I. Logical approach. J. Theor. Biol. 258 (4), 561–577. Alao, J.P., 2007. The regulation of cyclin d1 degradation: roles in cancer development and the potential for therapeutic invention. Mol. Cancer 6 (1), 1. Ballas, S.K., 2011. Defining the phenotypes of sickle cell disease. Hemoglobin 35, 511–519. Barabási, A.L., Gulbahce, N., Loscalzo, J., 2011. Network medicine: a network-based approach to human disease. Nat. Rev. Genet. 12 (1), 56–68. Barabási, A.L., 2007. Network medicine – from obesity to the “diseasome”. N. Engl. J. Med. 357 (4), 404–407. Basu, A., Haldar, S., 1998. The relationship between bci2, bax and p53: consequences for cell cycle progression and cell death. Mol. Hum. Reprod. 4 (12), 1099–1109. Bollman, D., Colón-Reyes, O., Orozco, E., 2007. Fixed points in discrete models for regulatory genetic networks. EURASIP J. Bioinform. Syst. Biol. 2007, 10.

60

C. Biane et al. / BioSystems 150 (2016) 52–60

Calzone, L., Tournier, L., Fourquet, S., Thieffry, D., Zhivotovsky, B., Barillot, E., Zinovyev, A., 2010. Mathematical modelling of cell-fate decision in response to death receptor engagement. PLoS Comput. Biol. 6 (3), e1000702. Campisi, J., 2003. Cancer and ageing: rival demons? Nat. Rev. Cancer 3 (5), 339–349. Carracedo, A., Pandolfi, P.P., 2008. The pten-pi3k pathway: of feedbacks and cross-talks. Oncogene 27 (41), 5527–5541. Chang, F., Lee, J.T., Navolanic, P.M., Steelman, L.S., Shelton, J.G., Blalock, W.L., Franklin, R.A., McCubrey, J.A., 2003. Involvement of pi3k/akt pathway in cell cycle progression, apoptosis, and neoplastic transformation: a target for cancer chemotherapy. Leukemia 17 (3), 590–603. Chettaoui, C., Delaplace, F., Manceny, M., Malo, M., 2007. Games network and application to PAS system. Biosystems 87 (2), 136–141. Ciliberto, A., Novak, B., Tyson, J.J., 2005. Steady states and oscillations in the p53/Mdm2 network. Cell Cycle 4 (3), 488–493. Corblin, F., Tripodi, S., Fanchon, E., Ropers, D., Trilling, L., 2009. A declarative constraint-based method for analyzing discrete genetic regulatory networks. Biosystems 98 (2), 91–104. Davidson, E.H., 2010. The Regulatory Genome: Gene Regulatory Networks in Development and Evolution. Academic Press. Delaplace, F., Klaudel, H., Cartier-Michaud, A.,2010. Discrete causal model view of biological networks. In: Proceedings of the 8th International Conference on Computational Methods in Systems Biology – CMSB’10. ACM Press, New York, NY, USA, pp. 4–13. Delaplace, F., Klaudel, H., Melliti, T., Sené, S.,2012. Analysis of modular organisation of interaction networks based on asymptotic dynamics. In: Computational Methods in Systems Biology. Springer, pp. 148–165. Deng, C.-X., 2006. Brca1: cell cycle checkpoint, genetic instability, DNA damage response and cancer evolution. Nucleic Acids Res. 34 (5), 1416–1426. Engelman, J.A., 2009. Targeting pi3k signalling in cancer: opportunities, challenges and limitations. Nat. Rev. Cancer 9 (8), 550–562. Frank, R., Hargreaves, R., 2003. Clinical biomarkers in drug discovery and development. Nat. Rev. Drug Discov. 2 (7), 566–580. Gasco, M., Shami, S., Crook, T., 2002. The p53 pathway in breast cancer. Breast Cancer Res. – BCR 4 (2), 70–76. Gibbs, J.B., 2000. Mechanism-based target identification and drug discovery in cancer research. Science (New York, NY) 287 (5460), 1969–1973. Ginsburg, G.S., 2009. Translational genomics: from discovery to clinical practice. In: Willard, H.F., Ginsburg, G.S. (Eds.), Genomic and Personalized Medicine. Academic Press, New York, pp. 262–274, Chapter 22. Grieco, L., Calzone, L., Bernard-Pierrot, I., Radvanyi, F., Kahn-Perles, B., Thieffry, D., 2013. Integrative modelling of the influence of MAPK network on cancer cell fate decision. PLoS Comput. Biol. 9 (10), e1003286. Gustafsson, M., Nestor, C.E., Zhang, H., Barabási, A.L., Baranzini, S., Brunak, S., Chung, K., Federoff, H.J., Gavin, A.-C., Meehan, R.R., Picotti, P., Pujana, M., Rajewsky, N., Smith, K., Sterk, P.J., Villoslada, P., Benson, M., 2014. Modules, networks and systems medicine for understanding disease and aiding diagnosis. Genome Med. 6 (10), 82. Hanahan, D., Weinberg, R.A., 2011. Hallmarks of cancer: the next generation. Cell 144 (5), 646–674. Higgins, M.J., Baselga, J., 2011. Targeted therapies for breast cancer. J. Clin. Invest. 121 (10), 3797–3803. Hynes, N.E., Lane, H.A., 2005. ERBB receptors and cancer: the complexity of targeted inhibitors. Nat. Rev. Cancer 5 (5), 341–354. Kanehisa, M., Goto, S., 2000. Kegg: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28 (1), 27–30. Le Roux, S., Lescanne, P., Vestergaard, R., 2008. Conversion/Preference Games. CoRR, abs/0811.0071. Lee, Y.-R., Park, J., Yu, H.-N., Kim, J.-S., Youn, H.J., Jung, S.H., 2005. Up-regulation of pi3k/akt signaling by 17ˇ-estradiol through activation of estrogen receptor-˛, but not estrogen receptor-ˇ, and stimulates cell growth in breast cancer cells. Biochem. Biophys. Res. Commun. 336 (4), 1221–1226. Livraghi, L., Garber, J.E., 2015. Parp inhibitors in the management of breast cancer: current data and future prospects. BMC Med. 13 (1), 1. Lu, J., Zeng, H., Liang, Z., Chen, L., Zhang, L., Zhang, H., Liu, H., Jiang, H., Shen, B., Huang, M., et al., 2015. Network modelling reveals the mechanism underlying colitis-associated colon cancer and identifies novel combinatorial anti-cancer targets. Sci. Rep. 5. Mayo, L.D., Donner, D.B., 2002. The PTEN, Mdm2, p53 tumor suppressor-oncoprotein network. Trends Biochem. Sci. 27 (9), 462–467. Milo, R., Shen-Orr, S., Itzkovitz, S., Kashtan, N., Chklovskii, D., Alon, U., 2002. Network motifs: simple building blocks of complex networks. Science 298 (5594), 824–827.

Mirnezami, R., Nicholson, J., Darzi, A., 2012. Preparing for precision medicine. N. Engl. J. Med. 366, 489–491. Mullan, P.B., Quinn, J.E., Harkin, D.P., 2006. The role of brca1 in transcriptional regulation and cell cycle control. Oncogene 25 (43), 5854–5863. Musgrove, E.A., Caldon, E.C., Barraclough, J., Stone, A., Sutherland, R.L., 2011. Cyclin d as a therapeutic target in cancer. Nat. Rev. Cancer 11 (8), 558–572. Nguyen, D., Zajac-Kaye, M., Rubinstein, L., Voeller, D., Tomaszewski, J.E., Kummar, S., Chen, A.P., Pommier, Y., Doroshow, J.H., Yang, S.X., 2011. Poly (adp-ribose) polymerase inhibition enhances p53-dependent and-independent DNA damage responses induced by DNA damaging agent. Cell Cycle 10 (23), 4074–4082. Osborne, M.J., Rubinstein, A., 1994. A Course in Game Theory. MIT Press. Oti, M., Snel, B., Huynen, M.A., Brunner, H.G., 2006. Predicting disease genes using protein–protein interactions. J. Med. Genet. 43 (8), 691–698. Sawyers, C., 2004. Targeted cancer therapy. Nature 432 (7015), 294–297. Schadt, E.E., 2009. Molecular networks as sensors and drivers of common human diseases. Nature 461 (September), 218–223. Sebastiani, P., Ramoni, M.F., Nolan, V., Baldwin, C.T., Martin, H., 2010. Genetic dissection and prognostic modeling of overt stroke in sickle cell anemia. Nat. Genet. 37 (4), 435–440. Shaw, R.J., Cantley, L.C., 2006. Ras, pi (3) k and mtor signalling controls tumour cell growth. Nature 441 (7092), 424–430. Sherr, C.J., McCormick, F., 2002. The rb and p53 pathways in cancer. Cancer Cell 2 (2), 103–112. Spiliotaki, M., Markomanolaki, H., Kallergi, G., Papadaki, M.A., Georgoulias, V., Mavroudis, D., Agelaki, S., 2012. Evaluation of proliferation and apoptosis markers in circulating tumor cells (CTCs) of women with early breast cancer who are candidates for tumor dormancy. Cancer Res. 72 (8 Supplement), 2397–2397. Stambolic, V., MacPherson, D., Sas, D., Lin, Y., Snow, B., Jang, Y., Benchimol, S., Mak, T.W., 2001. Regulation of pten transcription by p53. Mol. Cell 8 (2), 317–325. Sun, M., Paciga, J.E., Feldman, R.I., Yuan, Z.-q., Coppola, D., Lu, Y.Y., Shelley, S.A., Nicosia, S.V., Cheng, J.Q., 2001. Phosphatidylinositol-3-oh kinase (pi3k)/akt2, activated in breast cancer, regulates and is induced by estrogen receptor ˛ (er˛) via interaction between er˛ and pi3k. Cancer Res. 61 (16), 5985–5991. Tarjan, R., 1972. Depth-first search and linear graph algorithms. SIAM J. Comput. 1 (2), 146–160. Thomas, R., Thieffry, D., 1995. Dynamical behaviours of regulatory networks – II. immunity control in bactériophage lambda. Bull. Math. Biol. 57 (2), 277–297. Thomas, R., 1973. Boolean formalization of genetic control circuits. J. Theor. Biol. 42 (3), 563–585. Miyashita, T., Reed, J.C., 1995. Tumor suppressor p53 is a direct transcriptional activator of the human bax gene. Cell 80 (2), 293–299. Tsuruta, F., Masuyama, N., Gotoh, Y., 2002. The phosphatidylinositol 3-kinase (pi3k)-akt pathway suppresses bax translocation to mitochondria. J. Biol. Chem. 277 (16), 14040–14047. Moll, U.M., Petrenko, O., 2003. The mdm2-p53 interaction. Mol. Cancer Res. 1 (14), 1001–1008. Von der Heyde, S., Bender, C., Henjes, F., Sonntag, J., Korf, U., Beißbarth, T., 2014. Boolean erbb network reconstructions and perturbation simulations reveal individual drug response in different breast cancer cell lines. BMC Syst. Biol. 8 (1), 1. Ward, T.H., Cummings, J., Dean, E., Greystoke, A., Ho, J.-M., Backen, A., Ranson, M., Dive, C., 2008. Biomarkers of apoptosis. Br. J. Cancer 99 (6), 841–846. Welcsh, P.L., King, M.-C., 2001. Brca1 and brca2 and the genetics of breast and ovarian cancer. Hum. Mol. Genet. 10 (7), 705–713. Kaelin, W.G., 2005. The concept of synthetic lethality in the context of anticancer therapy. Nat. Rev. Cancer 5 (9), 689–698. Xu, X., Qiao, W., Linke, S.P., Cao, L., Li, W.-M., Furth, P.A., Harris, C.C., Deng, C.-X., 2001. Genetic interactions between tumor suppressors brca1 and p53 in apoptosis, cell cycle and tumorigenesis. Nat. Genet. 28 (3), 266–271. Yarden, Y., Sliwkowski, M.X., 2001. Untangling the ErbB signalling network. Nat. Rev. Mol. Cell Biol. 2 (2), 127–137. Zhong, Q., Simonis, N., Li, Q.-R., Charloteaux, B., Heuze, F., Klitgord, N., Tam, S., Yu, H., Venkatesan, K., Mou, D., Swearingen, V., Yildirim, M.A., Yan, H., Dricot, A., Szeto, D., Lin, C., Hao, T., Fan, C., Milstein, S., Dupuy, D., Brasseur, R., Hill, D.E., Cusick, M.E., Vidal, M., 2009. Edgetic perturbation models of human inherited disorders. Mol. Syst. Biol. 5 (321), 321.

Networks and games for precision medicine.

Recent advances in omics technologies provide the leverage for the emergence of precision medicine that aims at personalizing therapy to patient. In t...
1MB Sizes 0 Downloads 10 Views