Computational Verification of Large Logical Models—Application to the Prediction of T Cell Response to Checkpoint Inhibitors

At the crossroad between biology and mathematical modeling, computational systems biology can contribute to a mechanistic understanding of high-level biological phenomenon. But as knowledge accumulates, the size and complexity of mathematical models increase, calling for the development of efficient dynamical analysis methods. Here, we propose the use of two approaches for the development and analysis of complex cellular network models. A first approach, called "model verification" and inspired by unitary testing in software development, enables the formalization and automated verification of validation criteria for whole models or selected sub-parts. When combined with efficient analysis methods, this approach is suitable for continuous testing, thereby greatly facilitating model development. A second approach, called "value propagation," enables efficient analytical computation of the impact of specific environmental or genetic conditions on the dynamical behavior of some models. We apply these two approaches to the delineation and the analysis of a comprehensive model for T cell activation, taking into account CTLA4 and PD-1 checkpoint inhibitory pathways. While model verification greatly eases the delineation of logical rules complying with a set of dynamical specifications, propagation provides interesting insights into the different potential of CTLA4 and PD-1 immunotherapies. Both methods are implemented and made available in the all-inclusive CoLoMoTo Docker image, while the different steps of the model analysis are fully reported in two companion interactive jupyter notebooks, thereby ensuring the reproduction of our results.

Introduction

Recent technical developments have allowed scientists to study immunology and health-related issues from a variety of angles. For many diseases, especially for cancer, the current trend consists in aggregating data coming from different sources to gain a global view of cell, tissue, or organ dysfunction. Over the last decades, diverse mathematical frameworks have been proposed to seize a multiplicity of biological questions (Le Novère, 2015), including in immunology (Kaufman et al., 1985, 1999; Eftimie et al., 2016; Chakraborty, 2017). However, the increasing complexity of biological questions implies the development of more sophisticated models, which in turn bring serious computational challenges.

Among the mathematical approaches proposed for the modeling of cellular networks, the logical modeling framework is increasingly used. In particular, it has been successfully applied to immunology and cancer, leading to the creation of models encompassing dozens of components, some including many inputs components (Grieco et al., 2013; Abou-Jaoudé et al., 2014; Flobak et al., 2015; Oyeyemi et al., 2015). However, the large size of recent models hinders the complete exploration of their dynamical behavior through simulation, especially in non-deterministic settings.

To address these difficulties, we define and apply a model verification approach to systematically verify whether a model complies with a list of known properties. These properties are defined as model specifications, either at a local (i.e., for sub-models) or at a global level. This automated verification procedure fosters confidence during the development of a complex dynamical model and paves the way to the development of models with hundreds of nodes.

We further outline and apply a value propagation method, which enables the assessment of the impact of environmental or genetic constraints on the dynamical behavior of complex cellular networks.

These two complementary approaches can be applied to the development and analysis of large dynamical models, as illustrated in Figure 1. Noteworthy, they have been implemented in a multi-platform Docker image combining various complementary logical modeling and analysis tools (Naldi et al., 2018b). We further illustrate the power of these methods through the analysis of an original model. The different steps of analysis are fully reported in two companion interactive jupyter notebooks, available with the model on the GINsim website (http://ginsim.org/model/tcell-checkpoint-inhibitors-tcla4-pd1), thereby ensuring their reproducibility.

FIGURE 1

Figure 1. Description of the proposed workflow for the development and analysis of dynamical logical models. The novel methods described in this article are emphasized with blue fonts. Starting with the delineation of a molecular map integrating the available scientific knowledge, we derive a regulatory graph and logical rules to generate a logical model, and induce dynamical specifications serving as test cases to verify the model. Moreover, when the available knowledge is specific to a smaller part of the regulatory graph, a sub-model us extracted to perform local tests. We further implemented an analysis and visualization method, called Value propagation, to assess the impact of environmental and genetic perturbations. Figure 3 zooms into this part of the workflow and describes it in more details. The use of model verification, sub-model extraction and value propagation is illustrated in two reproducible and editable Jupyter notebooks, taking advantage of the CoLoMoTo Interactive notebook framework (Naldi et al., 2018b). This framework is available inside the CoLoMoTo Docker image together with packaged libraries for the analysis of dynamical logical models of biological networks.

Model Verification

A Software Engineering Framework for Logical Model Building

One of the main features determining the interest of a model is its ability to accurately recapitulate salient biological knowledge. More precisely, this knowledge can be used in two complementary ways during the model building process. On the one hand, it is used to define the model architecture, specifying which biological entities need to be included and which interactions between these entities need to be encoded. On the other hand, biological knowledge entails dynamical properties that must be achieved by the resulting model, whether transitory or asymptotic, to account for biological observations. These properties induce satisfaction criteria and must be clearly specified for rigorous model assessment or comparison with other models. Failures to reproduce such properties need to be carefully documented, thereby providing a basis for further model improvement.

In the domain of logical modeling applied to cellular networks, various formal methods have already been proposed to verify dynamical properties. For example stable states (or fixed points, characterized by all components being steady at the same time) tentatively correspond to asymptotic properties that can used to assess the reproduction of known persistent biological behavior. More complex asymptotic behaviors include cyclic attractors, which can be approximate by the computation of so called trap spaces. Also called stable motifs, trap spaces are hypercubes in the state space such that all successors of all states in the hypercube also belong to it (for synchronous and asynchronous updatings, or any other updating). These hypercubes then provide an approximation of complex attractors. Trap spaces and stable states can be defined as results of a constraint solving system, enabling their efficient computation (Klarner et al., 2018). Their reachability however must be assessed separately, often using model checking or stochastic simulations, which requires longer computations.

Model checking techniques have been successfully applied to specify and verify temporal constraints on a model behavior (Monteiro and Chaouiya, 2012; Miskov-Zivanov et al., 2016; Traynard et al., 2016; Wang et al., 2016).

In any case, whatever the formalism chosen, the building of a complex dynamical model is intrinsically iterative, as its establishment is usually incremental and requires continuous testing and adjustment with reference to a growing body of biological knowledge.

In the field of software engineering, the similar need to repeatedly assess criteria of success or failure of a software program led to the development of powerful software verification techniques, and in particular to software testing (Myers, 1979), which main goal is to assess whether a software meets a series of well-defined requirements. More importantly, such assessments must be repeated as soon as a new piece of code or specification is added. Software testing aims to check whether newly introduced modification might break any of the previous performances. In particular, software verification includes the notion of unit testing, where suites of tests describe the expected behavior associated with individual units composing a program. This idea can be transposed from computer science to model building and has been successfully applied in the context of other modeling frameworks (Hoops et al., 2006; Lopez et al., 2013; Sarma et al., 2016; Boutillier et al., 2018), but not yet to logical modeling.

Here, we transpose the unit testing approach to integrate a comprehensive series of verifiable criteria, from the early stages of model conception, in order to automate the dynamical evaluation of logical models. The core idea is to split the biological knowledge on which a model is based into individual verifiable criteria that can be formalized as specifications (Figures 1, 2). In this respect, individual units of knowledge, derived from the scientific literature or biological experiments, must be formulated into stable or dynamical properties. Each specification, coupling a property with an expected value, can serve as a basis to define a test case for a model. Testing such a specification amounts to compute an “observed” value based on the model and compare this value to the expected one.

FIGURE 2

Figure 2. Sub-model extraction for local model verification. When the available knowledge is fragmentary and covers the behavior of only a subset of components, verification becomes difficult at the global scale of the model. Based on this partial information, a series of specifications can still be defined for a sub-model which, after extraction using bioLQM's submodel() function, can then be rigorously tested.

In practice, the CoLoMoTo notebook environment (Naldi et al., 2018b) provides a Python API for several software tools, enabling the definition of a wide range of dynamical analyses for the computation of observed values. Individual test cases can be assembled into a library, also called testing suite. Existing tools and packages enabling software testing can then be applied to automatically assess whether a model satisfies (or not) a series of specifications. In this study, we used the python package “unittest,” taking advantage of its seamless integration into the CoLoMoTo interactive notebook. This unit-testing package is integrated by default into the recent versions of the Python standard library (http://python.org).

Local Verification of Sub-Models Can Cope With Sparsity of Biological Knowledge

Biological knowledge reported in the scientific literature is often insufficient to evaluate a comprehensive model, which may encompass hundreds of nodes. In particular, observations regarding component activity often relate to only a limited subset of nodes of the model. This greatly complicates the definition of specifications for the whole model.

Given a comprehensive model and a set of components of interest, one can extract a sub-model containing these core components, along with their associated logical rules. Components appearing in these logical rules but not in the selected set are considered as external inputs of the sub-model (Figure 2). This functionality has been implemented in the “submodel” function of the Java bioLQM library (Naldi, 2018) according to the following procedure.

Let M = (V, f) be a model, where V is the set of components, and f the update function. For each c in V, f_c is the logical function associated to the component c and R(c) is the set of its regulators (i.e., components that intervene in the logical rule). Given a list of selected components C ⊂V:

1. S = ∅

2. for each component c ∈ C: S = S + {c} + R(c)

3. create the sub-model M′ = (S, f′) such that for each component c in S: $f_c^{\prime }=\{\begin{array}{cc}{f}_c& \mathrm{\text{if }}R(c)\subset S\\ c& \mathrm{\text{otherwise}}\end{array}$

As shown in the application below, the delineation of such sub-models can greatly facilitate the definition and verification of local specifications.

Value Propagation Enables the Evaluation of the Impact of a Given Cellular Environment on Model Dynamics

The core idea of value propagation is presented in Figure 3. Given a set of logical rules and a cellular context, an iterative algorithm enables the computation of the dynamical consequences of the cellular context on all the components of the model.

FIGURE 3

Figure 3. The principle of logical value propagation analysis is illustrated with a simple example involving two core nodes, D and E, and four input nodes, A, B, C, and F. The value 1 is assigned to the node C and then propagated through the model. The assignment C:1 implies the evaluation of D to 1. Consequently, the function assigned to node E becomes “not F.” In other words, assigning the value 1 to node C activates node D independently from the value of its other inputs, while node E becomes completely dependent on the value of node F.

First, the cellular context is formalized by assigning constant values to some components of the model. Next, we apply a recent model reduction technique reported by Saadatpour et al. (2013). Briefly, for each constant node, the corresponding value is inserted into the logical rule associated with each of its target nodes. Each logical rule is then simplified using Boolean algebra. If the rule simplifies to a constant, this fixed value is further propagated into the logical rules of downstream nodes. This process is iterated until no further propagation or simplification can be made on the logical rules of the model. In contrast with the approach of Saadatpour et al., which aims at producing a reduced model, we focus principally on the outcome of the propagation of fixed values.

The result of value propagation can be very informative by itself. Indeed, the resulting stabilized values provide insights into the impact of a given (single or multiple) perturbation on the model, revealing which elements are consequently constrained to become activated or inactivated, vs. which elements keep some degree of freedom. Furthermore, this method greatly eases the comparison of the impacts of different biological contexts on network dynamics by performing a differential analysis of the corresponding lists and target values of fixed components. This method has multiple advantages when applied to complex networks, as it can be used efficiently on models with large numbers of components. It further simplifies the computation of attractors (stable states or even simple or complex cycles). Interestingly, Saadatpour and collaborators showed that this method conserves the stable states and complex attractors under the fully asynchronous updating assumption (Saadatpour et al., 2013).

This method was extended to multilevel models and implemented into the Java bioLQM library (Naldi, 2018). In this implementation, the fixed components are conserved during value propagation, enabling a direct comparison of the propagated effects of alternative perturbations.

The power of this approach is demonstrated on a concrete example in the following section.

Application: Assessing the Effect of Checkpoint Blockade Therapies on T Cell Activation

Biological Background

Over the last decades, immunotherapies have been the subject of intense studies and led to great advances in the field of cancer treatment. Through the years, it has then been recognized that T cells often display a reduced ability to eliminate cancer cells, and that expression of co-inhibitory receptors at their surface accounts for this compromised function. Receptors like Cytotoxic T-lymphocyte protein 4 (CTLA4, also known as CD152) (Walunas et al., 1994; Leach et al., 1996) and Programmed cell death protein 1 (PD-1, also known as PDCD1 or CD279) (Ishida et al., 1992) have been particularly studied in that context. Antibodies blocking the pathways downstream of these co-inhibitors (checkpoint blockade therapies) have become standard treatment for metastatic melanoma (Robert et al., 2011; Simpson et al., 2013) and other cancers (Ribas and Wolchok, 2018), including non-small cell lung cancer, renal cell carcinoma, Hodgkin's lymphoma, Merkel cell carcinoma and many others. The successes of these studies led to an increasing interest in T cell co-inhibitory receptors.

Nevertheless, a clear understanding of the mechanisms at work inside T cells remains elusive. Therapies targeting CTLA4 or PD-1 show different immune adverse effects (June et al., 2017), while the corresponding intra-cellular mechanisms remain to be clarified. Moreover, a rationale for the educated development of new immunotherapies focusing on other receptors or combinations of receptors is clearly needed. Co-inhibitory receptors are legions at the surface of T cells (Brownlie and Zamoyska, 2013) and biology of T cell activation or tolerance involves activation or repression of highly interconnected and complicated pathways (Baumeister et al., 2016).

Given the central role of T cells in many medical contexts, several mathematical frameworks have been applied to model T cell activation. Recent examples include rule-based approaches (Chylek et al., 2014), ordinary differential equations (Perley et al., 2014), and logical models (Oyeyemi et al., 2015; Rodriguez-Jeorge et al., 2019; Sanchez-Villanueva et al., 2019), considering different biomedical contexts as diverse as HIV infection or neonate vaccination. To our knowledge, none of them specifically focused on the impact of co-inhibitory receptors on T cell activation or tolerance.

In this study, we applied the logical framework to integrate current data on CTLA4 and PD-1 pathways and assess their impact on T cell activation. Our goal was triple. First, we wanted to create a comprehensive model building upon extensive knowledge encoded into a molecular map (see next section). Second, using model verification and a specific unit test suite, we aimed to firmly anchor the model at both the global and local scale into the collected biological knowledge. Third, using value propagation, we aimed to provide a tool for the comparative analysis of intra-cellular consequences when targeting CTLA4 vs. PD-1 T-cell co-receptors.

Comprehensive Molecular Mapping of T Cell Activation Network

Prior to mathematical modeling, knowledge about biological entities involved in T Cell activation was collected from available pathway databases, including Reactome (Fabregat et al., 2016), PantherDB (Mi et al., 2013), ACSN (Kuperstein et al., 2015), and WikiPathways (Slenter et al., 2018). Moreover, the scientific literature indexed in the PubMed database was further explored and carefully curated. Using the software CellDesigner (version 4.3.1) (Funahashi et al., 2008), this knowledge was encoded in a molecular map describing reactions between biological entities (either proteins, RNAs, genes, complexes, or metabolites). Each biological entity included in the map was annotated with a series of standard identifiers, including UniProtKB accession number, recommended and alternative names, gene name and synonyms, and cross-references to unique HGNC identifiers and approved symbols. The annotations also reference relevant scientific articles, including PubMed identifier, first and last authors, year of publication, and a list of observations extracted from these publications.

Our T cell activation map currently encompasses 726 biological entities, in different states (active/inactive, with or without post-translational modifications), and 539 reactions involving these entities (Supplementary Figure 1 and Supplementary File 1). Globally, the map currently integrates information from 123 scientific articles, which are cited in the annotations of the entities and reactions of the map.

Logical Modeling of T Cell Activation

Using the logical modeling software GINsim (version 3.0.0b) (Naldi et al., 2018a), we then manually derived a regulatory graph encompassing 216 nodes and 451 arcs (Figure 4) from the content of the molecular map. One by one, biological entities represented in the molecular map were re-created as components of the logical model. In most of the cases, the representation of entities having different states was further compressed into a single component summarizing their activity in the TCR signaling cascade. Furthermore, to obtain a dynamical logical model, a specific logical rule must be assigned to each node. In many cases, this can be achieved rather easily based on published data. For more complex situations, a default generic logical rule was initially considered, where all activators are needed for the activation of a component (using the AND operator) and where only one inhibitor is sufficient to repress it (using the OR and NOT operators), which served as a basis for further rule refinement. In some cases, however, in particular when a component is the target of various regulatory interactions or when metabolites are involved, finding direct support for a specific rule may be tricky or impossible. Hence, the delineation of consistent logical rules for a complex model is often the result of an iterative process, starting with generic rules and progressively correcting them based on the results of various analyses.

FIGURE 4

Figure 4. Regulatory graph of the T cell activation model. The global layout is similar to the molecular map (cf. Supplementary Figure 1 and Supplementary File 1), with ligands/receptors and proximal signaling at the top, and the nucleus-related events at the bottom of the graph. In between, the model encompasses interconnected pathways and signaling cascades related to cytoskeleton remodeling, the MAPK network, calcium fluxes, metabolic shifts, and NF-κB, to name a few. Boolean components are denoted by ellipsoids whereas rectangles denote ternary components. Green arcs denote activation events, red blunt arcs denotes inhibitions, while blue arcs denote dual regulations. The gray arcs represent interactions created during the translation of the molecular map into the regulatory graph, but that are not yet integrated at the dynamical level (i.e., not taken into account in the logical rule).

Hereafter, we demonstrate how we can take advantage of the methods presented in the previous sections to ease rule refinement by model verification. We first defined a series of properties expected for the model (see examples in Table 1). Next, stable states and/or trap spaces were computed and automatically compared with these properties (cf. first Jupyter notebook provided on the model web page at http://ginsim.org/model/tcell-checkpoint-inhibitors-tcla4-pd1). After some iterative runs of the notebook, manual refinements lead us to a set of rules complying with all the tests.

TABLE 1

Table 1. Global specifications used to assess the T cell activation model and example of local specifications for the calcium signaling module (Cf. Figure 2).

For example, the Endoplasmic Reticulum (ER) serves as a reservoir for calcium ions. This reservoir can be emptied through activation of the Inositol 1,4,5-trisphosphate receptor (IP3R1). When empty, this reservoir can be filled through activation of the Sarcoplasmic/endoplasmic reticulum calcium ATPase 2 (SERCA) pumps. A default logical rule for a node representing the presence of this Calcium quantity (Calcium_ER) is then “SERCA AND NOT IP3R1.” To check the behavior of the corresponding logical sub-model, we defined a test checking whether whenever Calcium_ER was evaluated to TRUE, SERCA was evaluated to FALSE (see test “test_calc_tp_rest_ER1_SERCA0”). However, consecutive model verification failed, allowing us to notice that the default rule implied that SERCA should be always TRUE for Calcium_ER to be TRUE. The rule was then corrected to take into account the fact that Calcium_ER should stay TRUE whenever it would reach this value in absence of IP3R1.

In the first Jupyter notebook provided as Supplementary Material, we include all the code enabling the verification of our final model, which encompasses 36 unit tests split in four test suites. On a MacBook Pro using macOS 10.13 High Sierra, with a 2.3 GHz Intel Core i7 and 16GB 1600 MHz DDR3, all the tests were run in 87s.

The four test suites cover the most complex parts of the model, some of them particularly difficult to define. These suites use sub-models, whose delineation was guided by known pathways and practical knowledge gained by the modeler during the assembly of the molecular map. The Calcium module test suite covers a sub-model related to the fluxes of Calcium ions between different cellular compartments, namely the endoplasmic reticulum, the cytoplasm, and the extracellular region. The LCK module test suite is centered on the Tyrosine-protein kinase Lck (LCK). This kinase is known to have multiple sites of phosphorylation, whose collective status determines the tridimensional conformation and thus the activity of the enzyme (Ventimiglia and Alonso, 2013). The Cytoskeleton module test suite covers the cytoskeleton remodeling events occurring during T cell activation, and has strong connections with the Calcium sub-model. Finally, the Anergy/activation/differentiation module covers a less documented module encompassing the nucleus compartment and gene transcription.

Comparison of the Impacts of CTLA4 and PD-1 Co-Inhibitory Receptors Through Value Percolation

Based on the model described in the preceding section, a comparative propagation analysis was performed to visualize the respective effects of CTLA4 and PD-1 receptor activation on model dynamics. Figure 5 displays the value propagation for each condition on a single regulatory graph, using a color code to distinguish the different situations (component inhibition/activation in one or both conditions). The value propagation for the two conditions are further shown separately in the second companion notebook (available at http://ginsim.org/model/tcell-checkpoint-inhibitors-tcla4-pd1). This analysis reveals that the activation of the CTLA4 receptor impacts most pathways of the model, impeding in particular the remodeling of the cytoskeleton and the metabolic switch associated with bona fide T cell activation. In contrast, the activation of the PD-1 receptor leads to more limited effects, predominantly freezing the components of the NF-κB pathway.

FIGURE 5

Figure 5. Visualization of the results of the propagation analyses for CTLA4 vs. PD-1 activation. Gray nodes correspond to inputs. Nodes in yellow are frozen OFF upon any of CTLA4 or PD-1 (PDCD1) activation. Nodes in orange are frozen ON (i.e., with level 1 or 2) for each of these conditions. Nodes in light blue are frozen OFF only for CTLA4 activation (i.e., they remain free upon PD-1 activation). Nodes in dark blue are frozen ON (i.e., with level 1 or 2) only upon CTLA4 activation (i.e., they remain free upon PD-1 activation). Upon PD-1 activation, the corresponding node (PDCD1) is the only one that gets specifically frozen (ON, shown in dark green). Nodes in white remain free for both conditions.

A more refined comparative analysis of value propagation from these two receptor activations entails the observation that the set of nodes frozen by the propagation of PD-1 activation is completely included inside the set of nodes frozen by the propagation of CTLA4 activation (see Table 2). Furthermore, the values of the components frozen in both propagation studies are the same. Interestingly, a set of nodes related to calcium influx from and to the endoplasmic reticulum remain unfixed by any of the propagation analyses. This could be an artfact of the positive feedback loops added on the nodes representing the Calcium ion levels in different compartments and would need to be further investigated. A more detailed biological interpretation of these results is proposed in the following section.

TABLE 2

Table 2. Quantification of the model nodes impacted by the propagation of CTLA4 or PD-1 persistent activation.

Conclusions and Prospects

In this study, we have implemented and applied two complementary methods enabling a specification-oriented model building approach, thereby easing the delineation and analysis of highly complex logical models. In this respect, the building of a knowledge base, e.g., in terms of a molecular map, is an important first step. In the molecular map (provided as the Supplementary File 1), we have integrated the most relevant biological references available on T cell activation and inhibition pathways.

This map is clearly due to evolve, in particular thanks to the generation and analysis of novel high-throughput data (see e.g., the recent extensive analysis of the TCR signalosome by Voisinne et al., 2019). But any modification needs to be manually propagated to the dynamical model. To date, methods to derive proper dynamical models from such molecular maps are still in their infancy. In the particular case of the Boolean framework, only one automated approach has been recently proposed (Aghamiri et al., 2020). However, a limitation of this approach is the generation of generic logical rules based on static knowledge. Hence, the methods presented here could be used to advantageously refine these rules, taking into account additional biological knowledge about the behavior of the system under study.

We used the information gathered in our T cell activation map to build a dynamical logical model encompassing over 200 components and 450 interactions. For such a complex model, defining the logical rules in concordance with biological knowledge is a difficult and error-prone process, usually involving iterative trial simulations, where failures are identified to suggest potential improvements. Hence, listing comprehensive and consistent model specifications is a crucial step for model construction. These specifications can be revised as the modeler deepens his understanding of the biological processes under study. Noteworthy, such systematic testing procures a sense of confidence during the development process.

In the unit tests developed for our model, the definition of sub-models was guided by biological knowledge and pathway definitions, while relying partly on the modeler intuition. This step could be improved by community analyses of the regulatory graph to improve their definition.

Model checking techniques have been previously applied to assess model behavior through systematic cycles of model refinements (see e.g., Traynard et al. (2016) and reference therein). Model verification, as defined here, is a generalization of this approach, as it can rely on any available analysis as long as its result can be compared to an expected outcome. In our hands, in the course of model building, the unit testing approach, strongly anchored to available knowledge, proved to be very efficient to assess and improve model consistency with respect to a list of biological specifications, without the need of time-consuming and costly simulations. Implemented in the CoLoMoTo Interactive notebook framework (Naldi et al., 2018b), this approach enabled us to define a model recapitulating the most salient properties observed in response to T cell activation, including quiescence, anergy, and differentiation.

The use of model checking techniques could be further extended to assess the sensitivity of model behavior to the choice of specific logical rules. Such extension is hindered by the exponential increase of the number of possible logical rule, as the number of regulators increases. We would thus need a rationale to explore the space of logical rules. A first step in this direction can be found in Abou-Jaoudé and Monteiro (2019).

The approach presented here could also be improved by taking into account and tracking uncertainty during model conception (Thobe et al., 2018), or yet by taking advantage of computational repairing methods (Gebser et al., 2010) to identify more precisely remaining inconsistencies with biological data. Furthermore, other software engineering techniques, such as code coverage, could be borrowed to further improve model building and verification. Code coverage computes how much of a program's code is covered by unit tests. Similarly, one could design a method computing the fraction of the components of a model that is effectively covered by specifications.

Value propagation analysis of our large and complex regulatory graph proved to be biologically insightful. Indeed, this straightforward approach enabled us to clearly contrast the respective impacts of CTLA4 and PD-1 on T cell activation in our model, providing some rationale for their differential effects in current therapeutic studies. Indeed, anti-CTLA4 immunotherapies are known for their strong adverse effects related to autoimmunity and immunotoxicity (June et al., 2017). Anti-CTLA4 immunotherapies are currently combined with anti-PD-1 immunotherapy, known for its milder impact on the immune system.

Interestingly, the state of the node representing the Interleukin 2 (IL2) cytokine activation illustrates the differences of action of these receptors. Activation of the IL2 gene depends mainly on the activation of three transcription factors: the Nuclear Factor of Activated T cells (NFAT), the AP1 complex, and the Nuclear factor NF-kappa-B (NF-κB) (Smith-Garvin et al., 2009). When NFAT and AP1 are both active, they form a complex and together bind a regulatory region of the IL2 gene. In absence of AP1, NFAT induces a different program leading to cellular anergy (Macian, 2005; Smith-Garvin et al., 2009): activation of Diacylglycerol Kinase (DGK) prevents DAG-mediated activation of RasGRP1, which regulates the threshold for T cell activation (Roose et al., 2007; Das et al., 2009).

Our comparative propagation analysis reveals that while the activation of the CTLA4 receptor leads to a general inactivation of the three transcription factors regulating IL2 production, activation of the PD-1 receptor leads only to the inactivation of NF-κB and FOS (a member of the AP1 complex), thereby preventing the formation of the NFAT/AP1 complex, but enabling the activation of DGK. This observation is consistent with the proposal to target DGK isoforms as a complement of checkpoint immunotherapy (Riese et al., 2016; Jung et al., 2018).

As a next step, new co-inhibitory receptors recently under study, such as the Hepatitis A virus cellular receptor 2 (also known as TIM3) or the Lymphocyte activation gene 3 protein (LAG-3) (Anderson et al., 2016), could be easily added to the model described here, provided sufficient information could be gathered regarding their interacting partners. Applying propagation analysis in this context would be greatly insightful for future therapy developments.

Data Availability Statement

All datasets presented in this study are included in the article/Supplementary Material.

Author Contributions

CH developed the T cell signaling molecular map and model under the supervision of MT-C and DT. CH and AN implemented the computational methods and applied them to the T cell model under the supervision of AN and DT. All co-authors contributed to the redaction of the manuscript and endorse its content. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by a grant from the French Plan Cancer (project SYSTAIM, 2015–2019).

Conflict of Interest

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.

Acknowledgments

We warmly thank Bernard Malissen, Romain Roncagalli, and Guillaume Voisinne for their thoughtful guidance during the construction of the T cell signaling model.

Supplementary Material

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

Supplementary Figure 1. Molecular map in good resolution or vectorial format Caption: This molecular map describes biological entities and reactions implicated in the process leading to activation of CD4+ T cells in humans/mouse. The cytoplasmic membrane and attached receptor proteins are placed at the top, while the cell nucleus is located at the bottom. Reactions represent current knowledge of various pathways related to cytoskeletal remodeling, calcium fluxes, metabolism, cell cycle or IL2 production, and are encoded using the CellDesigner software (Funahashi et al., 2008) (the CellDesigner file is provided as Supplementary File 1).

Supplementary File 1. Molecular map in CellDesigner/XML format (Funahashi et al., 2008).

References

Abou-Jaoudé, W., and Monteiro, P. T. (2019). On logical bifurcation diagrams. J. Theor. Biol. 466, 39–63. doi: 10.1016/j.jtbi.2019.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Abou-Jaoudé, W., Monteiro, P. T., Naldi, A., Grandclaudon, M., Soumelis, V., Chaouiya, C., et al. (2014). Model checking to assess T-helper cell plasticity. Front. Bioeng. Biotechnol. 2:86. doi: 10.3389/fbioe.2014.00086

PubMed Abstract | CrossRef Full Text | Google Scholar

Aghamiri, S. S., Singh, V., Naldi, A., Helikar, T., Soliman, S., and Niarakis, A. (2020). Automated inference of Boolean models from molecular interaction maps using CaSQ. Bioinformatics 1–13. doi: 10.1093/bioinformatics/btaa484. Available online at: https://academic.oup.com/bioinformatics/advance-article

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, A. C., Joller, N., and Kuchroo, V. K. (2016). Lag-3, Tim-3, and TIGIT: co-inhibitory receptors with specialized functions in immune regulation. Immunity 44, 989–1004. doi: 10.1016/j.immuni.2016.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Baumeister, S. H., Freeman, G. J., Dranoff, G., and Sharpe, A. H. (2016). Coinhibitory pathways in immunotherapy for cancer. Annu. Rev. Immunol. 34, 539–573. doi: 10.1146/annurev-immunol-032414-112049

PubMed Abstract | CrossRef Full Text | Google Scholar

Boutillier, P., Maasha, M., Li, X., Medina-Abarca, H. F., Krivine, J., Feret, J., et al. (2018). The Kappa platform for rule-based modeling. Bioinformatics 34, i583–i592. doi: 10.1093/bioinformatics/bty272

PubMed Abstract | CrossRef Full Text | Google Scholar

Brownlie, R. J., and Zamoyska, R. (2013). T cell receptor signalling networks: branched, diversified and bounded. Nat. Rev. Immunol. 13, 257–269. doi: 10.1038/nri3403

PubMed Abstract | CrossRef Full Text | Google Scholar

Chakraborty, A. K. (2017). A perspective on the role of computational models in immunology. Annu. Rev. Immunol. 35, 403–439. doi: 10.1146/annurev-immunol-041015-055325

PubMed Abstract | CrossRef Full Text | Google Scholar

Chylek, L. A., Akimov, V., Dengjel, J., Rigbolt, K. T. G., Hu, B., Hlavacek, W. S., et al. (2014). Phosphorylation site dynamics of early T-cell receptor signaling. PLoS ONE 9:e104240. doi: 10.1371/journal.pone.0104240

PubMed Abstract | CrossRef Full Text | Google Scholar

Das, J., Ho, M., Zikherman, J., Govern, C., Yang, M., Weiss, A., et al. (2009). Digital signaling and hysteresis characterize Ras activation in lymphoid cells. Cell 136, 337–351. doi: 10.1016/j.cell.2008.11.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Eftimie, R., Gillard, J. J., and Cantrell, D. A. (2016). Mathematical models for immunology: current state of the art and future research directions. Bull. Math. Biol. 78, 2091–2134. doi: 10.1007/s11538-016-0214-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Fabregat, A., Sidiropoulos, K., Garapati, P., Gillespie, M., Hausmann, K., Haw, R., et al. (2016). The reactome pathway knowledgebase. Nucleic Acids Res. 44, D481–D487. doi: 10.1093/nar/gkv1351

PubMed Abstract | CrossRef Full Text | Google Scholar

Flobak, Å., Baudot, A., Remy, E., Thommesen, L., Thieffry, D., Kuiper, M., et al. (2015). Discovery of drug synergies in gastric cancer cells predicted by logical modeling. PLoS Comput Biol. 11:e1004426. doi: 10.1371/journal.pcbi.1004426

PubMed Abstract | CrossRef Full Text | Google Scholar

Funahashi, B. A., Matsuoka, Y., Jouraku, A., Morohashi, M., Kikuchi, N., and Kitano, H. (2008). A versatile modeling tool for biochemical networks. Proc. IEEE 96, 1254–1265. doi: 10.1109/JPROC.2008.925458

CrossRef Full Text | Google Scholar

Gebser, M., Guziolowski, C., Ivanchev, M., Schaub, T., Siegel, A., Veber, P., et al. (2010). “Repair and prediction (under inconsistency) in large biological networks with answer set programming,” in Principles of Knowledge Representation and Reasoning: Proceedings of the Twelfth International Conference, KR 2010 (Toronto, ON: AAAI Press).

Google Scholar

Grieco, L., Calzone, L., Bernard-Pierrot, I., Radvanyi, F., Kahn-Perlès, B., and Thieffry, D. (2013). Integrative modelling of the influence of MAPK network on cancer cell fate decision. PLoS Comput. Biol. 9:e1003286. doi: 10.1371/journal.pcbi.1003286

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoops, S., Gauges, R., Lee, C., Pahle, J., Simus, N., Singhal, M., et al. (2006). COPASI - A COmplex PAthway SImulator. Bioinformatics 22, 3067–3074. doi: 10.1093/bioinformatics/btl485

PubMed Abstract | CrossRef Full Text | Google Scholar

Ishida, Y., Agata, Y., Shibahara, K., and Honjo, T. (1992). Induced expression of PD-1, a novel member of the immunoglobulin gene superfamily, upon programmed cell death. EMBO J. 11, 3887–3895.

PubMed Abstract | Google Scholar

June, C. H., Warshauer, J. T., and Bluestone, J. A. (2017). Is autoimmunity the Achilles' heel of cancer immunotherapy? Nat. Med. 23, 540–547. doi: 10.1038/nm.4321

CrossRef Full Text | Google Scholar

Jung, I. Y., Kim, Y. Y., Yu, H. S., Lee, M., Kim, S., and Lee, J. (2018). CRISPR/Cas9-mediated knockout of DGK improves antitumor activities of human T cells. Cancer Res. 78, 4692–4703. doi: 10.1158/0008-5472

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaufman, M., Andris, F., and Leo, O. (1999). A logical analysis of T cell activation and anergy. Proc. Natl. Acad. Sci. U.S.A. 96, 3894–3899. doi: 10.1073/pnas.96.7.3894

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaufman, M., Urbain, J., and Thomas, R. (1985). Towards a logical analysis of the immune response. J. Theor. Biol. 114, 527–561. doi: 10.1016/s0022-5193(85)80042-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Klarner, H., Siebert, H., Nee, S., and Heinitz, F. (2018). Basins of attraction, commitment sets and phenotypes of Boolean networks. IEEE ACM Trans. Comput. Biol. Bioinformatics 17, 1115–1124. doi: 10.1109/TCBB.2018.2879097

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuperstein, I., Bonnet, E., Nguyen, H.-A., Cohen, D., Viara, E., Grieco, L., et al. (2015). Atlas of Cancer Signalling Network: a systems biology resource for integrative analysis of cancer data with Google Maps. Oncogenesis 4:e160. doi: 10.1038/oncsis.2015.19

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Novère, N. (2015). Quantitative and logic modelling of molecular and gene networks. Nat. Rev. Genet. 16, 146–158. doi: 10.1038/nrg3885

PubMed Abstract | CrossRef Full Text | Google Scholar

Leach, D. R., Krummel, M. F., and Allison, J. P. (1996). Enhancement of antitumor immunity by CTLA-4 blockade. Science 271, 1734–1736. doi: 10.1126/science.271.5256.1734

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopez, C. F., Muhlich, J. L., Bachman, J. A., and Sorger, P. K. (2013). Programming biological models in Python using PySB. Mol. Syst. Biol. 9, 1–19. doi: 10.1038/msb.2013.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Macian, F. (2005). NFAT proteins: key regulators of T-cell development and function. Nat. Rev. Immunol. 5, 472–484. doi: 10.1038/nri1632

PubMed Abstract | CrossRef Full Text | Google Scholar

Mi, H., Muruganujan, A., and Thomas, P. D. (2013). PANTHER in 2013: modeling the evolution of gene function, and other gene attributes, in the context of phylogenetic trees. Nucleic Acids Res. 41, 377–386. doi: 10.1093/nar/gks1118

PubMed Abstract | CrossRef Full Text | Google Scholar

Miskov-Zivanov, N., Zuliani, P., Wang, Q., Clarke, E. M., and Faeder, J. R. (2016). “High-level modeling and verification of cellular signaling,” in 2016 IEEE International High Level Design Validation and Test Workshop, HLDVT 2016 (Santa Cruz, CA), 162–169.

Google Scholar

Monteiro, P. T., and Chaouiya, C. (2012). “Efficient verification for logical models of regulatory networks,” in 6th International Conference on Practical Applications of Computational Biology & Bioinformatics (Berlin; Heidelberg: Springer), Vol. 154, 259–267.

Google Scholar

Myers, G. J. (1979). The Art of Software Testing, 1st Edn. New Jersey, NJ: Wiley Publishing.

Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Naldi, A., Hernandez, C., Abou-Jaoudé, W., Monteiro, P. T., Chaouiya, C., and Thieffry, D. (2018a). Logical modeling and analysis of cellular regulatory networks with GINsim 3.0. Front. Physiol. 9:646. doi: 10.3389/fphys.2018.00646

PubMed Abstract | CrossRef Full Text | Google Scholar

Naldi, A., Hernandez, C., Levy, N., Stoll, G., Monteiro, P. T., Chaouiya, C., et al. (2018b). The CoLoMoTo interactive notebook: accessible and reproducible computational analyses for qualitative biological networks. Front. Physiol. 9:680. doi: 10.3389/fphys.2018.00680

PubMed Abstract | CrossRef Full Text | Google Scholar

Oyeyemi, O. J., Davies, O., Robertson, D. L., and Schwartz, J. M. (2015). A logical model of HIV-1 interactions with the T-cell activation signalling pathway. Bioinformatics 31, 1075–1083. doi: 10.1093/bioinformatics/btu787

PubMed Abstract | CrossRef Full Text | Google Scholar

Perley, J., Mikolajczak, J., Buzzard, G., Harrison, M., and Rundell, A. (2014). Resolving early signaling events in T-Cell activation leading to IL-2 and FOXP3 transcription. Processes 2, 867–900. doi: 10.3390/pr2040867

CrossRef Full Text | Google Scholar

Ribas, A., and Wolchok, J. D. (2018). Cancer immunotherapy using checkpoint blockade. Science 359, 1350–1355. doi: 10.1126/science.aar4060

PubMed Abstract | CrossRef Full Text | Google Scholar

Riese, M. J., Moon, E. K., Johnson, B. D., and Albelda, S. M. (2016). Diacylglycerol kinases (DGKs): novel targets for improving T cell activity in cancer. Front. Cell Dev. Biol. 4:108. doi: 10.3389/fcell.2016.00108

PubMed Abstract | CrossRef Full Text | Google Scholar

Robert, C., Miller, W. H., Francis, S., Chen, T.-T., Ibrahim, R., Hoos, A., et al. (2011). Ipilimumab plus dacarbazine for previously untreated metastatic melanoma. N. Engl. J. Med. 364, 2517–2526. doi: 10.1056/nejmoa1104621

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodriguez-Jeorge, O., Kempis-Calanis, L., Abou-Jaoude, W., Hernandez, C., Ramirez-Pliego, O., Thomas-Chollier, M., et al. (2019). Cooperation between T cell receptor and Toll-like receptor 5 signaling for CD4+ T cell activation. Sci. Signal. 12:eaar3641. doi: 10.1126/scisignal.aar3641

PubMed Abstract | CrossRef Full Text | Google Scholar

Roose, J. P., Mollenauer, M., Ho, M., Kurosaki, T., and Weiss, A. (2007). Unusual interplay of two types of Ras activators, RasGRP and SOS, establishes sensitive and robust Ras activation in lymphocytes. Mol. Cell. Biol. 27, 2732–2745. doi: 10.1128/MCB.01882-06

PubMed Abstract | CrossRef Full Text | Google Scholar

Saadatpour, A., Albert, R., and Reluga, T. C. (2013). A reduction method for Boolean network models proven to conserve attractors. SIAM J. Appl. Dynam. Syst. 12, 1997–2011. doi: 10.1137/13090537x

CrossRef Full Text | Google Scholar

Sanchez-Villanueva, J., Rodr-guez-Jorge, O., Ram-rez-Pliego, O., Rosas-Salgado, G., Abou-Jaoude, W. C. H., Naldi, A., et al. (2019). Contribution of ROS and metabolic status to neonatal and adult CD8+T cell activation. PLoS ONE 14:e0226388. doi: 10.1371/journal.pone.0226388

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarma, G. P., Jacobs, T. W., Watts, M. D., Vahid Ghayoomie, S., Larson, S. D., and Gerkin, R. C. (2016). Unit testing, model validation, and biological simulation [version 1; referees: 2 approved, 1 approved with reservations]. F1000Res. 5, 1–17. doi: 10.12688/F1000RESEARCH.9315.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Simpson, T. R., Li, F., Montalvo-Ortiz, W., Sepulveda, M. A., Bergerhoff, K., Arce, F., et al. (2013). Fc-dependent depletion of tumor-infiltrating regulatory T cells co-defines the efficacy of anti-CTLA-4 therapy against melanoma. J. Exp. Med. 210, 1695–1710. doi: 10.1084/jem.20130579

PubMed Abstract | CrossRef Full Text | Google Scholar

Slenter, D. N., Kutmon, M., Hanspers, K., Riutta, A., Windsor, J., Nunes, N., et al. (2018). WikiPathways: a multifaceted pathway database bridging metabolomics to other omics research. Nucleic acids Res. 46, D661–D667. doi: 10.1093/nar/gkx1064

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith-Garvin, J. E., Koretzky, G. A., and Jordan, M. S. (2009). T cell activation. Annu. Rev. Immunol. 27, 591–619. doi: 10.1146/annurev.immunol.021908.132706

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Traynard, P., Fauré, A., Fages, F., and Thieffry, D. (2016). Logical model specification aided by model-checking techniques: application to the mammalian cell cycle regulation. Bioinformatics 32, i772–i780. doi: 10.1093/bioinformatics/btw457

PubMed Abstract | CrossRef Full Text | Google Scholar

Ventimiglia, L. N., and Alonso, M. A. (2013). The role of membrane rafts in Lck transport, regulation and signalling in T-cells. Biochem. J. 454, 169–179. doi: 10.1042/BJ20130468

PubMed Abstract | CrossRef Full Text | Google Scholar

Voisinne, G., Kersse, K., Chaoui, K., Lu, L., Chaix, J., Zhang, L., et al. (2019). Quantitative interactomics in primary T cells unveils TCR signal diversification extent and dynamics. Nat. Immunol. 20, 1530–1541. doi: 10.1038/s41590-019-0489-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Walunas, T. L., Lenschow, D. J., Bakker, C. Y., Linsley, P. S., Freeman, G. J., Green, J. M., et al. (1994). CTLA-4 can function as a negative regulator of T cell activation. Immunity 1, 405–413. doi: 10.1016/1074-7613(94)90071-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Miskov-Zivanov, N., Liu, B., Faeder, J. R., Lotze, M., and Clarke, E. M. (2016). “Formal modeling and analysis of pancreatic cancer microenvironment,” in Computational Methods in Systems Biology, eds E. Bartocci, P. Lio, and N. Paoletti (Heidelberg: Springer International Publishing), 289–305.

Google Scholar