Received: 09 February 2018 Accepted: 04 September 2018 Published: 09 October 2018

#### Citation:

Thobe K, Kuznia C, Sers C and Siebert H (2018) Evaluating Uncertainty in Signaling Networks Using Logical Modeling. Front. Physiol. 9:1335. doi: 10.3389/fphys.2018.01335

**250**

system is gathered and validated on new data (**Figure 1B**). A main issue when building these models is that uncertainty cannot be included into the model, e.g., if an influence between two components is controversial in the literature. Since only one model is created, the modeler needs to make an assumption neglecting the uncertain information. A second popular strategy for modeling is to use a top-down approach where the model is inferred directly from data, but here prior knowledge about the system is neglected (De Smet and Marchal, 2010).

As a consequence, alternative approaches have become more popular, where uncertainty is included into the modeling process by either building more than one model or by adapting the model through training. For exploring a wide range of models that is able to show a certain dynamical behavior, called family or pool of models, there are different methods available. The group of Saez-Rodriguez et al. developed a software CellNOptR to train a candidate model to data, accounting for topological uncertainties (Terfve et al., 2012). The output is a family of models selected for an optimality criterion, but cannot guarantee completeness due to stochastic search. The software caspo by the group of Siegel et al. uses Answer Set Programming to infer a family of logical models from experimental data based on optimization, where a tolerance accounts for experimental noise. The resulting family of models then represents all optimal models that reproduce the data and the software provides several analysis tools to explore properties of the models, such as classification for input/output behavior or experimental design (Videla et al., 2017). A similar method using time series data for inferring a model pool showed to be more precise than caspo (Ostrowski et al., 2016). A different approach was developed by our group, where uncertainty in parameters of the model, such as an uncertain sign of an edge, is encoded into the model definition (Klarner et al., 2012) and all possible models that arise from this uncertainty are enumerated. Subsequently the models are tested for satisfiability for data without an optimality criterion (**Figure 1A**), which was implemented using efficient formal verification techniques in Tremppi (Streck et al., 2015) and in TomClass (Klarner, 2014). Even though we employ the software from our group for the analysis in this paper, one could apply different software along the pipeline for building the model pool or analyzing it.

While computing model pools and testing them for data sets is computationally challenging, the analysis of potentially hundreds or thousands of models is not straight-forward in terms of the biological interpretation. Thus we propose a hypothesisdriven approach for specific biological questions, where the use of model pools allows us to test multiple hypotheses at the same time and analyze their interdependencies. Mathematical models are artificial constructs used to help understanding biological processes. In order to receive meaningful results from a modeling study, the biology needs to be transferred into mathematics and the results need to be interpreted from a biological perspective. In this paper, we address this task of incorporating biological information into the formalism by expanding the workflow in **Figure 1** to a four-step pipeline. At first, the process of bottom-up model building formalizes the biological phenomena into a prior knowledge network, which we call system initialization. Here, the regulatory graph and the logical equations are derived from literature information. Then, the objective formalization includes the aim of the study into the model setup, e.g., by adding extra components or edges. After generating the model pool, the topdown filtering process uses biological data that is not restricted to be of a specific type such as steady-state or input-output behavior. However, it requires a data formalization step. Finally, the pool analysis examines the specific pool for new biological insight.

In previous work, we presented parts of this pipeline, i.e. the objective of investigating crosstalk between two signaling pathways in Thobe et al. (2014), as well as challenges for data discretization and analysis in Streck et al. (2015) in context of a specific software. Here, we generalize and expand this pipeline by two additional objectives and analysis methods. Especially in the context of signaling processes in cancer cells, the identification of driver mutations is of great interest (Bozic et al., 2010), thus one aim of our framework is to identify driver mutations by a change in the logical function. The second aim presented is testing the effect of drugs by introducing them as new inputs to the system. Analyzing pools containing possibly hundreds or thousands of models is challenging. Here, we show a classification analysis to structure the resulting models toward interesting features, as well as extracting minimal mechanisms for a more detailed view on the models.

We apply this pipeline to model two central signaling processes involved in cancer, the mitogen-activated protein kinase (MAPK) cascade and the mTOR pathway (Shaw and Cantley, 2006; Saini et al., 2013), in two renal cancer cell (RCC) lines. Both cell lines were treated with the Raf-inhibitor Sorafenib

yet displayed a differential response in terms of apoptosis induction (Kuznia, 2015). We hypothesized that the difference between these cell lines might be caused by distinct wiring of MAPK and mTOR signaling, which were shown to be connected via crosstalk (Mendoza et al., 2011; Aksamitiene et al., 2012). A rich dataset of time series measurement of key components in both pathways was generated using a high-throughput method, which was the foundation for the complex analysis presented in this study.

This paper is organized as follows. The Methods section first gives a brief introduction on the logical modeling framework and a detailed description on the pipeline we developed. In the Results section, the application on a signaling network is demonstrated, where first the model building process with the corresponding biological background is given, the data processing procedure is described and the results of the analysis are presented. Additionally, the biological interpretation is discussed and future experiments are suggested to wrap up the application section. Finally, the Discussion section exploits advantages and shortcomings of the method showing potential future extensions.

# 2. METHODS

# 2.1. Theoretical Background

The formalization of logical modeling for biological systems was introduced by Kauffman (1969) and further refined by Thomas (1991), which is the base for our work. However, we expanded this formalism to incorporate uncertain information leading to model pools (Klarner, 2014; Thobe et al., 2014, 2017; Streck, 2015; Streck et al., 2015).

#### 2.1.1. Logical Modeling

The topology of a biological system is defined as a directed graph R = (V, E, l), called interaction graph (IG), where the nodes V = {1, ..., n} represent the components of the system that are connected by edges e ∈ E ⊆ V × V called interactions, which represent a regulation of one component by another. The components adapt discrete values, called activity levels, and we consider Boolean networks (BN) with two levels assigned to each component B = {0, 1}, where 0 means inactive and 1 stands for active. By assigning activity levels to every component of the network, the state of the system s is defined by s:V → {0, 1}, ∀v ∈ V :s(v) ∈ B. Here, the notation of a state is specified as a sequence in the order of V. In our approach, we add information about the nature of a regulation to each interaction using edge labels l : E → {+, −, ¬+, ¬−} (adapted from Klarner, 2014). In application, the labels {+, −} are assigned to edges that represent well-known information, e.g., textbook knowledge, and are therefore required to be present in every model, which we call essential. In contrast, the labels {¬+, ¬−} are assigned to interactions that carry uncertainty, i.e., we not sure whether this interaction is present or not, which we call optional. However, we assume that the sign of an edge is known and exclude edges with unknown or ambivalent sign due to complexity.

Having defined the wiring of the network, the regulation of a component by its predecessors is defined by a logical function. The conditions describing when a component becomes active can be expressed using the logical operators ∨ (OR), ∧ (AND), and ¬ (NOT) in a formula fi for every component v ∈ V consistent with the edge labels. This means that variables j are literals in f for component i, if j → i is a possible edge. Then a positive edge label has to cause an increasing value in the target component at some point, whereas a negative edge label has to cause a decrease. For optional edges, the increase or decrease can occur or that value is constant. However, in case the regulation of a component is uncertain of a component has optional incoming edges more than one model can be build from the available information. Then the set or all logical equations that are consistent with the edge labels are created and form the so-called model pool. An example is given in **Figure 2**.

#### 2.1.2. Dynamical Behavior and Model Checking

In order to compare biological measurements with the dynamic behavior of the models, we need to define the transition from one state to another to generate the systems behavior over discrete time steps. For this aim, different update strategies have been developed, where some make assumptions on the timing of events, e.g., in synchronous update all components change in one transition, and others restrict the ordering of events, e.g., stochastic updates randomly update a component. Here, we employ asynchronous update, which is the least restrictive strategy at the cost of being computationally expensive (Thomas, 1991). In this strategy only one component can change its value per transition step, which means for fv(s) = ¬sv for a state s = (s1, . . . ,sv, . . . ,sn) denote with s v = (s1, . . . , ¬sv, . . . ,sn) the state which differs from s in the value of the component v. If no component changes f(s) = s a steady-state of the system is reached. This update schedule produces every possible trajectory emerging from a state, thus the dynamics are non-deterministic which can be visualized in the so-called state transition graph (STG). Here, the states are the nodes of the graph and the transitions are edges.

After building the model pool from the available information, we want to filter those model that are in agreement with observed experimental data. Depending on the utilized software, either the data can be implemented as continuous values (e.g., Terfve et al., 2012) or needs to be discretized. Here, we want to describe two different kinds of biological data: time-series measurements and steady state observations. For this aim, we use temporal logics that are able to describe an ordering or a sequence of events in time, where computation tree logic (CTL) can cope with nondeterministic sequences (Clarke et al., 1986) and are therefore suited to explore the STG. For time-series measurements we encode a series of states that should exist at some point in the future and for steady states we encode the state of a component(s) that should hold for every state in the future. These formulas are then tested on the STG of the models using model checking. This process can be computationally expensive, since the state space exponentially increases with the number of components, also the number of models can quickly add up to thousands of models. For this reason, an efficient model checking software should be employed, e.g., Tremppi (Streck et al., 2015) and TomClass (Klarner, 2014).

some path from an initial state the X holds true at some point. STGs in (F) of the three models in the pool show the process of model checking for the CTL formula indicated by green states and edges, where the second model is not in agreement with the data.

#### **<IP_ADDRESS>. Toy example**

In **Figure 2**, model definition to model checking is visualized for a toy example. Here, an IG with four components is given, where the regulation of components A, B, and D is known, indicated by the edge labels and the corresponding functions. Component C has an uncertain regulation by component D, therefore the edge is labeled as not inhibiting and the function for C is undefined (**Figure 2A**). The resulting model pool then contains three different models that arise from the edge label. The process of temporal encoding of data is shown for time-series measurements, which is discretized by a threshold and encoded as CTL formulas, where the CTL operator exists finally is used. This operator states that the measurements must lie on one path in the ordering of the measurements in the STG, where there is no restriction made on how many states are visited in between the measurements. The CTL language offers more operators that could be employed depending on the type of data (Klarner et al., 2012), e.g., reflecting that the measurement frequency was so

high that one can assume that all qualitative activity changes of each component have been captured. However, in this paper we used the most conservative form. Finally, the process of model checking is visualized in **Figure 2F** showing that one model is not in agreement with the data.

# 2.2. Pipeline for Modeling Uncertain Systems

Based on the very general workflow in **Figure 1**, we want to formalize in more depth how specific biological questions can be addressed using model pools (**Figure 3**). In comparison to the traditional workflow for building and analyzing one model, there are similarities and differences. The main difference is that the standard bottom-up approach creates one model to test a single hypothesis, which is validated using (formalized) data and subsequently analyzed. Using model pools, the system initialization and data formalization remains the same, but we can test multiple hypothesis at the same time. This leads to a higher complexity in both the formulation of the aim of the study and the analysis of the model pool for biological information. To this end, the workflow has been adopted and specified to address common analysis themes for signaling networks. The resulting pipeline, shown in **Figure 3**, contains four steps: system initialization, objective formalization, data formalization, and pool analysis.

In the following, we provide a detailed formal description for the objective formalization and pool analysis. We assume that for the system initialization, available information is gathered and classified according to the theory presented in Section 2.1 for components, edges, and edge labels. Moreover, insight on regulations of components can be included by defining logical functions. This first step results in the prior knowledge network (PKN), which forms the starting point for our analysis (Saez-Rodriguez et al., 2011).

#### 2.2.1. Objective Formalization

In the second step of our pipeline, we want to include the objective of the analysis into the PKN. In a simple setup, this could mean adding optional edges as hypotheses, but there are also more complex aims that require changes in the PKN. The first objective we identified, was to examine crosstalk between two pathways while preserving the dynamical properties of each pathway, presented in Thobe et al. (2014). Here, we want to present two different objectives: Finding driver mutations and drug testing.

#### **<IP_ADDRESS>. Finding driver mutations**

Cancer cells often accumulate mutations that are distinguished as either driver or passenger mutations. A lot of effort has been made to identify the driver mutations, since they are assumed to be a major cause for cancerous behavior (Greenman et al., 2007; Bozic et al., 2010). This abnormal behavior is due to the fact that the mutations affect the protein they are encoding in quality, changed sequence of the protein, or quantity, such as overexpression or knock-out of a gene. These effects cause changes in the regulatory network leading to an insensitivity of the component from its regulators, e.g., constantly active receptors. We aim to identify these changes in the regulation of a component in our approach, which we were able to confirm in previous work (Streck et al., 2016).

We account for mutations in components of a model with uncertain effect by setting the respective incoming and outgoing edges to optional even if these connections are textbook knowledge. In case more detailed information on the effect of the mutation is available, only a subset of edges can be set to optional. Formally, the network R = (V, E, l) is defined as in Section 2.1, where the set of mutated components is given by V m ⊆ V with edges E m = {(u, v) ∈ E | v ∈ V m ∨ u ∈ V m}. The labeling l m of edges in E m is set to:

$$\mathcal{E}^m(\mu, \nu) = \begin{cases} \neg + & \text{if } (\mu, \nu) \in E^m \text{ and } l(\mu, \nu) = -\\ \neg - & \text{if } (\mu, \nu) \in E^m \text{ and } l(\mu, \nu) = + \end{cases}$$

Thus, the affected edges are allowed to either stay the same or lose their function in the resulting model pool. If in the specific pool an incoming edge is not observable in any model, the mutated component becomes independent from its inputs and the function of the component can either be set to 0 or 1

l

indicating a loss-of-function or constitutively active mutation, respectively. A lost outgoing edge of the mutated component can indicate that the mutation affected the protein structure which can result in a dysfunctional protein. However, we do not account for gain-of-function mutations in this set-up, since this would require to add new edges to the model or change the sign of an edge. This would strongly increase the complexity of the study and should be addressed only based on suggestive data in a case by case way, which is no fit for this general set-up.

#### **<IP_ADDRESS>. Drug testing**

This objective aims to test qualitative effects of drugs on pools of models, without knowledge of the "true" network. Especially in cancer research, combinatorial therapies have become increasingly popular to enhance efficiency and overcome resistances (Ho et al., 2012; Manchado et al., 2016). Since we cannot represent concentrations or generally quantitative effects, the questions we want to address can be formulated as: where do we have to interrupt the signaling process to achieve a certain outcome. A similar study was done by Klinger et al. where they predicted the model structure and treatment quantitatively, however, the predictions resulting from the study were qualitative nature (Klinger et al., 2013).

For this approach, we introduce drugs as new components to the PKN and connect them with an inhibitory edge to their target, since they are supposed to suppress the activity of their target. For the network R′ = (V ′ , E ′ , l ′ ) with f ′ as given logical equations, an extended set of components V is given by V ′ ∪ V D where v D ∈ V D is a set of drug components. The interactions of the network are given by E = E ′ ∪ E D where new edges E D are added, which contain an edge for self-activation for each new component to create the drug as input and an inhibitory edge from v D to its target u, since the drug suppresses the activity of its target. Similarly, the set of labels is composed of the labels of the original network and the additional labels for the drug components, where known interactions are labeled with an essential label and uncertain effects with an optional label.

We can also include available information about the drug's mode of action into the logical function of its target. Usually drugs are selected to have a dominant influence on their target, for example through binding or modification it fails to interact with its former regulators. In case the logical equation of a drug target is known, we can directly translate this dominant effect on the target u in a new logical equation:

$$f\_u = f'\_u \wedge \neg \nu\_D.$$

However, if detailed information about the biochemical properties of the drug on the target and other regulators is missing, the logical equation of the target is not defined and all possible regulations are generated in the pool.

#### 2.2.2. Pool Analysis

After building the generic model pool from the PKN, this pool of models gets reduced for those models that are valid for data. Depending on the software, the data needs to be processed to apply it to logical models usually by discretization (Dimitrova et al., 2010; Gallo et al., 2015). As a result, we receive one or more specific model pools that need to be analyzed. For this aim, different kinds of analysis tools can be employed depending on the aim and the size of the resulting model pool such as statistical analysis (Thobe et al., 2014; Streck et al., 2016) or optimization (Terfve et al., 2012; Videla et al., 2017). Here, we want to present an analysis approach that allows a closer look at classes of models as well as single models.

#### **<IP_ADDRESS>. Classification**

Depending on the study, this pipeline can lead to specific model pools that contain too many models to analyze them by hand. This analysis step aims to get an intuition for commonalities or differences of models within the model pool, with respect to properties of interest. For this goal, properties such as validity for data or presence of an optional edge can be annotated to each model by e.g., using a database. Then, we can group sets of models into classes and compare them according to these properties to find difference between sets of models, for example we could observe that two optional edges are present in the model pool but occur mutually exclusive.

Here, the model pool is stored in a database and SQL queries are used to classify the models. For the queries, properties or a list of properties can be used as a classifier and are defined in the parameter Classes. Also we can restrict the pool to a subpool using the parameter Restriction, where we can select models for their property, e.g., only including all models that carry an optional edge. Mathematically, the analysis finds subsets of models that have a non-empty intersection and computes the cardinalities of these sets (Klarner, 2014). For this aim, an SQL query is generated using statements of the form:

SELECT DISTINCT Classes FROM models WHERE Restriction,

where SELECT DISTINCT computes all combinations of labels, i.e., subsets, of the selected Classes in the database models, possibly restricted using WHERE. Additionally, COUNT is used to determine the cardinality of each subset, i.e., the number of models in a class later denoted as size of a class. It is possible that classes are empty if there exists no model in the pool with a particular label combination.

#### **<IP_ADDRESS>. Minimal models**

While the classification gives a broad overview on the model pool, we also want to look at single models in the pool. The selection of models can be motivated by the objective and the biological background, by the classification analysis or by general criteria such as minimality. The criterion minimal can be interpreted in different ways: structural minimality in terms of number of edges, functional minimality in terms of shortest logical equations, or models that require least number of transitions to fulfill data. Each minimality can be interesting to regard separately or in combination. Structural minimality is a common biological assumption, where the system is assumed to have evolved in an energetically optimal way and is therefore sparse. Along with the number of interaction partners, the complexity of the regulation formulated as logical function can be assumed to be rather simple. Previous studies often used fixed rules for creating these

Technically, this analysis counteracts the problem of overfitting. In general, the more degrees of freedom are available to a system, the easier is can produce various dynamics, thus our method has a bias toward building dense models. It is therefore beneficial to identify minimal structures or functions.

#### **<IP_ADDRESS>. Interpretation of analysis results**

Finally, the results from the analysis need to be transferred and interpreted to gain biological insight, which is not straight forward. Since the models are qualitative, the level of abstraction is high and the fact that we are looking at pools of models increases the complexity. However, by specifying clear objectives and predefined analysis options, the pipeline guides the modeling process and can deliver valuable information for experimental design or further modeling steps (Streck et al., 2015; Thobe et al., 2017).

# 3. RESULTS

# 3.1. Application on Growth Factor Signaling in Renal Cancer Cells

After presenting a pipeline to build model pools for different objectives and analysis options, we wanted to apply this pipeline to model growth factor signaling in two renal cancer cell lines, MZ1851RC and MZ1257RC. Motivation for this study was an observation that cell line MZ1851RC showed apoptosis after being treated with the drug Sorafenib while MZ1257RC seemed to be resistant (Kuznia, 2015). Sorafenib was developed to inhibit pathways controlling proliferation and cell survival and was shown to have anti tumor activity in colon, breast, and nonsmall lung cancer (Wilhelm et al., 2004; Gadaleta-Caldarola et al., 2015). The multikinase inhibitor Sorafenib was designed to suppress activity of Raf kinases in the MAPK pathway (Liu et al., 2006), however, it also affects a wide variety of receptor tyrosine kinases (RTKs) (Wilhelm et al., 2004). Very recently, Sorafenib was shown to inhibit the IGFR in vitro (Yaktapour et al., 2013), which indicates that Raf comprises an uncertain drug target in the renal cancer cell lines tested with our approach.

A second uncertainty was introduced by a mutation in the component mTOR in cell line MZ1851RC, but the effect of this mutation is unknown (Kuznia, 2015). A third uncertainty was caused by crosstalk between the MAPK pathway and PI3K signaling (**Figure 4A**), which was shown to compensate drugging of one of the pathways (Mendoza et al., 2011; Aksamitiene et al., 2012). Thus, the overall aim of this study was to clarify if the deviating behaviors are caused by differences the cellular wiring, which effect the mutation has and which targets Sorafenib is affecting.

# 3.2. Objective Formalization

The objective of this study splits into three different aims: investigating crosstalk between MAPK and PI3K pathways, finding the target for Sorafenib, and clarifying the effect of the mutation in mTOR. The PKN was extracted from literature, where the MAPK model was based on work by Kholodenko (2000) and the PI3K model was adapted from Courtney et al. (2010), also it is an adaption from a previous study (Thobe et al., 2014). For investigating the wiring between MAPK and PI3K pathway, candidate crosstalks were added. In detail, strongly activated MAPK signaling was found to cross-activate PI3K signaling, i.e., Erk was observed to phosphorylate Tsc2 suppressing it and Erk was also shown to phosphorylate Raptor, where both crosstalks activate mTORC1 signaling similarly to Akt (Roux et al., 2004; Winter et al., 2011). For simplicity, we summarized this effect to one crosstalk. Moreover, a crossactivation of EGFR on PI3K through Ras was shown, which is downstream of EGFR and upstream of Raf (Wong et al., 2010). A study of Will et al. found that PI3K inhibition, but not Akt inhibition, causes rapid decrease in wild type Ras activity and in Raf/Mek/Erk signaling concluding that PI3K cross-activates the MAPK cascade (Will et al., 2014). For the PKN, the crosstalk edges were labeled as optional edges and the edges within a pathway were assumed to be essential, shown in **Figure 4B**.

In order to test the effect of Sorafenib, it was added as additional input to the system as well as optional edges to possible targets: Raf, EGFR, and IGFR. Note that EGFR as Sorafenib target is a hypothesis and not based on experimental data. Moreover, one cell line, MZ1851RC carries a mutation in mTOR with unknown effect for mTORC1, thus the outgoing edge to IGFR was set to optional. A full list of optional edges in given in **Figure 4D**, also for some components the logical function can be set, since they only have one regulator (**Figure 4C**). All other components have undefined logical functions, which gives rise to the generic model pool.

Moreover, components that were neither measured nor perturbed were excluded from the model to reduce the complexity. For example, Mek and Tsc were not considered in the model, since both were lined up in a cascade as components with single input and output, thus deleting them does not pose problems for the model dynamics.

# 3.3. Data Formalization

#### 3.3.1. Experiments Show Differential Behavior of Cell Lines

For our investigation, we used two different data sets: Western blot measurements of mTORC1 activity over time and a high throughput assay both published in Kuznia (2015). In the western blot measurements, the activity of mTORC1 was measured by its targets p70S6K (S6K) and S6RP in MZ1257RC and MZ1851RC cells. Here, the cells were either treated with DMSO or Sorafenib and the phosphorylation of the mTORC1 targets was measured over time. Regarding the measurements until 12 h, MZ1257RC cells showed a significant decrease in phosphorylation levels for S6K and S6RP. However, MZ1851RC cells only showed a reduction in S6RP phosphorylation for later time points, but the phosphorylation of S6K remained high. The 24 h time point is not considered, since we are only interested in signaling effects and this measurement is likely to be influenced by transcriptional effects. S6K was used as the read-out for the mTORC1 activity in the formal encoding of the Western blot data as CTL formulas WB.DMSO, WB1257Sora and WB1851Sora in the **Table 1**. Here, both cell lines show active mTORC1 for

DMSO treatment throughout the measurements, thus a steady state was assumed and encoded in the CTL formula accordingly. For Sorafenib treatment, cell line MZ1257RC shows a steady state with decreased S6K phosphorylation, therefore mTORC1 was set to 0. In contrast, cell line MZ1851RC has stable S6K phosphorylation, thus mTORC1 was set to 1.

After observing differences in the activity of mTORC1 in the Western blots toward Sorafenib treatment, we wanted to investigate where the differences in the upstream regulation of mTORC1 originate from. For this aim, a high throughput approach using the Bio-Plex R system was applied (Kuznia, 2015). Here, the cells were unstimulated and not starved but treated with Sorafenib or DMSO and measured at different time points over a total period of 36 h in two experiments. In detail, the activity of the PI3K/mTORC1 signaling pathway was measured by the phosphorylation of Akt, and p70S6K as well as the MAPK activity was determined through the phosphorylation of Erk. Moreover, the receptors EGFR, and IGFR were included into the experiment, since we were interested whether the receptors are targeted by Sorafenib and to account for the feedback processes. For the complete dataset see Kuznia (2015), processed data is listed in the **Supplementary Table 1**.

#### 3.3.2. Discretization of Time Series Data

In order to fit the models in to pool to the time series measurements, the data needs to be discretized. The choice of discretization method is influenced by the kind of data acquired and the experimental method used, for example with large data

#### TABLE 1 | Filtering model pool using model checking.

CTL formulas derived from Western blot and Bio-Plex® experiments. For denoting the CTL formulas, the following semantics are used: EF(X): is a CTL operator exists finally. This states that on some path from an initial state the X holds true at some point. AG(X): is a CTL operator all globally. This states that X has to hold for all future states, i.e., X is in a steady state. v=b: where v ∈ V, b ∈ B states that value of a component v is set to b. IS: declares the initial state and is a list of boolean constraints on the values of the components. A state is considered initial, if all the constraints are satisfied.

sets statistical methods provide good results, but with small data sets the choice is more difficult (Dimitrova et al., 2010). Here, we opted to show a simple approach by using the arithmetic mean as threshold. More specifically, the data was discretized by defining for each experiment e and component v, a threshold

$$\theta\_{\rm ev} = \frac{\sum\_{l} \sum\_{\chi} m\_{\rm evtx}}{|t||\chi|}$$

where, mevtx is the measured activity of component v in the experiment e with treatment x. The total number of time points is |t| and the total number of treatments is |x|.

Since the cells were cultivated and treated in parallel, the phosphorylated levels for both treatments were expected to be comparable. Thus, the threshold for e.g., Erk is the same mean value under both Sorafenib and DMSO treatment within each cell line for each experiment. Moreover, the standard deviation for each component was calculated in order to avoid the problem of discretizing a component that does not change over time. By looking at small standard deviations relative to the mean, IGFR measurements for MZ1257RC in both experiments were identified as problematic (see **Supplementary Table 1**). Comparing the IGFR levels between the cell lines, we decided to exclude this data.

Since we are interested in the signaling processes, only measurements until 8 h were included. The resulting CTL formulas are listed in **Table 1**, where all Bioplex measurements were encoded as transient states, due to the fact that they changed throughout the 8 h of measurement. An exception is the data set Bp1257DMSO, which was encoded as steady state (see **Supplementary Table 1** MZ1257RC-DMSO Exp1). Note that the discretization of data is not always straight-forward, thus we excluded data which was problematic mathematically (such as IGFR) or had poor quality in the measurements.

#### **<IP_ADDRESS>. Robustness of results**

As a basic test of robustness with respect to the discretization method being used, we additionally performed a discretization by median instead of mean value. This change in the discretization TABLE 2 | Number of models consistent with CTL formulas.

Filtering for CTL formulas gives pool sizes as the number of models in agreement. Rp.1851 and Rp.1257 are the cell line specific pools as the intersection of the data sets shown in (A,B), respectively.

threshold had a negligible effect on both cell lines: 5.7% of the boolean values changed for MZ1257RC and 7.1% for MZ1851RC. Furthermore, we repeated the subsequent analysis using the median discretization, and observed only minor changes in the size of the model pool for cell line MZ1257RC and no change in the resulting biological interpretation of that pool.

# 3.4. Pool Analysis

After deriving the PKN from the literature and including the objectives of the study, the generic model pool was created. As a result from combining of all optional edges and logical equations the pool contains 19,404 models. In order to find biologically relevant models, the third step of the pipeline generated the specific pool(s) by filtering the generic pool for those models that are able to simulate experimentally observed behavior for the two RCC cell lines.

#### 3.4.1. Cell Line Specific Model Pools

Each CTL formula has a non-zero pool size and is therefore feasible for our analysis (see **Table 2**). To determine the cell line specific models, we calculated the intersection of the different subpools for cell line MZ1257RC as Rp.1257 and for cell line MZ1851RC as Rp.1851:

Note that both pools are required to fulfill WB.DMSO, since this dataset was identical for both cell lines. Although the single CTL formulas resulted in relatively large pools, containing 5,000– 15,000 models, the intersection for the cell line specific pools shows a strong reduction with 1017 models for Rp.1257 and 293 models for Rp.1851 (see **Table 2**). Thus, there exists a cell line specific pool for each cell line. One interesting question is whether these cell line specific pools share any models, which we addressed by calculating the intersection between Rp.1257 and Rp.1851. The result is an empty set, which means the model pools Rp.1257 and Rp.1851 are disjoint. In the next step, we wanted to further characterize and explore these cell line specific pools for information on crosstalk and Sorafenib targets.

#### 3.4.2. Classification Shows Differences Between Cell Lines

Besides the sizes of the pools and the information about the intersection of subpools, we did not receive any information about the models within a pool yet. Since we were interested in the structure of the models, especially the wiring of Sorafenib and crosstalk edges, we selected the classification analysis from the pipeline. Here, we classified for the number of optional edges and the presence of an optional edges. As a result, all models within one class have the same interaction graph, thus only differ in their logical functions (see tables in online repository). Looking at these classes, we can state that for both specific pools there are no rejected edges, since each edge appears in at least one model. This also means that we cannot exclude any of the Sorafenib targets for both cell lines. We can visualize these results by showing the frequency of edges across a pool in **Figure 5**, where we define the frequency of an edge as the number of models containing this edge divided by the number of models in the pool. Here, the graphs (A) and (B) for the full pools of each cell line show differences in the frequency of the optional edges, especially in the Sorafenib targets and the feedback. While in MZ1257RC the frequency of Sora influences is high, in MZ1851RC they are low especially for IGFR and EGFR. Moreover, one of the objectives was to identify the effect of the mutation in mTOR on the feedback in cell line MZ1851RC, where in **Figure 5B** 100% of the models in the pool Rp.1857 contain this feedback, while in the other cell line this value only reaches 71%. Also we can observe that every optional edge is present, since no edge is missing.

We can also restrict the classification to a subset of models, e.g., shown in the **Figures 5C,D**. For cell line MZ1257RC, we selected all models that contain an activating edge from Sorafenib to Raf, where we can observe an enrichment in the crosstalk from PI3K to Raf compared to the full pool (**Figure 5C**). In contrast, if we filter for all models with an edge from Sora to EGFR in cell line MZ1851RC (**Figure 5D**), the connection between Sora and IGFR is lost, which means that there is no model containing both edges. Moreover, the frequency of the connection between EGFR and PI3K is reduced. However, since these effects are statistics across a pool of models, it is hard to draw any conclusions about single models.

#### 3.4.3. Minimal Mechanisms for Sorafenib Targets

Although the classification analysis provides a good overview and intuition about the cell line specific pools, the result shows that a more detailed view can provide more information. Looking at

the minimal structures or mechanisms of each pool, we wanted to extract more information on how the crosstalk might be linked to the Sorafenib mechanism. For this aim, we analyzed the cell line specific pools for two features: the number of Sorafenib targets and possible crosstalk mechanisms. Due to the large number of models in the pools, we separated the pool for three scenarios: no influence of Sorafenib, meaning that all three optional outgoing edges of Sora are not present, Sorafenib has one target only, Sorafenib has exactly two targets and Sorafenib has exactly three targets. In **Table 3**, the minimal models according to these scenarios for the pool Rp.1257 and in **Table 4** for the pool Rp.1851 are listed.

The specific pool for MZ1257RC shows that every combination of Sorafenib target from none to all is present in the pool, thus we cannot exclude any hypotheses in this cell line (see **Table 3**). However, we can see that every model contains at least one crosstalk edge and every edge appears in a minimal model. For one Sora target, there is always a basic crosstalk from the MAPK pathway to PI3K signaling either by EGFR on PI3K or by Erk on mTORC1. For Raf, additionally the crosstalk from PI3K on Raf becomes necessary, which we already identified in **Figure 5C**. For dual targets in **Table 3C**, IGFR/EGFR requires the cross-activation from EGFR on PI3K or the feedback, IGFR/Raf require (PI3K,Raf) in combination with any of the other crosstalk or the feedback and EGFR/Raf needs (PI3K,Raf) and one crosstalk. In case all three targets are affected by Sora, the (PI3K,Raf) edge and any of the crosstalks or the feedback are required.

The minimal structures in the second cell line MZ1851RC exclude two scenarios: IGFR/EGFR as dual targets (as shown in

Minimal models after classification of Rp.1257 for (A) no Sorafenib targets, (B) exactly one target, (C) exactly two targets, and (D) all three possible targets are affected.

TABLE 4 | Minimal mechanisms for Sorafenib targets and crosstalk in Rp.1851.

Minimal models after classification of Rp.1851 for (A) no Sorafenib targets, (B) exactly one target, (C) exactly two targets, and (D) all three possible targets are affected.

**Figure 5D**) and all three targets simultaneously. **Table 4** shows models in the pool that are not affected by Sorafenib. Compared to cell line MZ1257RC, there are similarities and difference in the model structures. Raf as a Sorafenib target again requires the edge from PI3K on Raf to be present and also the models always require a crosstalk from the MAPK pathway on PI3K signaling, but additionally the feedback is essential for every model.

#### 3.4.4. Interpretation of Analysis Results

The minimal models give an overview about how the system could compensate the influence of the inhibitor to fit the data for different levels of influence. For this aspect the cell lines show similarities and differences in their model structures, where two trends can be extracted from the minimal mechanisms. First, all models require at least one crosstalk edge to be able to produce trajectories that match the data we applied. Interestingly, adding more Sorafenib targets most often does not enforce more or different crosstalk edges, with the exception of Raf. Within Rp.1851 the mechanisms for every Sorafenib target and all combinations show the identical minimal mechanism, plus the edge for Raf. A possible explanation for this observation could be the symmetrical structure of the model, in particular when the feedback is active as it is in every model of Rp.1851. Both pathways consist of a cascade of activating edges with a negative feedback on the Sorafenib target. It would be interesting to apply data which breaks with this symmetry, e.g., with a PI3K inhibitor to block the crosstalk.

The second clear trend we can identify from the results is that Raf as Sora target requires the cross-activation from PI3K on Raf. Since Raf is the designated Sorafenib target, this result is interesting. Looking at the PKN structure and the data, we can see that Erk becomes active under Sorafenib treatment and the only activator for Erk is Raf. In the MAPK pathway, Raf is activated by EGFR, which itself is inhibited by Erk. Thus, if Erk should become and stay active over longer time periods as shown in the data, Raf needs another activator to compensate the inhibition through Sorafenib. However, Sorafenib was described to have a paradoxical effect on the MAPK pathway. While, in cell lines carrying a BRAF mutation the signaling was efficiently blocked, cell lines with WT-BRAF showed an activation of Erk (Hatzivassiliou et al., 2010; Heidorn et al., 2010; Poulikakos et al., 2010). Thus, further investigations are necessary to exploit whether this observation is an artifact of the model or has biological relevance. In detail, paradoxical activation by Sorafenib and the role the crosstalk from PI3K to Raf would need to be examined, which would require a refined model where the edge from Sorafenib on Raf could also be activating and more data, e.g., an experiment with a PI3K inhibitor would be interesting in this context.

#### **<IP_ADDRESS>. Overlap of Sorafenib targets**

Another general question is, whether we assume Sorafenib to have the same targets in both cell lines. One could argue that the cell lines could differ in their internal wiring meaning the crosstalk, but the biochemical targets of Sorafenib should be independent of cell lines. Assuming that all three targets, IGFR, EGFR, and Raf, are expressed in both cell lines, the intersection of the results in **Tables 3**, **4** would further narrow down possible targets. In that case, we could exclude the case of Sorafenib affecting all targets simultaneously, since in cell line MZ1851RC there are no models that have IGFR, EGFR, and Raf as targets. Moreover, the combination IGFR/EGFR is not present in Rp.1851, thus either Sorafenib targets either one of the receptors by themselves or additionally Raf in these cell lines. Even though these results are not clear, they can support and guide further studies, especially experiments where receptors are stimulated additionally to the drug treatment would be beneficial.

#### **<IP_ADDRESS>. Models without Sorafenib targets**

A surprising result of the analysis is the presence of models without a Sorafenib target. In cell line MZ1257RC, <1% of the models have no Sorafenib target, while the pool for MZ1851RC 16% of the models fall into this category. Since the data clearly shows an effect of the drug on components in this pathway, we expected all models to have at least one target of Sorafenib to be influenced. Thus, the data set from cell line MZ1851RC seems to be not restrictive enough for every model to require an interaction from Sorafenib. Since only a subset of components is measured, some models can match the data by specific initial states. Here, additional data would be beneficial to refine the results, especially measuring more components would reduce the degree of freedom for fitting the data.

# 4. DISCUSSION

In this paper, we present an alternative approach to standard modeling procedures. Instead of building and validating one model, we incorporate uncertain information or hypotheses to build a pool of models that is then filtered for data and analyzed using specific strategies. An advantage of this method is that we can test multiple hypotheses at the same time, but it comes at the cost of high complexity and challenging analysis. For this reason, we created a pipeline with specifically defined objectives and analysis templates that the modeler can select and combine. In addition to templates for objectives, data formalization and pool analysis presented in previous work (Thobe et al., 2014, 2017; Streck et al., 2015), we introduce two new objectives, namely finding driver mutations and drug testing, as well as two analysis options, namely classification and minimal models.

In the second part of the paper, the pipeline is applied to study the uncertain wiring and effect of a drug in cancer cells based on a rich data set. Two RCC cell lines, MZ1257RC and MZ1851RC, were observed to behave differently upon Sorafenib treatment, thus we tested possible drug targets in the MAPK and PI3K signaling and also investigated possible crosstalk between these pathways in a cell line specific manner, incorporating a mutation with uncertain effect as objectives. As a result, a substantial reduction from 19,404 for the initial pool to 1,017 for MZ1257RC and 293 for MZ1851RC was observable, and the empty intersection of both pools shows that the cell line specific models indeed have a different wiring. In order to cope with the complexity of having hundreds of models as outcome of the study, we developed different analysis tools. Here, we showed that classification of the pool can provide an overview on the models in the pool and give information on essential or neglected edges. In the case study, the classification showed that the feedback from mTORC1 on IGFR was active in both cell lines. We had set this edge to optional since the cell line MZ1851RC carries a mutation in mTOR and we hypothesized that this affects the feedback. As a result, all models in Rp.1851 show the feedback in their models and thus the mutation does not affect the function of mTORC1 toward IGFR. However, for cell line MZ1257RC, which does not carry a mutation in mTOR, only 71% of the models in the pool contain this edge. An explanation for this could be that we had to exclude the data for IGFR in the Bioplex experiment, since the variance of the data was too low to allow for a meaningful discretization, which could also be the reason for the larger model pool in comparison to MZ1851RC.

For the crosstalk and the Sorafenib mode of action the classification analysis showed no clear trend, since the results are complex and hard to interpret. For this reason, we listed the minimal models according to the number of Sorafenib targets and the required crosstalk to gain more detailed information on the simplest solution (**Tables 3**, **4**). Even though we cannot exclude any Sorafenib target and crosstalk in the analysis, we are able to identify patterns, where specific Sorafenib targets require different crosstalk edges to be present, e.g., Raf requires a crosstalk from PI3K on Raf. Another important observation from the classification was that there are models in the pools of both cell lines without any interaction between Sorafenib and its target. The conclusion from this result is that the data was not restrictive enough to exclude these models and further data is necessary resolve this issue. However, in case we would only fit one model to the data, we would have missed this lack of expressiveness.

The strength of underlying approach is based on its paradigm of considering possibly huge sets of models for testing and comparison. Consequently, it does not scale as well as single model approaches. The software utilized here is limited by its model checking tool NuSMV, where more than 50 components are not solvable within reasonable time. For the analysis presented in this paper, the program was run on a Ubuntu 17.10 workstation with a processor i7-7700, 3.6 GHz, and 32GB RAM. The script with 10 components took 143 minutes in TomClass, which included building the pool, model-checking, and classification of 19,404 models. Tools like caspo list running times of approximately 56 mins for models of size 45 generating a model pool with 384 models and thus can still handle medium sized models (Videla et al., 2017). The software Tremppi was shown to be able to handle a model pool of size 259,200 and perform model-checking of 40 data sets within 151–177 min depending on parameter settings on a similar workstation (Streck et al., 2016). In general, the kind of models that are feasible for this approach are a trade-off between number of components and number of uncertain edges, the latter of which affects the size of the model pool. This approach aims at exploring uncertainties in small to medium sized models, which is well-suited to represent interesting processes such as signaling pathways and regulatory modules.

While the generation and especially model-checking process is computationally expensive, the analysis and interpretation of these pools, which in our case are just large tables, is challenging from a biological perspective. Thus we propose to define clear objectives for designing the study as well as offer different analysis strategies to extract new information from the complex results. There are many possibilities for extensions especially one could think of further biologically interesting objectives, but also including different kind of data and more analysis options such as algorithms to find special patterns. Moreover, we are not limited to the Boolean setup, but are able to handle multivalued models (Streck et al., 2015). Finally, although the pipeline was developed for signaling networks, the approach can be applied to any related modeling problem.

#### DATA AND SOFTWARE AVAILABILITY

Python scripts, data sets and the classification of the model pools from the case study analysis are available on GitHub: https:// github.com/kthobe/RCC\_ModelPoolAnalysis.

The software TomClass is also available on GitHub: https:// github.com/hklarner/TomClass.

#### AUTHOR CONTRIBUTIONS

KT designed study, analyzed and interpreted data, implemented computations, and drafted the manuscript. CK designed and performed experiments, processed and interpreted data, revised the manuscript. CS participated in discussion, interpreted data, and revised the manuscript. HS participated in conceptualization and discussion, interpreted data, reviewed

#### REFERENCES

and revised manuscript. All authors read and approved the final manuscript.

# FUNDING

The work was partially funded by the German Federal Ministry of Education and Research (BMBF), grant no. 0316195.

# ACKNOWLEDGMENTS

This work was based on the dissertation of Thobe (2017). The authors would like to thank Hannes Klarner and Adam Streck for technical support.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2018.01335/full#supplementary-material

Supplementary Table 1 | Dataset from two Bioplex experiments. Each experiment was done for each cell line, where MZ1851RC is shown in the top table and MZ1257RC in the bottom. Both cell lines were treated with DMSO or Sorafenib at the indicated time before measurement and all measured at once, thus one experiment can be seen as one batch. Measurements marked in red were excluded from the analysis due to low variance.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Thobe, Kuznia, Sers and Siebert. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# BioLQM: A Java Toolkit for the Manipulation and Conversion of Logical Qualitative Models of Biological Networks

#### Aurélien Naldi\*

Computational Systems Biology Team, Institut de Biologie de l'École Normale Supérieure, École Normale Supérieure, CNRS, INSERM, PSL Université, Paris, France

Here we introduce bioLQM, a new Java software toolkit for the conversion, modification, and analysis of Logical Qualitative Models of biological regulatory networks. BioLQM provides core modeling operations as building blocks for the development of integrated modeling software, or for the assembly of heterogeneous analysis workflows involving several complementary tools. Based on the definition of multi-valued logical models, bioLQM implements import and export facilities, notably for the recent SBML qual exchange format, as well as for formats used by several popular tools, facilitating the design of workflows combining these tools. Model modifications enable the definition of various perturbations, as well as model reduction, easing the analysis of large models. Another modification enables the study of multi-valued models with tools limited to the Boolean case. Finally, bioLQM provides a framework for the development of novel analysis tools. The current version implements various updating modes for model simulation (notably synchronous, asynchronous, and random asynchronous), as well as some static analysis features for the identification of attractors. The bioLQM software can be integrated into analysis workflows through command line and scripting interfaces. As a Java library, it further provides core data structures to the GINsim and EpiLog interactive tools, which supply graphical interfaces and additional analysis methods for cellular and multi-cellular qualitative models.

Keywords: qualitative modeling, computational systems biology, biological networks, boolean networks, static analysis, model conversion

# 1. INTRODUCTION

Logical models are highly abstract dynamical models, which have been proposed to study biological regulatory systems in the late 60s (Kauffman, 1969; Thomas, 1973). This modeling framework has since gained popularity (Bornholdt, 2005; Saadatpour and Albert, 2013; Samaga and Klamt, 2013) and has been successfully applied to a wide range of regulatory and signaling systems (Saez-Rodriguez et al., 2007; Naldi et al., 2010; Helikar et al., 2013; Abou-Jaoudé et al., 2016).

In logical models, components are represented by discrete variables with a small range of possible values, representing qualitative differences in activity. Boolean components can only be active (1) or inactive (0), while multi-valued components define multiple activity levels. Regulatory effects are often represented as signed arcs between components in the **regulatory graph**. These effects are further formalized as logical rules (also called logical parameters or logical functions),

#### Edited by:

Yoram Vodovotz, University of Pittsburgh, United States

#### Reviewed by:

Noriko Hiroi, Keio University, Japan David McMillen, University of Toronto Mississauga, Canada

\*Correspondence:

Aurélien Naldi<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology

Received: 04 April 2018 Accepted: 25 October 2018 Published: 19 November 2018

#### Citation:

Naldi A (2018) BioLQM: A Java Toolkit for the Manipulation and Conversion of Logical Qualitative Models of Biological Networks. Front. Physiol. 9:1605. doi: 10.3389/fphys.2018.01605

**265**

specifying the target activity level of each component according to the current levels of its regulators (a subset of all model components). Interactive software for model definition such as GINsim (Naldi et al., 2018a) or The Cell Collective (Helikar et al., 2012) enable the definition of regulatory graphs and logical rules. However, these logical rules are self-contained and can be used to recover signed regulatory interactions. The relative simplicity of this formalism enables the definition of large models with dozens of components, without requiring precise knowledge of kinetic parameters. A formal definition of logical qualitative models is provided in **Appendix 1** in Supplementary Material.

The CoLoMoTo consortium was recently founded to facilitate model sharing and foster cooperation in the qualitative modeling community, building on the introduction of the SBML qual exchange format (Chaouiya et al., 2013; Naldi et al., 2015). The bioLQM toolkit presented here reinforces this effort by implementing a collection of model modification, format conversion, and dynamical analysis operations in an extensible architecture illustrated in **Figure 1**. On one hand, format conversions enable the integration of several software tools in complex analysis workflows. On the other hand, the core data structure and model modifications provide building blocks for the development of integrated modeling tools, which can add their own model edition and visualization capabilities. BioLQM is notably embed in the popular GINsim software (Naldi et al., 2018a), which provides a graphical interface to most of its features. It is also used as backend for model definition and computation of successor states in Epilog (Varela et al., 2013), as well as in the CoLoMoTo notebook for model conversion and some dynamical analysis features (Naldi et al., 2018b). Preliminary versions of this toolkit were mentioned as the "LogicalModel" library Chaouiya et al. (2013) and Naldi et al. (2015).

Section 2 introduces model loading, saving and converting operations. Section 3 introduces the simulation and dynamical analysis features. Section 4 introduces model modifications. Section 5 illustrates the use of these features through the command-line and scripting interfaces for the analysis of a small model of the p53-Mdm2 network controlling DNA repair.

#### 2. LOADING AND CONVERTING LOGICAL QUALITATIVE MODELS

The increasing use of qualitative models to study biological systems led to the development of various software tools for the logical formalism (Albert et al., 2008; Garg et al., 2008; Müssel et al., 2010; Terfve et al., 2012; Naldi et al., 2018a) and related qualitative approaches (Batt et al., 2012; Paulevé, 2017; Stoll et al., 2017). Most software tools use their own file format for the definition of models, hindering the delineation of analysis workflows combining different tools. The SBML qual exchange format (Chaouiya et al., 2013) has recently been proposed to improve interoperability between modeling tools. However SBML support is often missing from existing software and may not be a priority for newer ones.

To ease model exchange between software tools that do not all support the SBML qual format, the bioLQM toolkit provides an extensible list of format handlers connected to the internal model representation. Each format is described as a Java class providing annotations (name of the format, default file extension and multi-valued support) along with optional implementations of model import (loading a file into the internal representation) or export (saving the internal representation to a file) operations. These descriptor classes are available through service discovery to facilitate the addition of new formats.

The supported formats are listed in **Table 1** and in bioLQM documentation 1 . BioLQM uses JSBML (Rodriguez et al., 2015) to load and save SBML qual models. The other import parsers are based on the antlr parser generator (Parr and Quong, 1995). While some formats natively support multi-valued models, many are limited to the Boolean case. Multi-valued models can be exported to these Boolean formats through an implicit booleanization step, described in section 4.

#### 3. MODEL DYNAMICS AND SIMULATION

A **state** of a model is a vector giving the activity levels of all its components. As the activity level of each component is restricted to a finite range, the **state space** (containing all possible states) itself is also finite. However, the total number of possible states grows exponentially with the number of components. We say that a component is **called to update** in a given state if the evaluation of the associated logical rule is different from its current activity level: for example an inactive component can become active. **Stable states** (also called fixed points, or steady states) are states in which no component is called to update. Such stable states denote a qualitative equilibrium in which all components can maintain their current activity level.

The dynamics of the model (i.e., its evolution over time) is given by transitions between states of the model, controlled by the updating calls (i.e., by the logical rules of the model) and by **updating modes** which define the synchronization between concurrent updating calls. Various types of updating modes have been introduced, with most software tools focusing on a specific subset. BioLQM aims to provide an extensive choice of updating modes in a single toolkit. In the following subsections, we further distinguish deterministic and non-deterministic simulations and provide an overview of all updating modes implemented in bioLQM. While stable states, which have no transition toward other states, do not depend on the updating mode, reachability properties and cyclical attractors can be strongly affected by the choice of updating mode as illustrated in **Figure 2**. More formal definitions of the updating calls and updating modes are given in Appendix 2 in **Supplementary Material**.

#### 3.1. Deterministic Simulations

In a **deterministic** simulation, each state has a unique successor, except stable states which have no successor at all as we consider here that a successor must denote a change of state. Starting with an initial state, a deterministic simulation yields an ordered list

1 See http://colomoto.org/biolqm/

model modifiers to transform an input model into a modified model ; (iii) the tools module contains a collection of analysis tools. All these feature are accessible through a central service manager, which handles service discovery and serves as main entry point for the Java API. A simple command line launcher provides quick execution of simple workflows, while a scripting engine can be used for more complex use cases.

The Import/Export capabilities are listed in the corresponding columns (all formats can be exported). The formats natively supporting multi-valued models are also identified, other formats rely on implicit model booleanization.

of successive states, called a **trace**. Given a sufficient number of steps, all traces end in an **attractor**, which can be either a stable state or a **cyclical attractor** of length k in which the k-th successor of each state is itself. The trace tool, illustrated in section 5, uses an initial state and a deterministic updater to compute a simulation trace. The following deterministic updating modes are supported:

FIGURE 2 | Comparison of updating modes. State transition graphs obtained with the multi-valued model shown in the top part using various deterministic (left-side) and non-deterministic (right-side) updatings. Dashed arcs denote multiple transitions and node coloring emphasizes attractors. Note that the stable state is common to all updating modes. The sequential and priority updaters follow the implicit ordering of the components. The STG obtained with the complete updating contains all synchronous and asynchronous transitions, as well as additional transitions to leave the states encompassing more than two updating calls. These transitions are colored in red in the corresponding panel. Finally, the bottom-right panel contains the asynchronous STG obtained for the booleanized version of the model. In this STG, gray nodes and arcs on the left side correspond to non-admissible states and transitions between them. These states are unreachable from the admissible ones, and the transitions enabling to leave this set of states are highlighted.

(Robert, 1986). The definition of a block-sequential updater relies on an ordered partition of the model components.

• The **synchronous priority** updating is also based on a partition of components into blocks, but only the first block containing updated components will be considered. The set of possible updaters is a subset of the priority-based updaters introduced by Fauré et al. (2006).

#### 3.2. Non-deterministic Simulations

In a **non-deterministic** simulation, each state can have several successors. Starting with an initial state, a non-deterministic simulation can lead to a large number of alternative trajectories. This type of dynamics if often represented as a **State Transition Graph** (STG), where the nodes are states of the model, and arcs denote possible transitions between these states. Like in the deterministic case, all trajectories end in an attractor, but starting from an initial state, a non-deterministic simulation can lead to several alternative attractors. These attractors can be **stable states**, **cyclical attractors**, as well as sets of intertwined cycles called **complex attractors**. More formally, all attractors are terminal strongly connected components of the STG. State transition graphs can represent deterministic traces as well as more complex dynamical behaviors. Such a graph can cover several alternative initial states or even all possible states. The current version of bioLQM supports the definition of nondeterministic updaters, enabling the computation of the lists of successor states. However, it does not provide a complete engine for non-deterministic simulations, or a data structure for state transition graphs. GINsim (Naldi et al., 2018a) implements these features on top of bioLQM. The following non-deterministic updating modes are supported:

#### 3.3. Stochastic Simulations

Stochastic updaters enable the computation of a single successor, which is selected randomly among multiple possibilities and can thus change between calls. A stochastic updater can be derived from any non-deterministic updater by assigning identical probabilities to all transitions defined by the original updater. Alternatively, a custom updater can be constructed by defining individual probabilities.

BioLQM provides the random tool to compute single random trajectories using the above stochastic updaters. This tool is limited to the construction of individual trajectories and does not provide a complete stochastic analysis. As listed in **Table 1**, bioLQM enables the conversion of Boolean models to the format of the MaBoSS software, which uses the Gillespie algorithm to estimate the probabilities of Boolean states of a continuous time Markov process, and provides a collection of scripts to further analyze the simulation results (Stoll et al., 2017).

# 3.4. Identification of Attractors

The dynamical analysis of large regulatory networks through model simulation suffers from combinatorial explosion, especially in the non-deterministic case. BioLQM implements two published methods based on constraint-solving for the identification of attractors without explicit state enumeration.

While complex attractors are well estimated through stable patterns, their exact identification requires further analysis using external software tools, adapted to the selected updating mode. In the synchronous case, the BNS tool (Dubrova and Teslenko, 2011) identifies cyclical attractors of length k using constraint solving. This approach could be extended to other deterministic updatings, but can not handle non-deterministic cases. In contrast, BoolSim uses symbolic exploration for the identification of complex attractors in the synchronous and asynchronous case (Garg et al., 2008). While this approach scales better than simple simulation, it is more sensitive to combinatorial explosion than approaches based on constraint-solving. To perform the analysis provided by the BoolSim and BNS tools, bioLQM can convert models to their respective formats.

# 4. MODEL MODIFICATIONS

Several software tools propose to emulate biological **mutations** by constructing model variants in which one or several logical rules have been modified. In bioLQM, the various **model modification** tools enable the flexible definition of model variants. The resulting modified models can have a different set of components than the original model. Each modification can be described by a keyword (identifier of the type of modification) and some parameters. The model modifier API in bioLQM allows to chain several modifications before model conversion or analysis. The following describes the various types of model modifications implemented in bioLQM.

#### Perturbations

A **perturbation** (often called **mutation**) enables to change some of the logical rules of a model. BioLQM provides three types of "atomic perturbations" (fixed value, range restriction, and removal of a regulator) which modify a single logical rule. They are briefly described below, more formal definitions can be found in Appendix 3 in **Supplementary Material**. "Multiple perturbations" can then be used to combine several atomic perturbations. The definition of these perturbations is supported by a simple syntax, as illustrated in section 5.3 and described in the online documentation.

Perturbations are commonly used to model gene knockouts by **fixing the activity level** of the corresponding component to 0, or ectopic expressions by fixing it to 1. Multi-valued components can also be fixed to a higher activity level (inside their normal activity range).

**Restricting the activity range** of multi-valued components enables to account for a partially impaired activity ([loss of the higher activity level(s)] or to set a minimal activity level.

Lastly, it is possible to define the **perturbation of a single interaction**, i.e., to remove one of the regulators of a component. This type of perturbation enables for example the definition of the loss of a single binding site preventing the action of the source component on a subset of its targets. The removal of an interaction amounts to rewrite the logical rule of its target component. Note that the atomic perturbation describes the effect on a single target: a single "biological mutation" may correspond to a "multiple perturbation" in the model if several targets are affected by the loss of the same binding site. This type of perturbation is also convenient to evaluate the importance of an interaction representing an hypothetical effect.

# 4.1. Model Reduction

Model reduction aims to ease the analysis of models with a large number of components by constructing a smaller model involving fewer components, but exhibiting similar dynamical properties. BioLQM provides a model reduction method which updates the logical rules of the remaining components to emulate the effect of the removed components (Naldi et al., 2011; Veliz-Cuba, 2011). This reduction preserves key dynamical properties of the model, in particular the stable states and stable patterns. However, it can affect some dynamical properties, depending on the choice of reduced components.

This modifier usually relies on the specification of the set of components to reduce. Some types of reduction can be fully automated. In particular, bioLQM supports the reduction of output components, which was shown to preserve attractors and reachability properties (Naldi et al., 2012), as well as the propagation of fixed components, which has also been shown to preserve attractors (Saadatpour et al., 2013).

After reduction, the reduced components are not fully eliminated from bioLQM: they are no longer allowed to regulate other components, but they keep a logical rule to allow the computation of their expected value in the reduced model.

# 4.2. Boolean Mapping of Multi-Valued Models

As discussed above, some software tools and formats are limited to Boolean models, for example as they rely on specific theoretical results or data structures. To apply such software tools to the analysis of a multi-valued model, we can construct a Boolean model such that its dynamical properties can be transferred to the original multi-valued model.

This **model Booleanization** step is based on the Boolean mapping discussed by Didier et al. (2011). In this mapping, a multi-valued component with a maximal activity level m is replaced by m Boolean components, each denoting increasing activity. All possible states of the original model can then be associated to states of the Boolean model. The logical rules of the new model ensure that we obtain the same transitions between these states. However, some states of the Boolean model are not mapped to states of the original model. These additional states are called "non-admissible states."

The dynamical properties observed on the admissible states of the Boolean model can be transferred to the original model. The implementation proposed here further ensures that all synchronous and asynchronous simulations starting with a nonadmissible state can lead to an admissible state after a sufficient number of steps. This property ensures that no attractor contains any non-admissible state (see **Figure 2**).

Model Booleanization is used automatically when converting multi-valued models to formats supporting only Boolean models. It can also be performed explicitly, like other model modifications.

# 5. USE CASE: ANALYSIS OF THE P53-MDM2 NETWORK

The cellular response to DNA damage relies on the p53 transcription factor, which induces the synthesis of DNA repair proteins. The ubiquitin ligase Mdm2 blocks the transcriptional activity on p53 in the nucleus, while p53 activates the transcription of Mdm2 and inhibits its nuclear translocation. In this section, we use a logical model involving DNA damage, p53, Mdm2 in the cytoplasm and Mdm2 in the nucleus. See the recently published GINsim tutorial Naldi et al. (2018a) and the enclosed references for a more complete description of this system and its encoding into a logical model.

In the following, we define the model in a text file named p53.mnet, using a simple text format for the definition of multi-valued logical models (p53 and cytoplasmic Mdm2 are represented by ternary components). Each line of the file reproduced below assigns a logical function to one of the components of the model. The line starts with the the identifier of the component, separated from the function itself by a leftwards arrow (<-). The &, |, and ! symbols stand for the AND, OR, and NOT operations respectively. The colon character (:) is used to specify multi-valued thresholds, both for assigning the target component and inside the functions.

```
DNAdam <- DNAdam & !p53:2
p53:2 <- !Mdm2nuc
Mdm2cyt:1 <- !p53:2
Mdm2cyt:2 <- p53:2
Mdm2nuc <- Mdm2cyt:2 | (Mdm2cyt & !p53 & !DNAdam)
```
# 5.1. Install and Launch bioLQM

Documentation, source code and releases (under the LGPL v3 license) are available on http://colomoto.org/biolqm. BioLQM is distributed as a JAR file2 , which can be launched with the command java -jar bioLQM.jar. In this section, we will use the bioLQM command as shorthand.

## 5.2. Resting State and DNA Repair

We start by looking for the stable states (fixed points) of this model. For this, we launch bioLQM on the command-line, load the model from the p53.mnet file defined above, and run the fixpoints tool. The corresponding command line and its output are reproduced hereafter.

```
$ bioLQM p53.mnet -r fixpoints
DNAdam p53 Mdm2cyt Mdm2nuc
0011
```
In the output, bioLQM displays a first line with the list of components, followed by a line for each identified stable state, giving the activity level of each component in the same order. The p53-Mdm2 model has a single stable state corresponding to a resting state in absence of DNA damage. In this state, the basal activity of Mdm2 prevents p53 activation. This analysis shows all the stable states of the model, but does not identify more complex attractors. We can then use the trapspace tool to identify stable patterns, which provide a good approximation of complex attractors in practice.

```
$ bioLQM p53.mnet -r trapspace
DNAdam p53_b1 p53_b2 Mdm2cyt_b1 Mdm2cyt_b2 Mdm2nuc
0 0 0 1 0 1
```
In this output, the two multi-valued components of the model have been extended to four Boolean components. While this requires a careful interpretation, it provides fine-grained results for complex attractors in which multi-valued components can be restricted to a range of their possible activity levels. Here we obtain a single pattern corresponding to the previously identified stable state. Note that this result does not strictly rule out the existence of a complex attractor, but attractors which do not correspond to such stable patterns are rare in practice and often depend on subtle delay effects. In this model, the resting state is indeed the only attractor.

We can then evaluate the behavior of this network upon addition of DNA damage to this resting state. For this, we use the trace tool to perform a synchronous simulation, starting with an initial state (defined after the -i flag) obtained by adding DNA damage to the resting state.

```
$ bioLQM p53.mnet -r trace -i 1011
1011
1010
1110
1210
0220
0221
0121
0011
```
In this simulation trace, we see that the introduction of DNA damage in the resting state leads to the inactivation of Mdm2 in the nucleus, enabling the activation of p53. This triggers DNA repair and allows Mdm2 to accumulate in the cytoplasm. Finally, Mdm2 can enter the nucleus and inhibit p53, coming back to the resting state. In this simulation, we assume that all possible transitions happen synchronously in each state, which could lead to artefactual trajectories. Asynchronous simulations are widely considered as more reliable, but they lead to a large number of alternative branches and are not well suited for simple commandline simulations. We can however perform a random walk in the set of possible asynchronous trajectories using the random tool. In this case, all asynchronous trajectories eventually lead to the same stable state (not illustrated here).

# 5.3. Definition of Model Perturbation

We then apply a perturbation to study the impact of a p53 knockout on the list of stable states. The -m perturbation parameters trigger the construction of a modified model. The following parameters (up to the next flag starting with a minus sign) define the modified functions. Here p53%0 describes a loss of p53 activity.

```
$ bioLQM p53.mnet -m perturbation p53%0 -r fixpoints
DNAdam p53 Mdm2cyt Mdm2nuc
0011
1010
```
We see that the resting state is still valid in the p53 knockout, however a new stable state appears in which DNA damage could not be repaired.

Instead of a full knockout of p53, we then evaluate a more subtle perturbation in which only its ability to trigger the DNA repair machinery is impaired. This corresponds to the removal of the interaction between p53 and DNAdam in our model.

```
$ bioLQM p53.mnet -m perturbation p53:DNAdam%0 -r fixpoints
DNAdam p53 Mdm2cyt Mdm2nuc
0011
```

```
$ bioLQM p53.mnet -m perturbation p53:DNAdam%0 -r trapspace
DNAdam p53_b1 p53_b2 Mdm2cyt_b1 Mdm2cyt_b2 Mdm2nuc
0 0 0 1 0 1
1 - - 1 - -
```
Here we see that this perturbation does not affect the stability of the resting state, and does not create an additional stable state as in the full p53 knockout case. However, the trapspace tool reveals the creation of a complex attractor involving oscillations

2 It requires a Java Runtime Environment, see https://www.java.com

of p53 and Mdm2. Note that these oscillations exist transiently in the original model but lead back to the resting state after DNA repair.

# 5.4. Model Conversion Enables Interoperability

As discussed in section 2, the analysis of complex models can combine several software tools. After running the following command, the new p53.sbml file will contain the functions defined above in the SBML qual format. This format is suitable for use in several other tools, or for submission in the BioModels database (Chelliah et al., 2013).

**\$ bioLQM p53.mnet p53.sbml**

# 5.5. Definition of Complex Analysis as Scripts

More complex analysis tasks can use the integrated **scripting interface**. Based on the java scripting engine, it supports scripts written in javascript (as part of the java platform) or in another supported language by providing additional libraries (including python and lua). The following sample script generates all possible individual knockout perturbations, and saves each modified model.

This script can be launched using the -s flag, followed by the script file name. Additional arguments can be used to adapt the behavior of the script. In this example, we specify the name of the original model file.

**\$ bioLQM -s generate\_perturbations.js p53.mnet**

The recently introduced CoLoMoTo Docker image (Naldi et al., 2018b) provides a python API integrating several complementary software tools. This environment includes a dedicated python API for bioLQM, which plays a central role in model conversion.

# 6. SUMMARY AND DISCUSSION

The increasing use of logical models of biological regulatory networks led to the development of multiple complementary software tools for their analysis. The recent introduction of the SBML qual format (Chaouiya et al., 2013) and the formation of the CoLoMoTo consortium (Naldi et al., 2015) aims to facilitate the exchange of models between tools. The bioLQM toolkit enables the use of additional software tools through conversion to their native formats. It provides model conversion operations in the CoLoMoTo notebook (Naldi et al., 2018b), enabling the delineation of analysis workflow involving a series of different tools.

BioLQM can also be used to apply various perturbations to the converted models, enabling the study of model variants emulating a knockout, an ectopic activity, or the loss of an interaction. Model modifications include the booleanization of multi-valued models for analysis with tools restricted to a Boolean formalism, as well as model reduction, decreasing the number of components to ease the analysis of complex models.

Finally, bioLQM provides several internal tools for the dynamical analysis of logical models. Two of the included tools allow the construction of deterministic and stochastic simulation traces, based on a comprehensive collection of updating modes. BioLQM also implements non-deterministic updating modes, which can be used as core components of complete simulation engines, as done by the GINsim (Naldi et al., 2018a) and Epilog (Varela et al., 2013) software suites. Two other tools enable the efficient identification of stable states and the approximation of most complex attractors.

The features described above are organized in a flexible architecture to facilitate the addition of new modules (file formats, model modifications, analysis tools) and to provide a consistent API. In the next version, the configuration API of analysis tools will be further improved to improve their use through python scripts in the new CoLoMoTo notebook.

Hardware requirements strongly depend on the size and structure of models and the operations performed. The complexity of individual logical rules can be a limiting factor: components with tens of regulators could have intricate rules with high computational cost. Fortunately, such rules are seldom used in biological models. Any desktop computer should be able to load and convert most models, including large ones. However, detailed dynamical analysis of models beyond 30 components can rapidly fill the available memory. In bioLQM, the fixpoints and trapspace analysis tools rely on efficient constraint-solving methods, which can scale to hundreds of components. The trace and random simulation tools are designed to work on large models as well by avoiding to store all visited states and interrupting the simulation when a stable state is reached or after a limit on the number of steps. In future versions, these tools will further use the identified trapspaces to interrupt the simulation when reaching a complex attractor.

BioLQM uses decision diagrams to store the logical rules internally, which enforces a normalized representation of the function, depending on the ordering of components. It has the advantage of providing guarantees on the number of tests to perform to evaluate a function, but it replaces the original representation of the function, making it harder to manipulate by the user afterwards. Future version will include several alternative representations to preserve hand-crafted logical functions through conversion (when the output format allows it).

Logical models are non-deterministic when using the asynchronous updating, however individual logical rules are deterministic: they associate a single "target value" to each state of the system. The ability to lift this limitation is considered in the design of the new internal data structure, but is not an immediate goal: the next releases of bioLQM will remain focused on deterministic functions.

In logical models (and by extension in bioLQM), each component is associated to its own logical rule, however the Petri net and automata network formalisms separate components from transitions. This separation allows in particular the definition of transitions affecting several components simultaneously. Such behaviors could be emulated in logical models through the addition of synchronizing components. Proper support for this use case would require extensions of the SBML qual specification, as well as changes in the internal data structure.

Like most modeling tools, bioLQM is currently centered on logical rules, however a complete model may contain important additional information, such as annotations and graphical layout information. Model annotations are supported in SBML core (without additional extensions), however annotations can be defined in any format, hindering interoperability. Further discussions are needed within the community to delineate best practices and ensure that annotations can be shared efficiently. Graphical layout information can be stored along with SBML qual models using a dedicated extension. This information is currently supported by JSBML and GINsim, it will be integrated in future versions of bioLQM. JSON "sidecar" files could then be used to facilitate the integration of such additional information with file formats which do not support it directly.

The reproducibility of model analysis relies on sharing both the model itself and the definition of simulation parameters, in particular initial states and updating modes. A single initial state can be defined in the SBML qual file. Additional initial states and simulation parameters fall in the scope of the Simulation Experiment Description Markup Language (SED-ML) format (Bergmann et al., 2015), which does not yet support qualitative models. Ongoing discussions should lead to extensions of the SED-ML format and the Kinetic Simulation Algorithm Ontology

#### REFERENCES

(Courtot et al., 2014) to describe model modifications and simulation parameters. These extensions will then be integrated into bioLQM and other qualitative modeling software.

#### AUTHOR CONTRIBUTIONS

The author confirms being the sole contributor of this work and has approved it for publication.

#### FUNDING

AN acknowledges support from the French Agence Nationale de la Recherche (ANR), in the context of the project SCAPIN [ANR-15-CE15-0006-01].

#### ACKNOWLEDGMENTS

This work benefited from the feedback of members of the CoLoMoTo consortium (colomoto.org). Pedro Monteiro, Claudine Chaouiya and Denis Theffry provided feedback for the integration in GINsim and Epilog. Julien Dorier and Gautier Stoll helped with the BoolSim and MaBoSS formats respectively. Loic Paulevé provided feedback and implemented the Pint format. Francisco Plana implemented the block-sequential updater. Hannes Klarner implemented the CNET and bnet formats. Céline Hernandez tested the API and provided feedback.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2018.01605/full#supplementary-material

**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Naldi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Gene Regulatory Network Modeling of Macrophage Differentiation Corroborates the Continuum Hypothesis of Polarization States

Alessandro Palma1 , Abdul Salam Jarrah2 , Paolo Tieri3,4, Gianni Cesareni1,5 and Filippo Castiglione3 \*

#### Edited by:

Matteo Barberis, University of Amsterdam, Netherlands

#### Reviewed by:

Carlos Villarreal, Universidad Nacional Autónoma de México, Mexico Nathan Weinstein, Universidad Autónoma de Mexico, Mexico

> \*Correspondence: Filippo Castiglione<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Systems Biology, a section of the journal Frontiers in Physiology

Received: 31 July 2018 Accepted: 02 November 2018 Published: 27 November 2018

#### Citation:

Palma A, Jarrah AS, Tieri P, Cesareni G and Castiglione F (2018) Gene Regulatory Network Modeling of Macrophage Differentiation Corroborates the Continuum Hypothesis of Polarization States. Front. Physiol. 9:1659. doi: 10.3389/fphys.2018.01659 1 Department of Biology, University of Rome Tor Vergata, Rome, Italy, 2 Department of Mathematics and Statistics, American University of Sharjah, Sharjah, United Arab Emirates, 3 Institute for Applied Computing, National Research Council of Italy, Rome, Italy, 4 Data Science Program, Sapienza University of Rome, Rome, Italy, 5 Fondazione Santa Lucia Istituto di Ricovero e Cura a Carattere Scientifico (IRCCS), Rome, Italy

Macrophages derived from monocyte precursors undergo specific polarization processes which are influenced by the local tissue environment: classically activated (M1) macrophages, with a pro-inflammatory activity and a role of effector cells in Th1 cellular immune responses, and alternatively activated (M2) macrophages, with anti-inflammatory functions and involved in immunosuppression and tissue repair. At least three different subsets of M2 macrophages, namely, M2a, M2b, and M2c, are characterized in the literature based on their eliciting signals. The activation and polarization of macrophages is achieved through many, often intertwined, signaling pathways. To describe the logical relationships among the genes involved in macrophage polarization, we used a computational modeling methodology, namely, logical (Boolean) modeling of gene regulation. We integrated experimental data and knowledge available in the literature to construct a logical network model for the gene regulation driving macrophage polarization to the M1, M2a, M2b, and M2c phenotypes. Using the software GINsim and BoolNet, we analyzed the network dynamics under different conditions and perturbations to understand how they affect cell polarization. Dynamic simulations of the network model, enacting the most relevant biological conditions, showed coherence with the observed behavior of in vivo macrophages. The model could correctly reproduce the polarization toward the four main phenotypes as well as to several hybrid phenotypes, which are known to be experimentally associated to physiological and pathological conditions. We surmise that shifts among different phenotypes in the model mimic the hypothetical continuum of macrophage polarization, with M1 and M2 being the extremes of an uninterrupted sequence of

**275**

states. Furthermore, model simulations suggest that anti-inflammatory macrophages are resilient to shift back to the pro-inflammatory phenotype.

Keywords: macrophage, differentiation, phenotype, model, gene regulating network, polarization, immune system

# AUTHOR SUMMARY

fphys-09-01659 November 23, 2018 Time: 15:54 # 2

Macrophages are key players in the elicitation of an efficient immune response. Latest classification of macrophage functional types comprises the classically activated (M1) macrophages with a pro-inflammatory activity and the alternatively activated (M2) macrophages, with anti-inflammatory functions. The latter is further subdivided into at least three different subsets, namely, M2a, M2b, and M2c, which are characterized on the basis of distinct eliciting signals.

Accounting for the gene-related mechanisms of macrophage differentiation is a challenging task. We have used the methodology known as gene regulation network modeling on a newly constructed network of gene regulation originated from published experimental data. We have used computer simulations to explore the dynamical behavior of this network and derived conclusions about the hypothetical continuum of macrophage polarization with M1 and M2 being the extremes of an uninterrupted sequences of states. Our simulations also suggest that anti-inflammatory macrophages are resilient to shift to the pro-inflammatory phenotype.

# INTRODUCTION

Macrophages and neutrophils of the innate immune system represent the first line of defense against most common microorganisms. Indeed, macrophages can recognize and respond to a wide range of stimuli, expressing a great variety of surface and intracellular receptors that activate several signal transduction pathways and complex gene expression patterns. Macrophages respond to extracellular stimuli upon contact with different cell types via endocytic, phagocytic, and secretory functions (**Figure 1**). Their activity is modulated by contact synapsis established with proximal cellular and molecular entities, including microorganisms, chemical mediators, and other macrophages (Gordon et al., 2014).

The monocyte–macrophage differentiation pathway is known to exhibit plasticity and diversity (Mantovani et al., 2002; Bowdish et al., 2007; Gordon, 2008; Mantovani, 2008). Similar to the polarization process of helper T type 1 and 2 cells (Th1–Th2), two distinct polarized forms of macrophages have been recognized in the past: the classically activated (M1) macrophage phenotype and the alternatively activated (M2) macrophage phenotype (Biswas and Mantovani, 2010). Moreover, macrophages have also been observed in "M2-like" states, which share some features of both M1 and M2. Indeed, recent studies support the view that fully polarized macrophages (M1 and M2) are the extremes of a continuum of macrophage polarization (Mantovani, 2008). For example, various stimuli, such as immune complexes (IC) together with LPS or interleukin-1 beta (IL-1β), glucocorticoids, transforming growth factor-β (TGF-β), and interleukin-10 (IL-10), give rise to M2-like functional phenotypes that share properties with IL-4- or IL-13-activated macrophages [such as high expression of mannose receptor (MR) and IL-10, as well as TNFα, IL-1β, and IL-6] (Mantovani et al., 2004). Variations of the gene expression patterns corresponding to M1 or M2 are also found in vivo (e.g., in the placenta and embryo, and during helminthic infection, Listeria infection, obesity, and cancer) (Raes et al., 2005; Biswas et al., 2006; Kraakman et al., 2014).

The M1 and M2 phenotypes Kraakman et al., 2014 correspond to cell activation states driven by cytokines, which are typically secreted by Th1, Th2, and T-regulatory cells, but also basophils, mast cells, B lymphocytes, and eosinophils. The M1 phenotype is polarized by single or a combination of Th1 cytokines and pro-inflammatory mediators, including granulocyte-macrophage colony-stimulating factor (GM-CSF), tumor necrosis factor (TNF)-α, IL-6, IL-1β, IL-12, and various pathogen-associated molecules, such as lipopolysaccharide (LPS). By contrast, the M2 polarization is induced by macrophage colony-stimulating factor (M-CSF), IL-4 and IL-13, IC, IL-10, as well as glucocorticoid, TGFβ, and serotonin (Sang et al., 2015) (see **Table 1**).

Although there is a wealth of information about the different macrophage subsets in vitro, features such as plasticity, heterogeneity, and adaptability make them very difficult to study using conventional experimental tools. Furthermore, as many of the studies are done in different settings or for different goals, some literature reports are not conclusive and sometimes contradictory. It is not clear how robust the different macrophage subsets are to environmental changes. In particular, how does a modification of the cytokine environment affect the phenotype of macrophages? Which polarization state is most stable? Which possible gene knockouts can lead to a phenotypic change?

Macrophages polarization is essential in orchestrating the immune system response both in infectious and sterile immune settings. To shed light on this complex molecular process and address the questions above, we employed computational modeling of gene regulatory networks (GRNs) (Karlebach and Shamir, 2008).

Computational and mathematical modeling provide a means to assemble the known relevant molecules and their interactions into a network of pathways, with cross-talk between them. This allows, for examples, the test of whether the assimilated knowledge is sufficient to reproduce experimental results, and, furthermore, introduce cell-specific perturbations into the network to generate and test hypotheses in silico. For recent reviews, see Eftimie et al. (2016), Chakraborty (2017).

Computational models of GRNs have been shown to be a good approach to study how cells integrate several signals driving the cell phenotypic changes, especially for their ability to quantitatively and qualitatively describe a great variety of poorly

TABLE 1 | Summary of key molecules in macrophage polarization as taken from the literature.

Cytokines, chemokines, receptors, and genes involved in metabolism are represented for each specific macrophage phenotype (adapted from Foey, 2014). Dashes indicate missing/contrasting data in the literature.

characterized biological situations (Méndez and Mendoza, 2016). Computational models are used to describe immunological phenomena, to provide a better understanding of aspects of the immune response, and to produce outcomes coherent with available data, thus unraveling basic mechanisms of immunology and possibly leading to new hypotheses that can be tested experimentally in vivo or in vitro (Castiglione and Celada, 2015).

Discrete logical (Boolean or multi-state) models are usually the method of choice especially when the biological questions are of qualitative nature or when the available data (and knowledge) are mainly qualitative. Boolean networks and logical models have been used extensively to model many biological systems including immunological systems such as T-cell signaling and T helper cell differentiation (Naldi et al., 2010; Abou-Jaoudé et al., 2016; Méndez and Mendoza, 2016).

There are several computational models of some pathways that are involved in the pro and anti-inflammatory immune response, such as the NF-κB, TNF-α, IL-1, and IL-10 signaling pathways. Furthermore, there are computational models of T helper cell differentiation including continuous (Carbo et al., 2014), Boolean (Martinez-Sanchez et al., 2015), multistate logical (Naldi et al., 2010), and multi-scale (Santoni et al., 2008; Tieri et al., 2014). However, we are not aware of any GRN models of the molecular network describing macrophage differentiation. We have recently developed a multiscale model (Castiglione et al., 2016) of the immune response incorporating a minimalistic Boolean model of macrophages differentiation accounting for M1 and M2 polarization, but not for the subsets of M2. Maiti et al. (2014) presented an ODE model to describe the pro- and antiinflammatory signaling in macrophages toward understanding immune homeostasis.

In this paper, we present a novel logical model of the gene regulation underlying macrophage differentiation and polarization, where the regulatory interactions and logical rules are inferred from the literature. We then used the model to study the dynamical behavior of the network. The model not only was able to reproduce known experimental data but also provides the first computational evidence of the continuum hypothesis of phenotypes which was suggested by Sica and Mantovani (2012).

# MATERIALS AND METHODS

fphys-09-01659 November 23, 2018 Time: 15:54 # 4

#### Logical Models of Regulatory Networks

Gene regulatory network modeling aims at describing the way cells integrate extracellular stimuli to run cellular programs consisting of activations and inhibitions of genes (Kestler et al., 2008).

Logical network modeling was introduced by the geneticist R. Thomas (Thomas and D'Ari, 1990; Thomas and Kaufman, 2001) for the study of GRNs. Since then, they have been developed further, and have been used extensively to model many biological systems including cell-fate determination in A. thaliana (Espinosa-Soto et al., 2004; Benítez and Hejátko, 2013), E. coli metabolism (Samal and Jain, 2008), and the differentiation and plasticity of T helper cells (Naldi et al., 2010; Abou-Jaoudé et al., 2014), to name a few.

Gene regulatory networks are typically drawn from a mixture of literature, data mining and experimental data. Signal transducers, transcription factors and target genes in the activation of specific cellular programs (e.g., cell maturation or differentiation) are identified, as well as their relationships coded in terms of inhibition/activation. This data mining step produces a network (N, E) in which the nodes N are the molecules and the edges E = E− ∪ E + are the activations (edges in the set E +) and inhibitions (edges in E −) relationships. Gene activation levels (states) or molecular concentrations are represented either by a discrete and usually very small set of values (two levels, i.e., active/inactive, represents the most used one, called Boolean) or by a continuous range of activity levels. In this paper, we have used the discrete Boolean formulation.

Each node nk of the network N has a function Fk specifying how the state of that node may change in response to changes in the states of its neighbors (the nodes nj for which there exists an edge ejk ∈ E) in the network. The synchronous or asynchronous calculation of the functions F1,. . .,Fn, at each discrete step makes the network evolve from one macro-state to another. In the synchronous mode, all node states are updated at the same time, while in an asynchronous case, nodes are randomly updated at different time steps.

The Boolean model of a GRN is therefore defined as a discrete dynamical system which can then be studied for its dynamical properties. Since the space of all possible macro-states is finite, starting from any configuration, the repeated application of the functions F1,. . .,Fn, will lead the system to be in states that it has reached before. These states correspond to stable patterns of gene expression that can be reasonably regarded as real biological states characterizing a specific cellular function. Starting from any configuration and after a certain transient period, the network dynamics will either reach a state and stay there (such a state is called a steady state), or can keep cycling forever among the same set of states (such a set of states is called a limit cycle) (Guevara, 2003; Ortiz-Gutiérrez et al., 2015). The transient period before the network dynamics reaches a certain steady state or limit cycle is called the basin of attraction of that state or cycle.

The dynamics of the system is encoded by a graph, whose vertices are all configurations (states) of the network and directed edges where each such edge indicates the transition of the system from one state to the next.

We used the software GINsim (Naldi et al., 2009) for the development of the model and the analysis of the network, including the identification of all steady states (Karlebach and Shamir, 2008; Méndez and Mendoza, 2016), and the BooleanNet Python library (Albert et al., 2008) and BoolNet R library (Müssel et al., 2010) for the study of the dynamics of the system.

# RESULTS

# Molecular Basis of the Macrophage Polarization

During the inflammation process, several immune cells are involved in initiating and maintaining the inflammatory state. Macrophages, together with leukocytes, are the first cells recruited to the inflammation site. They start releasing proinflammatory cytokines (mostly IFN-γ and IL-1β), creating an inflammatory environment. The binding of those molecules to their specific receptors triggers a signal transduction cascade resulting in the release of other inflammatory molecules. This positive feedback mechanism allows the maintenance of the inflammatory state and reinforce the M1 polarized state.

The resolution of inflammation occurs by different mechanisms, such as the downregulation of pro-inflammatory molecules, the short half-life of the inflammatory mediators, and the production of anti-inflammatory molecules. In this context, macrophages are expected to switch to M2, and, consequently, produce anti-inflammatory mediators, such as IL-10, inhibiting M1-related transcriptional regulators, while a positive feedback loop provides the means to maintain their anti-inflammatory phenotype.

Interferon (IFN) receptors have multi-chain structures and interact with members of the Janus-activated kinase (JAK) family (Darnell et al., 1994). When IFN-γ binds to its cognate receptor, the activation of the receptor-associated JAKs occurs in response to rearrangement and dimerization of the receptor subunits, followed by auto-phosphorylation and activation of the associated JAKs. This process determines the activation of classical JAK–STAT (signal transducer and activator of transcription) signaling pathways, resulting in the transcription of target genes (Platanias, 2005; Mosser and Edwards, 2008; McLaren and Ramji, 2009). Among the STATs, a pivotal role is played by STAT1, which undergoes dimerization after its JAK-mediated tyrosine phosphorylation. Hence, STAT1–STAT1 homodimer binds to cis elements known as "gamma-activated sequences" (GAS) in the promoters of the genes encoding NOS2, the MHC class II transactivator (CIITA) and IL-12, among others (Darnell et al., 1994; Sadler and Williams, 2008; Lawrence and Natoli, 2011). The IFN-associated JAK/STAT pathway exerts its function in the regulation of several immune cells, including macrophages, with a great increase of IFN production, the

synthesis of several cytokines, such as interleukins IL-1β, IL-6, IL-12, IL-18, IL-23, and TNF-α, and nitric oxide (NO), as well as reactive oxide intermediates (ROI) and enzymes required for tissue remodeling.

Toll-like receptors (TLRs) mediate the immune response to a great variety of infectious agents and facilitate transcription of many pro-inflammatory genes (Sheikh et al., 2014). LPS is a component of the Gram-negative bacteria cell wall and induces expression of a wide variety of genes that constitute the innate immune response to bacterial infections. LPS signals through TLR4 on the cell surface of many cell types, including macrophages (Kawai and Akira, 2010, 2011). Signaling through TLR4 induces rapid activation of two distinct intracellular signaling pathways: one is the MyD88-dependent pathway, which leads the cascade through interferon regulatory factor (IRF)- 3, and the other is the MyD88-independent signaling pathway, which acts through TIR-domain-containing adapter-inducing interferon β (TRIF). These pathways converge to activate the transcription of NOS2; the inducible NO synthase (Kawai et al., 2001; Doyle et al., 2002).

The M1 phenotype can also result from differentiation in the presence of GM-CSF, with increased expression of IL-12 and proinflammatory cytokines, the ability to activate Th1 cell immune responses and decreased expression of IL-10 (Krausgruber et al., 2011).

M2 macrophages exhibit a functionally distinct phenotype to that of M1s, originally via the ability of IL-4 to induce MR expression, followed by IL-13, which is another Th2 cytokine. IL-4/IL-13 and TGFβ/IL-10 have been described to be associated with priming M2 macrophage subsets (M2a and M2c, respectively). The role of IL-4- and IL-13-mediated signaling in M2 macrophage polarization has been well established both in vitro and in vivo (Gordon, 2003; Martinez et al., 2009; Gordon and Martinez, 2010). Mice with a myeloid cell-specific knockout of IL-4 receptor-α (IL4Rα) were found to lack M2 macrophage development in mouse models of helminth infection and in Th2 cell-mediated inflammation, where IL-4 has a major role (Lawrence and Natoli, 2011). It is well established that IL-4 and IL-13 are associated with Th2-type responses, which have well-defined effects on macrophages, other cells and immune functions. IL-4 and IL-13 are produced particularly in allergic, cellular, and humoral responses to parasitic and extracellular pathogens. IL-4 and IL-13 upregulate expression of the MR and MHC class II molecules by macrophages, which stimulates endocytosis and antigen presentation, and they induce the expression of selective chemokines (Gordon, 2003; Gordon and Martinez, 2010). IL-4 and IL-13 act through a common receptor chain – IL-4Rα – through signal transducer and activator of transcription 6 (STAT6).

Interleukin-1 beta and IC, together with TLR4-signaling inducers (i.e., LPS), drive the macrophage to an M2b phenotype. IL-1β not only plays a pivotal role in the initiation and maintenance of the inflammatory response but also modulates immunosuppressive mechanisms through the process of macrophages endotoxin tolerance. IL-1β is also produced in response to LPS, emphasizing a collaborative interplay between M1 and M2b macrophages in eliciting and maintaining the inflammatory response (Sato et al., 2012).

Interleukin-10 acts on a distinct plasma membrane receptor to those for IL-4 and IL-13 (Riley et al., 1999; Moore et al., 2001; Deng et al., 2012), and its effects on macrophage gene expression are different, involving a more profound inhibition of a range of antigen-presenting and effector functions, together with the activation of selected genes or functions. T cells themselves are more heterogeneous than was thought originally, including not only Th0-, Th1-, and Th2-type cells but also regulatory and possibly Th3-type cells, some of which secrete TGF-β and IL-10 (Gordon, 2003). TGFβ and IL-10 have been described to be associated with priming M2-like macrophages subset polarization. TGFβ and IL-10 modulate macrophage polarization and functional plasticity to that of an M2c subset which exhibits a characteristic cytokine phenotype of IL-10hi, IL-12lo, IL-23lo, and TGFβ + which is associated with anti-inflammatory responses, scavenging, immune regulation, tissue repair, and tumor promotion. Both TGFβ and IL-10 directly suppress immune activation via the down-regulation of the expression of MHC II and pro-inflammatory cytokine production, with an indirect effect through cross-regulation of M1-derived cytokines and functionality (Gordon and Martinez, 2010; Lawrence and Natoli, 2011; Sica and Mantovani, 2012). IL-10 is a potent STAT3 dependent inhibitor of pro-inflammatory cytokine production and NO release, after challenge with LPS. IL-10-deficient mice develop widespread inflammatory cell infiltrates, including in the bowel, and transgenic animals that constitutively overexpress ILmetricconverterProductID10 in10 in macrophages suffer from septic shock and over-activity of pro-inflammatory cytokines (Lang et al., 2002b). The upregulation of expression of IL-4Rα by IL-10 correlates with increased IL-4-dependent expression of arginase-1. IL-10 also synergizes with LPS to increase the expression of arginase-2. Therefore, IL-10 increases the total level of arginases in macrophages in many ways (Lang et al., 2002a,b).

Phenotypes depending on complex regulatory logic can be effectively studied by using mathematical and computational approaches, such as GRN models.

# A Logical Network Model of Macrophage Differentiation

We have constructed a logical regulatory network model (**Figure 2** and **Supplementary File S1**) that describes macrophage polarization using experimental data and knowledge derived from literature (see **Table 2**) and a curated database of causal relationships between biological entities (Perfetto et al., 2016). The network comprises 30 nodes and 49 interactions among them. Interactions can be either positive (activations) or negative (inhibitions) (**Figure 2**). **Table 2** shows a list of the molecules, interactions, and references from the literature supporting each interaction, while **Table 3** shows logical rules for each molecule.

Nodes are of four kinds, depending on cellular location and function (**Figure 2**): seven input nodes, which represent the extracellular stimuli (IFNγ, GM-CSF, IL-1β, LPS, IC, IL-4, and IL-10), seven receptors (IFNγR, CSF2Ra, IL-1R, TLR4, FcγR, IL-4R, and IL-10R), 14 internal regulators (STAT1, STAT5, NF-κB,

#### TABLE 2 | Interactions in the macrophage polarization network.

fphys-09-01659 November 23, 2018 Time: 15:54 # 6

(Continued)

Source and target nodes are reported as well as the sign of the interaction between them (positive: source molecule activates target molecule; negative: source molecule inhibits target molecule) and the references. Each input node is annotated with an "\_e" suffix, which stands for external stimulus, as well as an "\_out" suffix which stands for output.

PPARγ, STAT6, JMJD3, STAT3, IRF3, ERK, KLF4, SOCS1, IRF4, IL1β, and IFN-β), and two main products of each distinct type of macrophage (IL-12 and IL-10). The input nodes represent the main intercellular molecular stimuli that drive macrophage polarization, as reported in the literature. Each external molecule (input) is connected to its specific receptor, and this binding elicits a signaling cascade, involving intracellular transducers and transcription factors (mostly STAT factors). Each specific transcription factor binds the promoter of a target gene, resulting in the production of IL12 or IL10 depending on the macrophage polarized form.

Interactions among nodes are derived from experimental data available in the literature as shown in **Table 2**. All interactions have been deposited in SIGNOR (Perfetto et al., 2016), a public database of causal interactions between biological entities. Each node is associated to a logical function which determines the activation level of the node based on the activation levels reached by its source nodes in the previous time step. The logical function of each node is inferred from the available literature (see **Table 3**).

The network encompasses several pathways. Different cell fates, i.e., macrophage phenotypes, are defined by steady or stable states (also called fixed point attractors) of gene expression, and described in this dynamic model as multiple, specific, and stable configurations of activated/deactivated nodes. In other words, stable states are configurations toward which the system tends to evolve, for a wide range of starting conditions. Thus, according to the network, its starting configuration, and the initial external stimuli, the pathways lead to a configuration that resembles a specific cell state in terms of the given gene expression pattern. In this regard, we assumed that the sum of the sizes of the basins of attraction of the steady states characterizes the likelihood of finding the cell in a specific differentiation state. In other words, the probability that the cell, stimulated by cytokines, will switch to the certain differentiation state is proportional to the size of the subspace of all network configurations eventually reached by the network dynamics.

Inhibitory pathways among M1 and M2 phenotype-related transcription factors are particularly interesting, because they allow a mutual exclusivity of transcription factors and, therefore, of the macrophage phenotypes, as reported in literature (Lawrence and Natoli, 2011). Notably, among the interactions describing the network and reported in tables above, the inhibition of TLR4 and NF-κB signaling by FcγR activation were added. These relationships allow the inhibition of M1 polarization in the presence of IC, that together with LPS and IL-1β, drives the otherwise absent M2b polarization.

To analyze the dynamics of the network under different conditions we used GINsim [Gene Interaction Network simulation1 ; (Chaouiya et al., 2012)], a software tool for modeling and simulation of genetic regulatory networks (Chaouiya et al., 2012). In some cases, for further confirmation or additional details, we used the BooleanNet Python (Albert et al., 2008) as well as the BoolNet R library (Müssel et al., 2010).

The fate of a macrophage strongly depends on the local biochemical microenvironment. To reproduce these different microenvironments that influence the cells, we defined a set of inputs to run the simulations. Hence, we could discriminate among steady states with a real biological meaning.

The starting expression state of the network corresponds to the naïve macrophage M0 (unstimulated/not-activated) phenotype, in which the state of each node in the network is set to "0" (i.e., low expression).

1www.ginsim.org

In our simulations, we found that our model has five sets of steady states fitting the following five specific macrophage phenotypes markers according to literature (**Figure 3**):

We computed the steady states of macrophage polarization network using a synchronous update. The system reached 1056 states, 1040 of which are steady states and 16 are cycles made

#### TABLE 3 | Boolean functions in the macrophage polarization network.

fphys-09-01659 November 23, 2018 Time: 15:54 # 9

Based on the available literature (third column), a Boolean function (second column) is associated to each target node of the network (symbols ∧, ∨, and q indicate logical operators AND, OR, NOT, respectively).

of three different states. Among the 1040 unique steady states (**Supplementary Table S1**), 228 can be mapped to the five canonical macrophage phenotypes reported via experimental studies in the literature. The frequencies of these 228 steady states are reported in **Figure 4**. The remaining steady states do not characterize the macrophage in any of the known canonical phenotypes. These states, for which there is a lack of experimental knowledge, could correspond to input conditions not existing among in vivo inflammation settings or even be artefacts of the modeling approach. Alternatively, they could correspond to hybrid phenotypes (O'Carroll et al., 2013) resembling gene expression patterns of two or more canonical phenotypes (discussed below). It is worth to note that a higher number of steady states does not imply a corresponding higher probability of polarization, since the final outcome depends on the combination of external stimuli. In other words, the number of steady states indicates the propensity of the network logic to lead the cell to the specific phenotypes yet driven by environmental cues.

The most frequent polarized state is the M2a followed by M2c and then M1. This is consistent with the pivotal role of macrophages in inflammation (M1), and in the resolution of inflammation (M2a and M2c). On the other hand, according to our analysis, M2b is the least frequent state, which might be consistent with the lack of knowledge of M2b-related pathways which is reflected in the network. This behavior of the model is consistent with observed data (Sica and Mantovani, 2012).

A closer look at the dynamics of the model (**Figure 5**) is obtained by performing several rounds of asynchronous simulations by using the BooleanNet Python library. We observed that any combination of stimuli among IFN-γ, IL-1β, LPS, and GM-CSF keep the polarization of the M1 macrophage. Once macrophages have polarized into an M1 form, the steady states are taken as initial conditions to polarize macrophages into the three different forms of M2 macrophage. IL-4 input is activated (i.e., IL-4 binding by IL-4RA) to polarize M2a macrophage, IL-10 is activated to polarize the M2c macrophages, and IC in combination with either IL-1β or LPS is activated to polarize M2b macrophages, according to the available literature on macrophage polarization stimuli (Gordon and Martinez, 2010).

The M1 polarization is simulated starting from an M0 (i.e., all non-input signal nodes set to zero) cellular environment and switching on all input nodes, as reported in literature. Following the typical cellular response to inflammation, starting from an M1-like configuration, and M2-related external stimuli (i.e., IL-4 for M2a, IL-10 for M2c and IC in combination with LPS or IL-1β for M2b macrophages), the dynamics of transcription factors and secreted molecules (i.e., IL-12 and IL-10) show the macrophage moves from pro- to antiinflammatory states, as reported in literature. The M2-related polarizations from an M0 initial state have been also performed to check the ability of the system to simulate the situation in which new monocyte-derived macrophage populations are recruited to the inflammation site during the resolution of inflammation, in addition to M2 macrophages polarized

FIGURE 3 | Gene expression markers of macrophage polarization according to literature. Each row, associated to one of M0, M1, M2a, M2b, and M2c, indicates the expression of the 10 marker genes determining the polarization fate. White dots represent inactive genes; yellow dots indicate expressed genes.

from the pro-inflammatory M1 state (see **Supplementary File S2**).

We also tested in silico the "plasticity" of the polarized phenotypes, i.e., the capability to revert the state from inflammatory to anti-inflammatory and vice versa. In order to proceed, we run a set of numerical experiments in which macrophages, starting from the four polarized states M1, M2a, M2b, and M2c, were challenged with the four characteristic stimuli (i.e., pro-M1, -M2a, -M2b, and -M2c) resulting in 16 possible couples "initial condition/stimuli." Each of those simulation settings was repeated 104 times using the asynchronous updating scheme and averages were computed. After that, we used the steady states obtained as initial states for other simulations, giving each input from the input set (see **Figure 5**).

We focused on M1-related initial states, since a normal immune response begins with an inflammation state, followed by anti-inflammatory environment settings.

With an M0 steady state as initial condition, several stimuli were applied for each simulation. To represent the M1 polarization we gave a combination of random M1-related stimuli (LPS, GM-CSF, IFN-γ, and IL-1β). The initial state for each node of the network are those related to the M0 steady state (no active nodes at all).