Constraining fuzzy logic
Fuzzy logic is a highly flexible methodology to transform linguistic observations into quantitative specification of how the output of a gate depends on the values of the inputs 
. For example, in the simplest, ‘Sugeno’ form of fuzzy logic, one specifies the following: ‘membership functions’ designating a variable number of discrete categories (“low, medium, high', etc.) as well as what quantitative value of a particular input belongs either wholly or partially to these categories; ‘rules’ designating the logical relationships between the gate inputs and outputs; AND and OR ‘methods’ designating the mathematical execution of each logical relationship; ‘weights’ designating the credence given any rule; and ‘defuzzification’ designating a scheme for determining a final output value from the evaluation of multiple rules 
. This flexibility is important in industrial process control 
, which aims to use uncertain and subjective linguistic terms to predict how a controller should modulate a process variable to achieve the desired output.
However, our goal is to train models on quantitative biological data that are inevitably incomplete in the sense that (i) measurements are not obtained under all possible conditions and (ii) available data are not sufficient to constrain both the topology and quantitative parameters of the underlying networks. Accordingly, we sought to develop a fuzzy logic system that minimizes the number of parameters to avoid over-fitting and simplifies the logic structure to facilitate model interpretability. Because we aim to represent relationships among proteins in enzymatic cascades, mathematical relationships should be biologically relevant. We therefore use a simple Sugeno fuzzy logic gate with a defined form (see Text S1
) based on transfer functions (mathematical functions describing the relationship between input and output node values) that approximate the Hill functions of classical enzymology.
Our ‘constrained’ fuzzy logic (cFL) framework uses a simplified fuzzy logic gate that is best described by the mathematical representation in . The value of an output node of a one-input positive interaction is evaluated using a transfer function. In this paper ‘input-output’ refers to the nodes of a specific cFL logic gate, where ‘nodes’ are molecular species. We use the terms ‘model inputs’ and ‘model outputs’ to denote the overall relationship between model inputs such as ligand stimulation of cells and the collective output of the network (protein modifications or phenotypic states in our application). The transfer function underlying cFL gates is a normalized Hill function with two parameters: (1) the Hill coefficient, n
, which determines the sharpness of the sigmoidal transition between high and low output node values and (2) the sensitivity parameter, k
, which determines the midpoint of the function (corresponding to the EC50
value in a dose-response curve, ). A negative interaction is represented similarly, except that the transfer function is subtracted from one, effectively inverting it (). Varying these parameters allows us to create a range of input-output transfer functions including linear, sigmoidal and step-like (). Moreover, this transfer function is biologically relevant: protein-protein interactions and enzymatic reactions can be described by Hill function formulations to a good approximation 
Construction of gates with constrained fuzzy logic (cFL).
In some cases, use of a normalized function is too restrictive for practical application. For example, if model inputs are purely binary (values of either zero or one), the output of a normalized function would also be zero or one, making it impossible for a cFL gate to achieve intermediate states of activation. Accordingly, our cFL method allows for alternative transfer functions. For example, although the method is not limited to binary model inputs, the ligand inputs of our current work are binary (either present or not). If we used normalized transfer functions to relate these model inputs to downstream outputs, all model species would also be either zero or one. Thus, for these transfer functions, we used a constant multiplied by the binary ligand input value (see Materials & Methods
If more than one input node influences an output node, this relationship is categorized as either an “AND” or “OR” interaction. An AND gate is used when both input nodes must be active to activate the output node, whereas an OR gate is used when either input node must be active. Mathematically, we represent AND behavior by evaluating each input-output transfer function and selecting the minimal possible output node value (i.e.
, applying the “min” operator, ) whereas we select the maximal value (“max” operator; ) to evaluate an OR gate. Finally, if both AND and OR gates are used to relate input nodes to an output node, our formalism evaluates all AND gates prior to OR gates. This order of operations corresponds to the disjunctive normal or sum of products form 
Use of cFL to understand experimental data in the context of a prior knowledge network: CellNOpt-cFL
The process of training a cFL network (CellNOpt-cFL) has two starting requirements. The first is a prior knowledge network (‘PKN’; , box A). A PKN depicts interactions among the nodes as a signed, directed graph (such as a PSN) and can be obtained directly from the literature. Alternatively, a large number of commercial (e.g.
, Ingenuity Systems: www.ingenuity.com
; GeneGo: www.genego.com
) or academic (e.g.
, Pathway Commons: www.pathwaycommons.org
, reviewed in 
) pathway databases as well as integrative tools (e.g. 
) can be utilized to construct a PKN. The second requirement is a dataset describing experimental measurements characterizing node activities following stimulation of and/or perturbations in upstream nodes (ligand and inhibitor treatment in our example; , box B). CellNOpt-cFL is then used to systematically and quantitatively compare the hypothesized PKN to the experimental dataset.
CellNOpt–cFL workflow and application to toy model.
In practice, available experimental data is usually insufficient to fully constrain both the parameters and topology of the cFL models, and CellNOpt-cFL recovers many models that describe the data equally well. Due to this typical absence of firm structural and parametric identifiability 
, we examine families of models that fit the data equally well rather than attempting to identify a single global best fit. Specifically, we examine interactions in the PKN that were either retained or consistently removed by training. We also use individual models to predict input-output characteristics. This treatment allows us to calculate both an average prediction as well as a standard deviation, which we show below can be useful for discrediting inaccurate predictions.
Our method comprises three main stages (): first, structure processing converts a PKN into a cFL model; second, model training trains the model to experimental data; and third, model reduction and refinement simplifies trained models. To illustrate CellNOpt-cFL, we examine a simple toy problem of training a PKN of the phospho-protein signaling network response to TGFα and TNFα () to in silico data of activation of several downstream kinases in response to these ligands in the presence or absence of PI3K or MEK inhibition ().
In the first step, we streamline the network to contain only measured and perturbed nodes as well as any other nodes necessary to preserve logical consistency between those that were measured or perturbed (
; , Step 1), resulting in a compressed
PKN ( box C). In our example, many nodes that were in the original PKN were neither measured nor perturbed experimentally. Because these nodes could be removed without causing logical inconsistencies, they were not explicitly included in the compressed network ().
In the second step, we expand the network into the multiple logical relationships (combinations of AND and OR gates) that can relate output nodes to their input nodes (, Step 2). For example, our toy PKN was expanded to include all possible two-input AND gates governing the response of nodes with more than one possible input node ().
In the third step, we train the cFL models to the data (, Step 3). We start by limiting the possible parameter combinations to a subset of discrete parameter values that specify seven allowed transfer functions as well as the possibility that the input does not affect the output node (i.e.
the cFL gate is not present). A discrete genetic algorithm determines transfer functions and a network topology that fit the data well by minimizing the mean squared error (MSE, defined in Materials & Methods
) with respect to the experimental data.
Due to the stochastic nature of genetic algorithms, multiple optimization runs return models with slightly different topologies and transfer function parameters that result in a range of MSEs. Models with an MSE significantly higher than the best models are simply eliminated from further consideration. Models with similar MSEs but different topology and parameters result from the insufficiency of the data to constrain the model such that each model fits the data well albeit with slightly different features. We consider each individual in this group as a viable model, and all are included for subsequent analysis. Thus, after multiple independent optimization runs using the discrete genetic algorithm to train the expanded PKN against the data, a family of models with transfer functions chosen from a discrete number of possibilities is obtained.
For each of these models, we generate unprocessed models (, box F) by removing all cFL gates that are logically redundant with other cFL gates (e.g., in the gate “(B AND C) OR B activate D”, the AND gate is logically redundant with the “B activates D” gate). These gates are removed because they increase model complexity by using multiple logic gates to encode a relationship that can be specified by a simpler gate.
In our toy example, a family of twenty unprocessed models was obtained by training the expanded map () to in silico data (.ii.) using the discrete genetic algorithm. The unprocessed models from different optimization runs had similar topologies with the exception of the gate describing the relationship of MEK to its input nodes: TGFα and Akt (, brown and green dashed gates). Sixteen of the unprocessed models described the activation of MEK as depending only on TGFα (brown, dashed gate) whereas four described activation using the AND NOT gate (green, dashed gate).
Model reduction and refinement
In the model reduction and refinement stage (Steps 4–6), we determine which gates can be removed altogether as well as AND gates that can be replaced with one-input cFL gates without significantly affecting the MSE. We implemented the non-exhaustive heuristic search procedure described below on each unprocessed model and illustrate its application to our toy example ().
Reduction of trained cFL models.
In the fourth step, we remove or replace all gates for which the alteration does not increase the MSE of the unprocessed model over some threshold, which we term the ‘reduction threshold’. We use a range of reduction thresholds such that each unprocessed model results in several models, one for each reduction threshold used. Following this step, the resultant models are considered reduced models.
In the fifth step, we fix the model topology to that obtained during Step 4 and treat the transfer function parameters in each reduced
model (, Step 5) as continuous parameters rather than the discrete set of transfer function parameters required for use of the discrete genetic algorithm. We use a Sequential Quadratic Programming method (Text S1
) to refine the model parameters and further improve the fit of the models to the experimental data. The resulting models are termed reduced-refined
models, which have a range of MSEs depending on the reduction threshold used ().
In the sixth and final step, we specify a reduced-refined model to represent each unprocessed model (, Step 6). For each unprocessed model, we choose the reduced-refined model that has the fewest number of fitted transfer function parameters without increasing the MSE above a defined ‘selection threshold.’ The selection threshold is chosen by comparing the average number of parameters in the family of models to the average MSE of the models (). The net result is a set of reduced-refined-filtered models (hereafter referred to as filtered models, , Box G).
In our toy example, the filtered
models have identical topology and in no case does Akt inhibit MEK activation (). This topology is, in fact, the topology from which the in silico
data was derived. The ability of cFL to fit intermediate values made it possible to recover the correct model topology, whereas BL did not identify the correct model, and a gate linking TGFα to PI3K was consistently missing (, dashed arrow). Specifically, BL was unable to return the correct topology because nodes downstream of PI3K (Akt and JNK) were partially activated (0.32 and 0.19, respectively) under conditions of TGFα stimulation, and a BL model that included the TGFα to PI3K gate had a higher error (MSE
0.56) than a model that omitted the interaction (MSE
0.07). In contrast, the improved ability of cFL to model graded activities made it possible to recover the true network topology.
Adjusting the complexity of CellNOpt-cFL model training
While the expansion step (, step 2) captures the many possible combinations of AND and OR logic relationships between nodes, it also increases the complexity of the network, resulting in an increase in the size of the optimization problem. Depending on the biological network of interest, some or most of these AND gates might not be biologically relevant. For example, it is unlikely that six receptors must be active in order to activate another species, as would be the case for a six-input AND gate (instead, it is more likely to be a OR gate). A profusion of AND gates also makes the resultant networks difficult to interpret because most AND gates are in only a few models whereas the majority of models contain single-input and OR gates. Thus, the AND gates can effectively appear as system “noise”, interfering with visual assessment as well as computational analysis of the model topologies. Because of these potential complications, the expansion step can be limited to include only AND gates with a few inputs, depending on the complexity one would like to capture with the trained network models.
In the current paper, we have limited the search in the discrete genetic algorithm to a set of seven transfer functions. Use of more or fewer transfer functions is possible, but we found that seven transfer functions allowed us to represent a variety of input-output relationships without unduly increasing problem complexity to the point that the discrete genetic algorithm no longer consistently returned models that fit the data well (see Materials & Methods
Applying CellNOpt-cFL to protein signaling data from HepG2 cells
To test the ability of cFL modeling to analyze real biological data, we modeled a set of measurements describing the response of the HepG2 hepatocellular carcinoma cell line to various pro-survival, pro-death, or inflammatory cytokines in the presence or absence of specific small molecule kinase inhibitors. This dataset was used to construct a recent BL model 
. Here we ran an independent analysis using the cFL approach and compare the results to the BL previously reported. The dataset comprises measurement of phosphorylation states as markers of activation of 15 intracellular proteins before and 30 minutes after stimulation by one of six cytokines in the presence or absence of seven specific small molecule kinase inhibitors (, Figure S1
). The measurements were normalized to continuous values between zero and one using a routine implemented in the MATLAB toolbox DataRail 
, as previously described (
, see Text S1
Initial analysis of cFL models trained to HepG2 dataset.
The HepG2 dataset was trained to several related PKNs which are enumerated in and Figure S2
. These PKNs were derived, with various extensions, from the Ingenuity Systems database (www.ingenuity.com
) with manual addition of literature data about IRS1 that was obviously missing 
. The first PKN, termed PKN0 was identical the one used previously for BL modeling 
. In the course of our analysis, we found it necessary to search the literature for interactions missing in PKN0 but supported by the data, resulting in several PKNs (). Furthermore, we limited the manner in which the PKNs were expanded in two ways: (1) expansion into all possible two-input AND gates or (2) expansion into a two-input AND gate only when one input was inhibitory. In the second case, the expansion of inhibitory gates was necessary because, in logic terms, an inhibitory gate indicates that the output node is active when the input node is not active. In biological networks, this is true if the output node is constitutively active, which was not observed in the normalized HepG2 data. Thus, in order to accurately model the inhibitory effect, it had to occur in conjunction with activation by some other input node, which is captured by an AND gate. If a PKN was processed with both types of expansion, we include a superscript to differentiate between the two cases – i.e.
for the expansion of all gates and PKN1i
for the expansion of only the inhibitory case.
Prior knowledge networks trained to HepG2 dataset.
CellNOpt-cFL training of PKN0
PKN0 was expanded to include all possible two-input AND gates and trained to the HepG2 dataset with CellNOpt-cFL (Figure S2
). The 90 unprocessed
cFL models obtained after training showed that PKN0 exhibited a poor fit to IL1α-induced protein phosphorylation (Figure S3
), a result we had also observed with BL analysis 
, confirming that the poor fit of BL was due to errors in the topology of PKN0 and not the inability of Boolean logic to fit intermediate values.
An inspection of systematic model/data disparity (Figure S3
) immediately indicated that the models did not fit IL1α-induced phosphorylation of IRS1, MEK and several species known to be modulated by the MEK pathway. In PKN0, no paths between IL1α and MEK or IRS1 were present. Based on careful reading of the literature, we added two links to PKN0: a TRAF6 → MEK link 
, and an ERK → IRS1 link 
. These links had been inferred by the BL framework 
and were supported by further literature evidence. To add a link that provided a path between IL1α and MEK in the absence of BL inference results, for simplicity one should first consider links from species that IL1α is already known to activate. In this case, TRAF6 is the most upstream species which experimental evidence suggests can activate MEK 
. In the case of IRS1 signal activation, the specific phosphorylation site measured should be considered. Our data included measurements of phospho-S636/639, and S636 is a known phosphorylation site of ERK2 
A novel finding from CellNOpt-cFL analysis of the HepG2 data was that IL6 treatment led to phosphorylation of several downstream proteins. Similarly to the links just considered, PKN0 included no paths between IL6 stimulation and these downstream proteins, resulting in an inability to fit this pattern of phosphorylation. Importantly, however, BL analysis would not have recognized this partial activation due to its inability to fit intermediate values (as illustrated in our earlier toy example). Because IL6 was observed to partially activate Akt in the data and known mechanisms exist for this activation 
, we added a prospective IL6R → PI3K link to the PKN, thus providing an extended PKN (PKN1) that we use below for subsequent CellNOpt-cFL analysis.
CellNOpt-cFL training of PKN1
PKN1 was expanded to include all possible two-input AND gates (PKN1a) for a total of 170 discrete parameters corresponding to 105 logic gates. The resultant network was trained to the HepG2 data. Reduction of the PKN1a–derived models indicated that almost all AND gates could be removed or replaced by single-input gates. Since the AND gates appeared to add unnecessary complexity to the cFL models, we also expanded PKN1 to only include AND gates if an input node was inhibitory (PKN1i; ), resulting in only 60 discrete parameters corresponding to 56 logic gates. We then compared the PKN1a- and PKN1i-derived cFL models.
The comparison of these two PKN-derived model families revealed a clear tradeoff between model fit and complexity. The more complex PKN1a
-derived models were able to fit the data slightly better than the PKN1i
-derived models (average unprocessed
model MSE of 0.032±0.002 compared to 0.035±0.002, p
<0.001). However, the more complex PKN1a
-derived models contained many more parameters than the PKN1i
-derived models both before and after optimization (170 compared to 60 discrete parameters before optimization and an average of 72.8±4.9 compared to 66.6±3.9 continuous parameters after optimization (p
<0.001); Figure S4
). The simpler PKN1i
-derived models used fewer initial and final parameters to arrive at a fit to the data only 9% worse than PKN1a
-derived models. Since the 9% deviation is in the range of error in the normalized data (error estimated to be 10% by comparing similar stimulation conditions), we focused subsequent analysis on the simpler PKN1i
-derived models. For completeness, we include the results of PKN1a
-derived models as supplemental information (Figure S5
Statistical significance of cFL models trained to PKN1i
To determine the statistical significance of our results, we compared the family of 243 unprocessed
models with unprocessed
models obtained from either training PKN1i
to randomized data or training a randomized PKN1i
to the data (Table S1
). Data was randomized by pairwise exchange of all data values while network topologies were randomized either by generation of an entirely random topology or by random pairwise exchange of gate inputs, gate outputs, or nodes' inputs 
. When compared to the results of all types of randomization, models trained to the real data and PKN1 were highly significant (P-value <0.001, Table S1
), indicating that the family of trained cFL models fit the data better than expected by random chance.
To probe the dependence of the CellNOpt-cFL training process on the quality of the PKN used, we randomly added links to or removed links from the PKN and trained the resultant PKN to the data. As expected, the models derived from PKNs with links randomly removed had a poorer fit to data than those derived from the complete PKN1i
(, solid line). Conversely, when links were randomly added to the PKN, cFL-CellNOpt effectively removed the links (Figure S6
), resulting in models with similar goodness of fit as models derived from PKN1i
(, dashed line). We thus conclude that an incomplete PKN degrades the ability of CellNOpt-cFL to fit the data whereas models derived from a PKN with extraneous links retain this ability.
As an initial investigation of model predictive capacity and a check for over-fitting, we performed a ten-fold cross-validation by randomly dividing the HepG2 data into ten subsets and, for each subset, reserving one as a test set while training with the remaining nine data subsets. The similar fits of the training and test data provided evidence that the family of models obtained from this procedure were predictive, and the difference in test and training MSEs did not depend on selection threshold, a measure of model size, suggesting that the models were not over-fit ().
Analysis of this cross-validation result combined with a plot of average filtered model size and fit (MSE) as a function of selection threshold () suggested that a selection threshold in the range 1×10−3 – 1×10−2 would result in a family of models that contain slightly fewer number of parameters than lower thresholds (, dashed line) while retaining the ability to fit the data well (, solid line). We used a threshold of 5.0×10−3 for the remainder of our analysis unless otherwise noted.
Finally, we obtain a family of 243 filtered models for further analysis (). By taking note of which cFL gates are removed during the CellNOpt-cFL training and reduction processes, one can generate hypotheses regarding these gates. summarizes a set of biological hypotheses readily suggested by our cFL model topologies.
Structure of family of cFL models resulting from training PKN1i to HepG2 dataset.
Biological hypotheses about signaling network operation suggested by gates removed during CellNOpt-cFL analysis.
Validated biological hypothesis 1: Crosstalk from TGFα to the JNK pathway
Analysis of error between the family of cFL models and experimental data (Figure S7
) highlighted consistent error in TGFα-induced partial activation of c-Jun. Both PKN0 and PKN1 allowed for TGFα-induced activation of c-Jun by the JNK pathway via crosstalk from Ras or PI3K to MAP3K1. In the BL methodology, this crosstalk was removed due to the inability to fit partial activation, and no BL model allowed for activation of c-Jun after TGFα stimulation. However, we found that a subset of cFL models accounted for this c-Jun partial activation by including crosstalk between Ras or PI3K and MAP3K1. These models also partially activated JNK after TGFα stimulation, a feature that was inconsistent with the training data (Figure S8
). Thus, these models predict that JNK was actually phosphorylated under conditions of TGFα stimulation, but our measurements did not detect it.
To test this prediction directly, we undertook de novo
measurement of JNK and c-Jun phosphorylation following stimulation with different doses of TGFα (). These new data show that JNK does indeed become phosphorylated upon stimulation of HepG2 cells with TGFα. Thus, the cFL models containing crosstalk from Ras or PI3K to MAP3K1 were the correct models. Combined with , this analysis highlighted the partial activation of the JNK pathway after TGFα stimulation as a singular instance of crosstalk from a pro-growth ligand to an inflammatory pathway. In support of the significance of our finding here, we note that TGFα-induced JNK activation has been shown to be important for hepatic regeneration 
and stimulation of DNA synthesis 
in primary rat hepatocytes.
Validation of cFL crosstalk predictions.
Validated biological hypothesis 2: Mechanism of IL6-induced protein phosphorylation
As previously mentioned, PKN0 was unable to fit IL6-induced protein phosphorylation (a feature of the data unappreciated by the BL methodology). Because Akt was observed to be partially phosphorylated under these conditions and we found literature evidence for a prospective IL6R → PI3K link, we added the link to PKN1. However, the media-only condition also induced partial phosphorylation of Akt. Discovery of the partial activation of Akt in the media-only control led us to consider that perhaps the IL6-induced phosphorylation of Akt was simply an assay artifact. Thus, we inserted an Assay → PI3K link into the PKN. This “Assay” node represents cell stress arising from changing environmental conditions during the assay (media change, etc.); it is postulated to activate PI3K because only Akt is consistently active in the untreated control. Having accounted for the potential that IL6-induced partial phosphorylation of Akt was an artifact, we undertook a series of computational experiments to determine the mechanism of IL6-induced phosphorylation of downstream proteins.
Upon exposure to IL6, SHP2 has been reported to bind to gp130, a subunit of the IL6 receptor complex. SHP2 is then phosphorylated in a JAK1-dependent manner. This phosphorylation can lead to PI3K/Akt pathway activation through interactions with Gab-1 or IRS1 or Ras/MEK/ERK pathway activation through Grb2 or Gab1 
. Thus, our computational experiments were designed to infer which pathway (PI3K/Akt or Ras/MEK/ERK) was mediating the IL6-induced protein phosphorylation. Four families of 150 filtered
models were examined, all of which were obtained after training a new PKN to the normalized HepG2 dataset (, PKN2A – PKN2D). The inability of PKN2A-derived cFL models with only the Assay → PI3K link to fit well the IL6-induced protein phosphorylation data suggested that some other link was necessary to fit this data. In our trained networks, the IL6R → PI3K link was present in only a fraction of the relevant trained models (PKN2B and PKN2C), but the IL6R → Ras link was present in more than 90% of relevant trained models (PKN2C and PKN2D). Additionally, models with IL6R → Ras links were better able to fit the IL6-induced protein phosphorylation. Consequently, our cFL results supported the hypothesis that IL6R activates downstream proteins through the Ras/Raf pathway. This hypothesis is supported by an independent dataset 
, where the IL6-induced protein phosphorylation response was more robust than in the training data (Figures S1
). Inhibition of MEK either alone or in combination with other inhibitors resulted in ablation of downstream protein activation whereas inhibition of PI3K did not (). Thus, we infer that IL6-induced protein phosphorylation was not an assay artifact and was instead mediated by the Ras/Raf pathway.
Results of cFL training of various prior knowledge networks for the investigation of IL6 crosstalk.
Predicting node-to-node transfer functions
CFL relates nodes in a network with transfer functions that describe quantitative input-output relationships between protein species represented as network nodes. To investigate the ability of the cFL models to predict these transfer functions, we simulated the PKN1i-derived, filtered cFL models to determine the activation state of a specified node under many theoretical combinations of its input nodes. We then plotted the model predictions of quantitative input-output relationships. As one instance, shows the predicted average and standard deviation of the quantitative values of CREB phosphorylation as a function of the activation of upstream nodes, p38 and MEK1/2. The resulting plots indicated that we were able to predict the activation response of CREB to the entire range of p38 and MEK1/2 although training set measurements were limited to a few values of these nodes (, black circles).
Transfer functions predicted by trained cFL models.
We tested this prediction using a set of data with combinations of ligands and inhibitors not present in the training data (
, Figure S9
). Roughly 20% of the test conditions were also present in the training data set, allowing us to control for differences between both data sets. When we compared this dataset to the predicted transfer functions, we observed that most of the data fell within one standard deviation of the predicted value (, green diamonds) with exception of overestimation under conditions of TGFα stimulation. This overestimation is expected, as a comparison of common conditions between the training and test dataset indicated that the normalized experimental values of CREB in the validation dataset were 38±4% lower than that in the training set.
This result demonstrates the ability of the trained cFL models to predict the quantitative relationship between nodes in the network. We also found that the family of cFL models was able to fit the phospho-protein signaling response in the validation dataset well, which we demonstrate as supplementary information (Figure S9
Predictive capability of a cFL model family
We performed a series of nineteen cross-validation experiments to further investigate the ability of our methodology to predict the signaling response under conditions that were not represented in the training data. For each experiment, we used training data from which we had removed the phosphorylation data of a specific protein signal, s
, under a single ligand stimulation condition and all inhibitor treatments. Nineteen signal/stimulation combinations were chosen to be test sets according to two criteria: (1) s
is at least partially activated under the stimulation condition of interest and (2) s
is at least partially activated under some other stimulation condition (Table S2
). These criteria ensured that the remaining training data contained some information regarding the activation of s
but it did not contain information regarding the activation of s
under the stimulation condition of interest. This procedure is a more stringent test for predictive capability than a random cross-validation procedure because training sets from which random data is removed might retain other data with the same information as the removed data (e.g.
, based on the network topology, Akt phosphorylation in the absence of MEK inhibition is the same as Akt phosphorylation with MEK inhibition, so removing only one of these data points is not a stringent test of predictive capacity).
We examined the ability of models trained on reduced training sets (n>45 for each case) to predict phosphorylation of the test protein signals. Because we used each individual in the family of models to predict the test signal, we could determine if the models were constrained in their predictions by examining the coefficient of variance (CV; standard deviation divided by mean) of the prediction. If the CV was high, the models were not constrained to a specific prediction (i.e. the prediction was imprecise), and the average prediction should be discounted. Thus, for these cross-validation results, we compared the precision (CV) and accuracy (MSE) of the models' predictions, where precise and accurate predictions exhibited both a low CV and low MSE ().
Accuracy vs. precision of cross-validation experiments.
We found that the families of models trained on these reduced training sets were able to precisely predict phosphorylation of the test protein signals in twelve of the nineteen cases (, green field). In six of the test sets, the models did not agree, although their average prediction was reasonably accurate (, yellow field). We observed no test sets for which the training sets agreed about an inaccurate prediction (, orange field). In one case (prediction of Iκb signaling under TNFα stimulation), the predicted phosphorylation state was highly inaccurate (MSE >0.20). However, this prediction was also very imprecise (CV >0.25), indicating that the average prediction was unreliable (, blue field). Thus, by taking the precision of the models' predictions into account, we were able to discredit an inaccurate prediction. This result underscores the importance of considering consensus among the family of models rather than examining the results of only one cFL model.
Using cFL models to relate phospho-protein signaling to cell phenotypic response
The ability to quantitatively model protein signal activation with cFL offers the prospect of predicting phenotypic response upon exposure to stimuli and inhibitors. To investigate the ability of cFL to model phenotypic data, we turned to data describing cytokine release three hours after stimulation under the same conditions as the phosphorylation data 
. As a first approach, we linked the output of our family of cFL models to a partial least squares regression model 
obtained by regressing normalized data of release of five cytokines (IL1β, IL4, G-CSF, IFNγ, and SDF1α) to the normalized protein phosphorylation measurements (see Text S1
The cFL models linked to a PLSR model were able to model phenotypic response with an accuracy of R2
0.79, near that of the PLSR model (R2
0.81; see Figures S10
). However, we found that the correlation indicated by regression coefficients did not lead to easily interpretable insights about phenotype because proteins in the same pathway were also highly correlated with each other.
To obtain a more interpretable model, we utilized a second approach where we included nodes specifying cytokine release in the PKN and linked them to a few protein signaling nodes. These nodes were chosen based on principle component analysis: if protein signals in a pathway clustered together in principle component space, the signal most downstream in the pathway was linked to cytokine release. Based on this analysis, the following protein signaling nodes were linked to each cytokine release node: MEK1/2, CREB, GSK3, c-Jun, Hsp27, Iκb, and STAT3 (, PKN3). We then trained a family of cFL models to the normalized dataset comprised of cytokine release at three hours and protein signaling at thirty minutes.
The resultant models were able to fit the cytokine release data reasonably well (R2
0.78 for the average predicted by a subset of best-fitting models, Figure S11
). Furthermore, the low frequency of several gates in the resultant family of cFL models (Figure S12
, Table S3
) indicated that, although the promoters of several of the modeled cytokines contained binding sites of transcription factors are known to be modulated by the MEK1/2, GSK3, and CREB pathways (Table S4
), activation of these nodes did not predict cytokine release. Thus, we altered our previous PKN by removing the links between these protein signaling and cytokine release nodes and trained it to the data. The resultant family of cFL models () indicated that STAT3 activation explained cytokine release after IL6 stimulation and other signals (Iκb, c-Jun, and Hsp27) explained cytokine release three hours after TNFα or IL1α stimulation.
Trained cFL models linking ligand cues, phospho-protein signals, and cytokine release phenotypic responses.