Fuzzy modeling and global optimization to predict novel therapeutic targets in cancer cells.

MOTIVATION: The elucidation of dysfunctional cellular processes that can induce the onset of a disease is a challenging issue from both the experimental and computational perspectives. Here we introduce a novel computational method based on the coupling between fuzzy logic modeling and a global optimization algorithm, whose aims are to (1) predict the emergent dynamical behaviors of highly heterogeneous systems in unperturbed and perturbed conditions, regardless of the availability of quantitative parameters, and (2) determine a minimal set of system components whose perturbation can lead to a desired system response, therefore facilitating the design of a more appropriate experimental strategy.
RESULTS: We applied this method to investigate what drives K-ras-induced cancer cells, displaying the typical Warburg effect, to death or survival upon progressive glucose depletion. The optimization analysis allowed to identify new combinations of stimuli that maximize pro-apoptotic processes. Namely, our results provide different evidences of an important protective role for protein kinase A in cancer cells under several cellular stress conditions mimicking tumor behavior. The predictive power of this method could facilitate the assessment of the response of other complex heterogeneous systems to drugs or mutations in fields as medicine and pharmacology, therefore paving the way for the development of novel therapeutic treatments.
AVAILABILITY AND IMPLEMENTATION: The source code of FUMOSO is available under the GPL 2.0 license on GitHub at the following URL: https://github.com/aresio/FUMOSO. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

1 Introduction

Cells are complex heterogeneous systems, whose functioning is governed by a finely regulated interplay between various types of molecules involved in gene expression, signal transduction and metabolic pathways, altogether resulting in different cellular phenotypes. Dysfunctional processes caused by events occurring at the molecular level can induce a cascade of local and global damages in cells, tissues, organs and, possibly, in the whole organism. Therefore, understanding molecular regulations at a mechanistic level is indispensable either to prevent or control the onset of many diseases. In this context, the integration between experimental data and computational methods facilitate the definition of predictive mathematical models, whose simulations can elucidate the emergent properties of the biological system in physiological and pathological conditions, reveal possible counter-intuitive mechanisms and envisage new hypotheses that can be tested in the laboratory (Faeder and Morel, 2016; Kitano, 2002).

The mathematical description of biological systems can be realized with different approaches, such as mechanism-based (Wilkinson, 2009) or logic-based modeling (Le Novère, 2015). Mechanism-based models provide a detailed description of the underlying biochemical reactions (Chylek ; Szallasi ). However, these models require the knowledge or the inference (Chou and Voit, 2009) of quantitative parameters (e.g. kinetic rates, molecular amounts) that are often difficult to be measured, especially in vivo and for large-scale systems, therefore hampering the effectiveness of many computational analyses (Somogyi ; Tangherloni ). Moreover, biological systems and their components are often described in a qualitative way, by using natural language terms, such as ‘moderately active’ or ‘highly expressed’, which highlight biological uncertainty and experimental measurement limitations. Given these issues, logic-based models may be used as a reliable alternative to study cellular systems (Morris ; Samaga and Klamt, 2013; Wynn ), since they are more suitable than mechanistic models when only qualitative data are available (Flobak ; Fumiã and Martins, 2013; Zañudo and Albert, 2015). In particular, fuzzy logic represents a powerful extension of Boolean logic to model complex systems, since it generalizes the binary formalization of the variables states to deal with any uncertainty related to the system (Yen and Langari, 1999). Fuzzy logic-based models are defined by: (i) a directed graph representing the set of system components (i.e. linguistic variables) and their mutual positive or negative regulations; (ii) a set of linguistic terms for each variable, necessary to give a qualitative description of intracellular concentrations or functional activities. Noteworthy, linguistic terms bypass the necessity of a precise parameterization but, at the same time, they allow to provide a quantitative representation of the variables states thanks to the so-called membership functions; (iii) a set of fuzzy logic rules that specify the state that each linguistic variable will assume over time, according to the states of the variables by which it is regulated.

So far, fuzzy logic has been applied to different strategies in the field of cellular biology, for instance, to model and simulate signaling pathways (Aldridge ) or gene regulatory networks (Küffner ), overcome the lack of kinetic parameters (Bordon ; Liu ), automatically infer network models (Keller ; Liu ; Morris ; 2016) or implement regression models (Schmidt-Heck ). In this work, we propose a novel and general-purpose computational method, based on fuzzy logic, designed to facilitate the modeling and analysis of heterogeneous systems. In particular, our method was designed to provide a mathematical description of complex biological systems whose components range from single molecules to whole cellular processes and observable cell phenotypes, together with their mutual regulations. One of the main advantage of this method is that, although it does not require the availability of quantitative (kinetic) parameters and exact values of the state or the abundance of cellular components, it is a dynamical modeling approach that allows to simulate and predict the temporal evolution of the system in both unperturbed and perturbed conditions. We also show that our fuzzy modeling approach, coupled with an optimization algorithm, automatically identifies a potential (minimal) set of system components whose perturbation can maximize, or minimize, a desired system response. This automatic identification represents one of the main novelties of our computational approach, which can largely facilitate the design of new laboratory experiments by yielding putative perturbations able to drive the behavior of an arbitrary complex system. Although several fuzzy logic tools and libraries are available in the literature, none of them was specifically designed to support the modeling, the dynamical simulation and the optimization of the heterogeneous systems that we aim to investigate. For this reason, our methodology was implemented from scratch and, in particular, we developed a novel user-friendly software named FUMOSO (FUzzy MOdel SimulatOr), which supports the definition, editing, export and simulation of heterogeneous fuzzy models of complex dynamical systems.

To show the potentiality of this novel computational method, we investigated a complex, heterogeneous system consisting of oncogenic K-ras cancer cells—characterized by the so-called ‘Warburg effect’—grown in a progressive limiting amount of glucose, in order to understand the glucose-dependent mechanisms driving cancer cells to death or survival. The Warburg effect, or aerobic glycolysis, is a metabolic hallmark of malignancy. Accordingly, numerous cancer cells, grown either in low glucose availability or in free glucose, are strongly susceptible to cell death as compared to normal cells. However, it has also been observed that not all cancer cells undergo cell death upon a really harsh environment, such as in glucose starvation, since some of them might acquire the ability to survive in this new environmental condition by activating compensatory signaling pathways (Huang ; Palorini ) and alternative metabolic routes (Ye ; Zaugg ). Worthy of note, metabolic rewired cancer cells, which often are more aggressive (Endo ), can be selected by chemotherapy, by therapies exploiting synergism between chemotherapeutic treatments and anti-metabolic drugs or by genetic and pharmacological ablation of oncogenic pathways (Elgendy ; Viale ; Zhao ). In addition, strategies to directly inhibit glycolysis in cancer patients have been partially dumped since they might also likely damage normal healthy tissues (i.e. smooth and skeletal muscle, and normal viscera). In this scenario, the combination of therapies targeting aerobic glycolysis, adaptive mechanisms (i.e. increased autophagy) and well-established cancer-specific targets (i.e. tyrosine kinase signaling pathways) represent a potential approach to be explored in cancer cure.

Here we present a fuzzy logic model of programmed cell death and survival under progressive glucose depletion, defined on the basis of an extensive prior knowledge of the main components involved in K-ras-transformed cells grown in this perturbed condition. In particular, we considered cancer cell death occurring upon glucose starvation along two major pathways: (i) a pathway centered on mitochondria (reactive oxygen species (ROS), adenosine triphosphate (ATP) depletion, calcium (Ca2+) overloading) (Elmore, 2007; Taylor ) and (ii) an endoplasmic reticulum (ER)-stress pathway associated with reduction of N-glycosylation and cell attachment, and a consequent activation of the unfolded protein response (UPR) leading to cell death (Hetz, 2012; Hetz and Papa, 2018). In contrast, we indicated two major mechanisms as survival routes: (i) mitochondrial activity rewiring and (ii) autophagy. The model was validated against data obtained from mouse fibroblasts transformed by oncogenic K-ras expression (NIH3T3 K-ras cells) and a human K-ras-mutated breast cancer cell line (MDA-MB-231), both grown in unperturbed and different perturbed conditions. The optimization analysis allowed to automatically search for and detect the combination of perturbations that maximize pro-apoptotic processes in cancer cells for the purpose of guiding the development of novel therapeutic treatments.

2 Materials and methods

2.1 Dynamic fuzzy rules-based modeling

A dynamic fuzzy rules-based model (DFM) is a computational paradigm to describe and analyze the emergent behavior of heterogeneous complex systems characterized by uncertainty. In DFMs, a linguistic variable and a set of linguistic terms (e.g. Low, Medium and High) are associated with each component of the system to provide a qualitative description of all the possible states that component can assume in time (Aldridge ; Yen and Langari, 1999). A DFM handles the intrinsic uncertainty of the state of the variables by means of the membership functions associated with each linguistic term. Linguistic variables and terms are used to define a set of fuzzy rules, which provide a qualitative description of the mechanisms (e.g. feedback regulation) driving the overall behavior of the system. Fuzzy rules are conditional statements generally written in the form ‘IF x IS a THEN y IS b’. The antecedent of a fuzzy rule is a predicate involving variables (i.e. a system component x) and their associated linguistic terms (i.e. a). The consequent of a fuzzy rule can either be a fuzzy set, a constant or a function (i.e. b) that is assigned to an output variable (i.e. y), according to the specific fuzzy inference engine employed.

The variables belonging to a DFM can be partitioned into two sets: outer and inner variables. The set of outer variables contains input and output variables, which can only appear as antecedents and consequents of fuzzy rules, respectively. Namely, input variables correspond to the components that trigger the dynamic evolution of the system, while output variables represent the components of interest for the analysis of the system (i.e. some experimentally measurable component). On the contrary, inner variables can appear on both sides of fuzzy rules, and they are used to represent mutual regulations among the system components.

The state of input variables is set (or perturbed) over time using appropriate user-defined functions, which induce the evolution in time of the whole system. On the contrary, the state of all other variables change as a result of the synchronous application of the fuzzy rules. The inference engine here exploited for the variables’ state update is the zero-order Sugeno method (Sugeno, 1985), in which the outputs of the rules are constant crisp values. To evaluate the next state of a variable over time, this method performs an aggregation of the output values produced by the rules, weighed according to the membership degrees of the antecedents of each rule. Additional information about the Sugeno inference method can be found in Supplementary Section S2. After the application of the Sugeno inference to all inner and output variables, their states are updated and the dynamic simulation of the DFM proceeds to the following time step. The process iterates until the simulation time reaches the maximum time step tmax set by the user. Given a set of linguistic variables, a perturbation π can be applied to the DFM as follows. For each variable belonging to π, a list of time intervals is specified. During each time interval, the Sugeno inference is disabled for that variable and the state update is performed by using a time-dependent arbitrary function.

2.2 FUMOSO: a general-purpose simulator of DFMs

FUMOSO is a novel open source and cross-platform software specifically designed and implemented for the definition, simulation and analysis of DFMs (see implementation details in Supplementary Section S1). FUMOSO is provided with an intuitive Graphical User Interface, realized to guide the users through the steps required for the creation of a DFM, that is, the definition of linguistic variables, linguistic terms, output crisp values, membership functions and fuzzy rules (see Supplementary Fig. S1). FUMOSO allows the user to enter all additional information required to simulate the dynamics of a DFM: (i) the simulation interval ; (ii) the functions that drive the dynamics of input variables; (iii) the initial state of all variables; (iv) the perturbation functions (if any). Optionally, for each variable involved in a perturbation π, the user can specify a list of time intervals in which that perturbation becomes active over that variable.

Given a set O of linguistic variables chosen as observed targets for the perturbation π, the effect of π on each variable is calculated as the difference between the state of o at the end and at the beginning of any perturbation interval. The overall effect of π on the set of all observed targets in O, during all time intervals that characterize π, can be calculated as a user-defined function .

Once the simulation is completed, FUMOSO plots the dynamics of any chosen system component, along with the shape of the membership functions and the weights of the firing rules involved in the update of that component at any arbitrary time step, according to the Sugeno inference method. The analysis of the simulation outcome is facilitated thanks to the possibility of creating groups of components, that is, subsets of variables whose dynamics are shown in the same plot.

2.3 Global optimization of DFMs

Given a DFM, FUMOSO allows to automatically explore the emergent behavior of the system in different scenarios, where the state of one or more linguistic variables is varied to simulate the effect of perturbations (e.g. drugs) in obtaining a desired system behavior. To this aim, a perturbation π can be defined by the user within a chosen time interval , by setting the state of (a subset of) inner variables to a specific value, which either belongs to the term set of each perturbed variable or is equal to the ‘unperturbed’ value. The latter corresponds to the state of each perturbed variable at the current time step, as evaluated by the Sugeno method. Given a set O of linguistic variables chosen as observed targets, a user-defined function can be defined to quantitatively assess the effectiveness of π in obtaining the desired system behavior. A sampling time instant , such that , can be set by the user to evaluate in between and tb. is used in FUMOSO as the fitness function for the global optimization of the DFM.

Global optimization allows to realize an effective and efficient exploration of the huge search space of possible perturbations of a DFM, whose dimension grows exponentially with the number of perturbed variables. To automatically investigate the search space and find out the optimal perturbation, FUMOSO is coupled with Simulated Annealing (SA) (Kirkpatrick ). The SA procedure integrated in FUMOSO starts from an initial perturbation π0—where all the variables are set to the ‘unperturbed’ state—and explores the neighborhood of potential perturbations. During the ith iteration of SA, a new putative perturbation is calculated, where η is a neighborhood function that randomly modifies the current perturbation by changing at most p variables. To be more precise, the perturbation π is modified by randomly changing the state of a selected variable taking a new linguistic term from the set of perturbable values of that variable. In the case of maximization problems, if then is accepted as the new perturbation and the process iterates; otherwise, is accepted if , where rnd is a random number sampled from the uniform distribution in , and T is a parameter that starts from an initial value T0 and linearly decreases to 0 during the iterations, progressively reducing the exploration capabilities of SA. The output produced by SA is the perturbation π characterized by the best fitness value . In this work, we considered and p = 2.

2.4 Experimental protocol

Mouse K-ras-transformed NIH3T3-derived cell line 226.4.1 and human breast cancer MDA-MB-231 (obtained from ATCC, Manassas, VA, USA) were cultured in DMEM containing 4 mM l-glutamine, 100 U/ml penicillin and 100 mg/ml streptomycin and pyruvate free (complete medium), supplemented with 10% newborn calf serum (mouse cells) or 5% fetal bovine serum (human cells). Cells were grown and maintained according to standard cell culture protocols. All reagents for media were purchased from Thermo Fisher Scientific. For analysis, cells were plated at a density of 3000 cells/cm2 in complete medium. After 18 h cells were washed with phosphate buffer saline (PBS) and incubated in growth medium (time 0) supplemented with 25 or 5 mM glucose (high glucose, HG) or 1 mM glucose (low glucose, LG). Cells were then treated and collected for analysis as described in the figure legends. To measure cell proliferation, harvested cells were counted using the Burker chamber. Where indicated, cell viable count was performed using Trypan Blue Stain 0.4%. Supplementary Section S3 contains detailed information about the used compounds, flow cytometric analyses, confocal microscopy, adhesion assays, Western blot analysis, transcriptome analysis and statistical analysis of the experimental data.

3 Results

3.1 Fuzzy logic model

The DFM defined in this work describes the cellular components that govern death and survival in K-ras cancer cells under progressive glucose depletion. This growth condition, influenced by highly interconnected processes with multiple levels of regulations (i.e. protein–protein interaction and modification, positive and negative feedback, etc.), and able to promote opposite effects on cancer cells (i.e. sensitivity versus resistance to chemotherapy), here is regarded as a complex and heterogeneous system to which apply our DFM approach, to the aim of efficiently identifying novel therapeutic treatments that maximize apoptosis over survival in cancer cells in such growth condition. The model can be depicted as a graph consisting of 25 nodes corresponding to heterogeneous components—e.g. proteins, small molecules and metabolites, biochemical pathways, cellular processes, output phenotypes—while edges between nodes indicate the known positive or negative regulations among these components (Fig. 1). To describe cell death processes as a result of glucose starvation in cancer cells, we considered:

Fig. 1.

Interaction network of the model of cell death and survival. Yellow circles represent metabolites and ions, green rectangles represent proteins, red rectangles represent pathways or cellular processes, light blue hexagons represent the system phenotypes related to cell death. Positive and negative regulations are pictured as arrows and blunt-ended arrows, respectively. Glucose, Ras-GTP and PKA are the input variables; survival, autophagy, apoptosis and necrosis are the output variables, while the remaining are inner variables. (Color version of this figure is available at Bioinformatics online.)

the main cellular processes and components involved in energy production, such as glucose, glycolysis, ATP, mitochondria and autophagy;

different mitochondrial processes and components, such as mitochondrial potential variation, ROS generation, mitochondrial complex I (CI) activity and B-cell lymphoma 2 (Bcl2) expression and activity;

some processes and components related to the ER, such as UPR, Ca2+, C/EBP-homologous protein (CHOP) and c-Jun N-terminal kinase (JNK);

processes and proteins involved in cellular adhesion, such as hexosamine biosynthesis pathway (HBP), N-glycosylation, attachment and Src;

proteins and processes involved in the regulation of cell death and survival mechanisms, such as Ras, extracellular-signal-regulated kinase (ERK), death-associated protein kinase 1 (DAPK), Bcl2, Beclin1 (BCN1), caspase 3 (Casp3), protein kinase A (PKA) and the phenotypes related to apoptosis, necrosis and survival of cells;

the protein Ras-GTP to mimic the hyperactivation of K-ras inside tumor cells displaying both the Warburg effect and PKA as a key node involved in cancer cells survival to glucose starvation.

Since not all of these components (e.g. apoptosis, survival, UPR, HBP, etc.) can be formally represented by a quantitative variable, and some interactions cannot be specified by means of kinetic reactions, we relied on fuzzy logic modeling to handle this heterogeneity and the lack of parameters by defining suitable linguistic variables that represent general concepts like concentration, activation or presence of a component. For each model component, we defined a linguistic variable and an associated set of linguistic terms to describe all the possible states of the variable (see, e.g. Supplementary Fig. S2, left side, for the ROS component shown in Fig. 1). To formalize the interactions existing among the components considered in the model, we exploited literature and expert knowledge. The complete list of linguistic variables and the corresponding output crisp values of their associated linguistic terms are available in Supplementary Table S1 and Table S2, respectively. A set of fuzzy logic rules was defined for each linguistic variable, for a total of 252 rules (see Supplementary Material). As an example, Supplementary Figure S2, right side, shows the fuzzy rules that describe the regulations acting on the ROS component given on the left side, according to the linguistic terms (e.g. Low, Medium and High) defined for the three components CI, DeltaPsi and ROS.

Three variables in the DFM of cell death and survival in glucose depletion represent the input of the system, namely, glucose, Ras-GTP and PKA: (i) Glucose is formally regulated by a custom update function that simulates glucose consumption; (ii) Ras-GTP is constantly kept to the High state to mimic the oncogenic hyperactivation of K-ras in cancer; (iii) PKA is either kept to the High or Low states to mimic the ability of cancer cells either to survive or die under glucose starvation (see Supplementary Section S4). Apoptosis, necrosis and survival represent the observable output variables of the DFM.

3.2 Model validation

The DFM of death and survival of K-ras cancer cell grown on progressive glucose depletion, defined on the basis of empirical data and the manual curation of literature, was experimentally validated against data obtained from cell cultures grown in different glucose availability, and in presence or absence of different protein and process modulating molecules (see Supplementary Section S5). Experimental data were obtained on NIH3T3 K-ras cells and a human model of breast cancer, the MDA-MB-231 cell line, which carries the oncogenic K-ras gene. The two cell lines were used in an exchangeable manner, since previously published results indicated that they have a comparable behavior as regard to numerous parameters associated with growth in 25 or 5 mM glucose (HG growth), and in 1 mM glucose (LG growth) (Gaglio ; Palorini et al., 2013a, b, 2016).

3.3 Perturbation analysis

A perturbation analysis of the DFM of programmed cell death and survival was performed in order to identify potential stimuli leading to both a reduced survival and an increased death by apoptosis in cancer cells. The total number of possible perturbations for the DFM defined in this work is , where 6 and 9 are the numbers of perturbable variables that can assume 3 and 4 possible states, respectively (see Supplementary Table S3). Since an exhaustive search of this huge space of perturbations would not be feasible, we exploited SA to automatically infer the minimal set of perturbations that allows to maximize apoptotic death in K-ras cancer cells. The optimization analysis was carried out both in the case of low and hyperactivated states of PKA, during different glucose availability. Among the solutions identified by SA, a set of promising single/double perturbations was selected and tested in the laboratory to assess their effectiveness in inducing death by apoptosis in MDA-MB-231 cells. A list of perturbations found by SA—validated either in this work or confirmed by previous experimental evidences—is given in Supplementary Tables S4 and S5, ranked by their fitness values. Details about the fitness function exploited to evaluate the effectiveness of the perturbations are given in Supplementary Section S6.

Figure 2 and Supplementary Figure S10 show the comparison between the simulated dynamics and the experimental data for single perturbations predicted by SA. In particular, we assessed the effects of UPR activation (Fig. 2) achieved by using 10 nM thapsigargin (thap), or CI inhibition (Supplementary Fig. S10) achieved by using 10 nM rotenone (rot) or 20 nM piericidin (pier) in high and low glucose availability, respectively (the experimental scheme for both perturbations is shown in Supplementary Fig. S9c). Both perturbations were chosen taking into account their rank and the existence of previous data indicating that the pharmacological over-activation of UPR pathway, as well as CI inhibition, can lead to cancer cell death (Cubillos-Ruiz ; Galluzzi ). For both perturbations, high or low activation of PKA was also experimentally analyzed. In particular, to avoid the endogenous activation of PKA under our experimental conditions, low PKA was achieved by cell treatment with the known PKA inhibitor H89 (Chijiwa ) (experimental details are available in Supplementary Fig. S9).

Fig. 2.

Assessment of the effects of perturbations (UPR activation) predicted by the global optimization algorithm. (a, b) Simulation outcome of the three main model output (apoptosis, necrosis and survival) upon UPR activation, either in (a) PKA Low state or (b) PKA High state. The perturbation was applied from time tb = 0 to the end of the simulation, and evaluated after a.u. (shaded area, see also Supplementary Section S6). (c) MDA-MB-231 cells, grown in HG, were daily treated with 10 μM FSK mimicking the PKA High state, or 5 μM H89 mimicking the PKA Low state and, upon 24 h, also with 10 nM thap (single treatment). Samples were evaluated for cell death at 48 and 72 h post-treatment by using trypan blue exclusion method. The experimental scheme is shown in Supplementary Figure S9c. All data represent the average of at least three independent experiments (±SD); *P < 0.05 (Student’s t-test). (Color version of this figure is available at Bioinformatics online.)

It is interesting to notice that the model properly predicted the experimental data. Indeed, an enhanced cell death was observed in the simulations (orange line for apoptosis and magenta line for necrosis in Fig. 2a and b) as well as in the experimental data upon both treatments (Fig. 2c). This increase is evident when moving from a high availability of glucose (left side of the simulation plots) to a situation of glucose starvation (right side of the simulation plots), and especially in the PKA Low state with respect to the PKA High state, suggesting an important role of PKA in cancer cell survival in acute UPR activation or in glucose starvation. These computational results are consistent with the experimental measurements, as represented by the graph bars in Figure 2c, which show a higher level of cell death in glucose starvation (72 h), with the fold change being higher when PKA is not activated. Similar considerations can be done for the condition of chronic CI inhibition, as shown in Supplementary Figure S10.

Analogously, Figure 3 shows the comparison between the simulated dynamics and the experimental data for the double perturbations predicted by SA. In particular, it shows the effects of UPR activation coupled with autophagy inhibition (Fig. 3a–c), N-glycosylation and HBP inhibition (Fig. 3d–f), N-glycosylation and autophagy inhibition (Fig. 3g–i). Both the simulation outcomes and the experimental data show an increase in cell death, especially when moving towards a state of glucose starvation, and the protective role of PKA. This result is evident by comparing, for example, the final states reached in each condition by apoptosis (orange line) and necrosis (magenta line) in Figure 3a, d and g, with the states reached by the same variables in Figure 3b, e and h. These computational results are consistent with the experimental measurements shown by the graph bars in Figure 3c, f and i.

Fig. 3

Assessment of the effects of perturbations (UPR activation and autophagy inhibition) predicted by global optimization. (a, b) Simulation outcome of the three main model output (apoptosis, necrosis and survival) upon UPR activation and autophagy inhibition, either in (a) PKA Low state or (b) PKA High state. The perturbation was applied from time tb = 0 to the end of the simulation, and evaluated after a.u. (shaded area, see also Supplementary Section S6). (c) MDA-MB-231 cells, grown in HG, were daily treated with 10 μM FSK mimicking the PKA High state and, upon 24 h, also with 10 nM thap and 20 μM chloroquine (CQ) (single treatment of both). Samples were evaluated for cell death at 48 and 72 h post-treatment by using trypan blue exclusion method. The experimental scheme is shown in Supplementary Figure 9d. (d, e) Simulation outcome of the three main model output (apoptosis, necrosis and survival) upon HBP and N-glycosylation inhibition, either in (d) PKA Low state or (e) PKA High state. (f) MDA-MB-231 cells, grown in HG, were daily treated with 10 μM FSK mimicking the PKA High state and, upon 24 h, also with 1 μM aza and 50 ng/ml tuni (single treatment of both). Samples were evaluated for cell death at 48 and 72 h post-treatment by using trypan blue exclusion method. The experimental scheme is shown in Supplementary Figure 9d. (g, h) Simulation outcome of the three main model output (apoptosis, necrosis and survival) upon N-glycosylation and autophagy inhibition, either in (g) PKA Low state or (h) PKA High state. (i) MDA-MB-231 cells, grown in HG, were daily treated with 10 μM FSK mimicking the PKA High state and, upon 24 h, also with 50 nM tuni and 10 μM CQ (single treatment of both). Samples were evaluated for cell death at 48 and 72 h post-treatment by using trypan blue exclusion method. The experimental scheme is shown in Supplementary Figure 9d. All data represent the average of at least three independent experiments (±SD); *P < 0.05, **P < 0.01 (Student’s t-test). (Color version of this figure is available at Bioinformatics online.)

To analyze the combinatory effect on cell survival we used, when possible, sub-toxic concentrations of the compounds, previously determined experimentally on this cell model (see, e.g. tuni and aza). In addition, in these sets of validation experiments, we evaluated also whether the model could predict the endogenous PKA behavior, since no H89 inhibitor was used (experimental scheme in Supplementary Fig. S9d). Combination effects of selected compound pairs exceeded the effects of single compounds for the majority of combinations, and well fitted with the model outcome also in the absence of exogenous PKA inhibition. Thus, these results corroborated the capability of the DFM, coupled with the optimization algorithm, in predicting the system’s response in perturbed conditions.

4 Conclusion

Complex biological systems are characterized by emergent, non-linear dynamic behaviors that arise from negative and positive feedback regulations among a huge number of different molecules and processes in cells. The elucidation of the mechanistic interactions that govern the correct functioning of cells or that, in contrast, can induce the onset of a disease, is an extremely challenging task from both the experimental and the computational perspectives. This is strikingly true when the biological system under investigation is poorly understood, or it cannot be easily subject to accurate experimental measurements, or the interplay between its intrinsic negative and positive control mechanisms is so complex that the system can be analyzed only by considering a subset of (usually, homogeneous) components, such as metabolic traits, gene expression, etc.

In this context, we presented a novel computational method able to address two open issues in the field of biochemical system modeling, that is, taking into account the heterogeneous nature of all biological systems, and dealing with the lack of quantitative parameters that generally prevents the precise characterization of most cellular components and processes. This modeling approach may be adopted to incorporate data generated from different ‘platforms’ using knowledge-based rules (e.g. transcriptional and proteomic data, microscopy imaging), and can integrate multiple data types together, i.e. signaling protein activation/inhibition, post-translational modifications, small molecule variation, end-process activation or inhibition, etc. As such, our computational method not only solves some relevant issues related to the modeling and dynamical simulation of heterogeneous systems, but it is also able to provide valuable predictions that facilitate our understanding in controlling cellular systems. Our method, based on the coupling between fuzzy logic modeling and a global optimization algorithm, allows for the definition of models in a readable and simple format, through the use of linguistic variables. Here, in particular, we showed that our approach was able to predict the behavior of K-ras-transformed cells grown either under progressive glucose depletion or in different perturbed conditions, as well as to identify possible novel cancer therapeutic treatments.

The model simulations showed that targeting the HBP, its downstream route controlling protein N-glycosylation, the ER processes or the mitochondrial CI, alone or in combination with PKA inhibition, cause a significant increase in cancer cell death. Noteworthy, all computational results were experimentally substantiated here, or matched literature results. For instance, it has been recently reported that HBP inhibition in HG growth induces cell death in breast cancer cells and xenograft mice (Ricciardiello ), as also shown in this work by computational and experimental data. Moreover, we demonstrated that inhibition of N-glycosylation has a profound effect on breast cancer cell survival, validating other recent observations (Chiaradonna ; Serrano-Negrón ; Wang ; You ). In both cases, the model correctly predicted the outcomes of new experiments that were not used in model construction, thus confirming the reliability of our computational method.

It is noteworthy that the relevance and the validity of our method is supported by its performance in predicting not only recognized cell death-inducing stimuli, but also unrecognized stimuli equally leading to cell death. For instance, a sub-toxic amount of tuni, coupled with both a sub-toxic amount of the HBP inhibitor aza or the autophagic inhibitor chloroquine, enhances cancer cell death. Remarkably, both treatments are strongly attenuated by exogenous PKA stimulation, implying the involvement of this pathway in ER stress response, at least in our experimental conditions. Confirming results were shown also about the protective role of PKA upon ER stress induction by thapsigargin that, in combination with PKA inhibition, induces cancer cell death, an effect that is strongly impaired by exogenous activation of PKA. Previous data indicated that PKA activation protects cancer cells from death induced by glucose starvation (Palorini , 2016). Here, we revealed a protective role of PKA in cancer cells under acute ER stress. While this protective mechanism has been shown to be active in mouse embryonic fibroblasts (Aguileta ) and hepatocytes (Li ), to the best of our knowledge it has never been described in cancer cells, further supporting the predictive value of our method. Interestingly, also simultaneous CI and PKA inhibition induces massive cancer cell death, an effect that is prevented by PKA activation. Therefore, the model also predicted that PKA activation is involved in mitochondrial CI function and, more in general, in OXPHOS activity, corroborating previous results (García-Bermúdez ; Lark ; Ould Amer and Hebert-Chatelain, 2018; Papa ). Altogether, our results provide the first evidence of a protective role for PKA against several treatments mimicking cellular stress conditions, such as ER stress, N-glycosylation inhibition, mitochondrial CI inhibition and glucose starvation in cancer cells. These findings will hopefully pave the way for the use of new and more specific PKA inhibitors in cancer therapy.

We further underline that the coupling of fuzzy logic modeling with optimization algorithms is a promising tool to uncover new potential therapeutic targets by assessing, in an automatic way, the response of the system to an extensive number of perturbations, therefore both reducing the costs and facilitating the design of laboratory experiments. The dynamic nature and predictive power of DFMs could prove useful in assessing the effects of different types of perturbations, such as drugs or mutations, on the behavior of the system under examination in many application fields (e.g. medicine, pharmacology, etc.). In the future, this novel computational method could be exploited to preliminary assess the effects of FDA-approved drugs, especially in combination with other metabolic drugs, on the survival of resistant cancer cells. The great flexibility of fuzzy logic could also be exploited to integrate different model formalisms to define complex hybrid models, able to both represent different layers of biological complexity (at the functional, temporal, or phenomenological level) and leverage precise kinetic information when available (Spolaor ). Namely, with this novel computational method we aim at opening the way to filling the gap between quantitative (mechanism-based) models and qualitative (logic-based) models, in order to simultaneously exploit the peculiar advantages provided by each modeling approach.

Funding

This work was supported by grants from the Italian Ministry of University and Research (MIUR) [SYSBIONET-Italian ROADMAP ESFRI Infrastructures to L.A. and F.C., PRIN2008 to F.C.] and A.I.R.C. IG2014 [Id.15364 to F.C.]; SYSBIONET-Italian ROADMAP ESFRI Infrastructures to G.V. and H.D.V. (fellowships); Italian Ministry of University and Research to R.P. and F.R. (fellowships); CAPES 9281-13-4 to H.D.V. (fellowship); and A.I.R.C. IG2014 [Id.15364 to G.V. (fellowship)]. Financial support from MIUR through grant ‘Dipartimenti di Eccellenza-2017’ to University of Milano-Bicocca, Department of Biotechnology and Biosciences is also acknowledged.

Conflict of Interest: none declared.

Click here for additional data file.