Systemic modeling myeloma-osteoclast interactions under normoxic/hypoxic condition using a novel computational approach

Interaction of myeloma cells with osteoclasts (OC) can enhance tumor cell expansion through activation of complex signaling transduction networks. Both cells reside in the bone marrow, a hypoxic niche. How OC-myeloma interaction in a hypoxic environment affects myeloma cell growth and their response to drug treatment is poorly understood. In this study, we i) cultured myeloma cells in the presence/absence of OCs under normoxia and hypoxia conditions and did protein profiling analysis using reverse phase protein array; ii) computationally developed an Integer Linear Programming approach to infer OC-mediated myeloma cell-specific signaling pathways under normoxic and hypoxic conditions. Our modeling analysis indicated that in the presence OCs, (1) cell growth-associated signaling pathways, PI3K/AKT and MEK/ERK, were activated and apoptotic regulatory proteins, BAX and BIM, down-regulated under normoxic condition; (2) β1 Integrin/FAK signaling pathway was activated in myeloma cells under hypoxic condition. Simulation of drug treatment effects by perturbing the inferred cell-specific pathways showed that targeting myeloma cells with the combination of PI3K and integrin inhibitors potentially (1) inhibited cell proliferation by reducing the expression/activation of NF-κB, S6, c-Myc, and c-Jun under normoxic condition; (2) blocked myeloma cell migration and invasion by reducing the expression of FAK and PKC under hypoxic condition.


Supplementary files list:
Text S1 State consistent rules.
Text S2 Computational procedure: Inference of cell-specific pathways with DILP.
Text S3 Computational procedure: fitting precision of data (goodness of fit).
Text S4 Computational procedure: addition of missing edges into a single optimal solution.
Text S5 Computational procedure: prediction of drug treatment effects by state transition analysis.

Table S1
The fitting precision of the proposed DILP model on the processed RPPA data.

Table S2
The expressions of all the proteins involved in the OC-mediated MM-specific pathways in normoxia (with the presence of OC) before and after treatment with PI3K inhibitor.

Table S3
The expressions of all the proteins involved in the OC-mediated MM-specific pathways in hypoxia (with the presence of OC) before and after treatment with combination of PI3K and integrin inhibitors.

Table S4
The details of the antibodies.

S1. State consistent rules
In our study, signaling pathway network was represented as a set of signaling proteins = , , , , , + + , − * , + 3 2 , ≤ 3 (8) , To infer cell-specific pathway networks, we applied our DILP approach with above constraints to minimize the differences between experimentally measured and predicted values of signaling proteins, as well as to obtain a minimized sub-network of original generic pathway map by optimizing the following objective function shown in formula (25): In this objective function, the first term represents the fitting error between experimentally measured and predicted values and the second term denotes a set of reactions from the original generic signaling network. The measured and predicted values of -th protein at time point were denoted as , , , ∈ {−1,0,1}, respectively. When , and , are equal, the value of term ( , − , ) 2 will be 0; otherwise it is either 1 or 4. Hence, optimization of the above objective function might induce localoptimal solution because of the non-uniform distribution of the term ( , − , ) 2 . In order to address this bias, binary variable , (taking value of 0 or 1) was designed as the difference between , and , as following constraints (26-27): , will be 1 if , is not equal to , ; the minimum of formula (28) will automatically set , as 0 if , is equal to , .
Then the term ( , − , ) 2 in above objective function was replaced by , .Therefore, the objective function (25) can be simplified to formula (28).
The negative constant in formula (28) is used to obtain a minimum sub-graph of the generic pathway maps as the finalized cell-specific pathways (here, we have− 1 | | < < 0), in which | | is the number of reactions in the network.
In addition, minimizing the optimal network topology by edge removals might eliminate some reactions, leading to some phosphor-signals can't be transduced into downstream proteins (such as two examples shown in Supplementary Fig. S6 and S7, respectively). After obtaining a single minimized subgraph of generic pathway network via DILP approach, we designed a strategy for searching the missing edges. The missing edges were retrieved to the optimal network obtained from formula (28) Where the binary variable , indicates the difference between the measurement , and predicted value , of -th protein at the time point . and are the total number of measured proteins and time points, respectively. The value of fitting precision ( ) is in the range from 0% to 100%.

S4. Computational procedure: addition of missing edges into a single optimal solution
Minimizing the optimal network topology by edge removals might eliminate some reactions, leading to some of phosphor-signals can't be transduced from upstream to downstream. Here, we firstly represent Supplementary Fig. S6-S7 as two examples to describe the situation which might occur in the optimization.
Supplementary Fig. S6(A-B) represent a generic pathway map and the inferred cell-specific pathway map, respectively. The states of and were measured as up-regulation, and was unknown (unmeasured). The goodness of fit reaches to 100% after optimization; however, minimizing the optimal network topology by edge removals eliminate the reactions 1 and 2 , which lead to the phosphor-signal from can't be transduced into downstream proteins . The predicted states of and in the network shown in Supplementary Fig S6(B) were up-regulated ("1") and was up-changed ("0"). The reasonable inferred cell-specific pathway network should be as shown in Supplementary Fig. S6(C).
In Supplementary Fig. S7(A), the downstream protein is activated if and only if at least one of three upstream proteins , , and is activated. Given a scenario that , , and were measured as up-regulation, and is unknown, the optimization of formula (28) might deliver one of two candidate optimal solutions shown in Supplementary Fig. S7(B-C). Supplementary Fig. S7(B-C) indicate that two candidate optimal network topology which both have the best fitting precision (100%) to measurements. However, minimization of network topology in this case also result in one upstream real effect cannot be transduced to downstream. Therefore, we consider the solution shown in Supplementary Fig. S7(D) might be more reasonable which consistent with the measurement. As to the node , its predicted state can be one of three possible states (up-regulation, down-regulation, and nochange) because un-measured nodes would not be included in the objective function (28).
In this study, we proposed a method to detect the missing edges whose addition would not reduce the goodness of fit. After computing a single minimal sub-graph of generic pathway network, we consider a greedy strategy for searching potential missing edges and keeping the same goodness of fit as that in the first run of formula (28). The pseudo-code for adding the missing edges was described as follows: Step 1: Get an initial optimal solution ( 1 ) from formula (28) with a negative (minimize the topology).
Step 2: Based on 1 : For each reaction ∈ , At the time point , , and ℎ , are the predicted values of parent and child node involved in reaction , respectively. Here, we still took Supplementary Fig. S7(A) as an example to illustrate our greedy search strategy for searching and adding the missing edges. At first, we might get a solution such as Supplementary Fig.   S7(B) and the edge → was a missing edge. After step 2, an extra constraint set for → was added and went into step 3. And then we obtained another optimal solution shown in Supplementary   Fig. S7(C) which also induced the same fitting precision as the solution in Fig. S7(B). Hence, we consider two cases in Figure S7(B-C) are both the candidate solutions; however, they both missed an edge. In step 4, our approach found the link from → can be removed because is un-measured and a minimal sub-graph of the generic pathway network is required. After step 4, a subset of reactions was obtained which eliminated the link → . Last step, we consider the required optimal solution will restrict in the scope of , hence, we run the optimization with formula (28) again to search a maximal topology with the same fitting precision and covered the missing ages. Finally, we got an adjusted solution shown in Supplementary Fig. S7(D), which also have the same goodness of fit (100%) as Supplementary Fig. S7 (B-C). Similarity, Supplementary Fig. S6(C) represents that the addition of two reactions ( 1 and 2 ) still keep the goodness of fit as 100%.

S5. Computational procedure: prediction of drug treatment effects by state transition analysis
Based on the established network topological structure and transfer functions (Boolean operation), state transition analysis in Boolean networks is a kind of approach to predicting the future state of each node from the current states of its parental nodes [3][4][5] . The assumption of state transition is that: the state of protein at time point + 1 ( , +1 ) is associated with the states of its parental proteins at time point 5,6 . Thus, we predicted the state of pathway network at time point + 1 from the state at time point through state transition, following perturbation of the cell-specific pathway network with drugs.
Given the inferred cell-specific pathway network topology and the states of all the proteins involved in the network at time point : = [ 1, , 2, , … , , ] (n is the total number of proteins in the network), we can obtain the states of these proteins at time point + 1 using the following formula: where = 1,2, … , , and , {−1,0,1}. is a set of transfer functions to change the signaling network from one state to another. In Boolean networks, transfer functions are denoted by using logical expressions via Boolean operators 5,7 . However, it is difficult to represent the transfer functions using an established mathematic expression, if each signaling protein has three possible states (1, -1, or 0) 2 . In our study, the transfer functions were represented by a set of integer linear constraints.
Let's use 0 = { 1,0 , 2,0 , … , ,0 , … , ,0 } ( ,0 ∈ {−1,1,0}) to denote a measured state of signaling network without any intervention or treatment. When cells are treated by an inhibitor, we assume protein is targeted by this inhibitor, and then the state of above signaling network after treatment is changed as 0 = { 1,0 , 2,0 , … , −1, … , ,0 }. We used 0 as the initial state to predict the performance of the inhibition or perturbation. We then used formula (30) to calculate the next state ( 1 ) of 0 and repeat this process until the signaling network reached to steady state 7 . We eventually generated a set of states { 1 , 2 ,…, }. is regarded as the finally effects of drug treatment in cell-specific pathway network.
In here, the constraints involved in the transfer functions were represented as following constraints: , → +1 + + 3 1 , ≤ 3 , → +1 + − * , + 3 2 , ≤ 3 (34)          Table S2. The expressions of all the proteins involved in the OC-mediated MM-specific pathways in normoxia (with the presence of OC) before and after treatment with PI3K inhibitor. The expressions of all the proteins before treatment (the second column in Table S2) were collected (at 24h) after we inferred the MM-specific pathways by fitting the measured data with the generic pathway map (see Fig.  S4). The expressions of all the proteins after treatment (the third column in Table S2) were predicted by the state transition analysis with our DILP approach. In normoxia, treatment on myeloma cells with Pi3K inhibitor induced the changes of some key proteins (orange color).  Table S3. The expressions of all the proteins involved in the OC-mediated MM-specific pathways in hypoxia (with the presence of OC) before and after treatment with combination of PI3K and integrin inhibitors. The expressions of all the proteins before treatment (the second column in Table S3) were collected (at 24h) after we inferred the MM-specific pathways by fitting the measured data with the generic pathway map (see Fig. S5). The expressions of all the proteins after treatment (the third column in Table S3) were predicted by the state transition analysis with our DILP approach. In hypoxia, treatment on Myeloma cells with Pi3K inhibitor induced the changes of some key proteins (orange color).