Dynamical modeling and analysis of large cellular regulatory networks

12 Pages·2017·0.93 MB·English

Dynamical modeling and analysis of large cellular
regulatory networks
D Bérenguier, C Chaouiya, Pedro T. Monteiro, A Naldi, E Remy, Denis
Thieffry, L Tichit
To cite this version:
D Bérenguier, C Chaouiya, Pedro T. Monteiro, A Naldi, E Remy, et al.. Dynamical modeling and
analysis of large cellular regulatory networks. Chaos: An Interdisciplinary Journal of Nonlinear Sci-
ence, 2013, 23 (2), ￿10.1063/1.4809783￿. ￿hal-01261849￿
HAL Id: hal-01261849
https://hal.science/hal-01261849
Submitted on 25 Jan 2016
HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est
archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
Dynamical modeling and analysis of large cellular regulatory networks.
D. B´erenguier,1 C. Chaouiya,2,a) P.T. Monteiro,2,3 A. Naldi,4 E. Remy,1,b) D.Thieffry,5,c) and L.Tichit1,6
1)Institut de Math´ematiques de Luminy, Marseille, France
2)Instituto Gulbenkian de Ciˆencia, Oeiras, Portugal
3)Instituto de Engenharia de Sistemas e Computadores - Investiga¸ca˜o e Desenvolvimento (INESC-ID), Lisboa,
Portugal
4)Center for Integrative Genomics, University of Lausanne, Lausanne, Switzerland
5)Institut de Biologie de l’Ecole Normale Sup´erieure, Paris, France
6)Aix-Marseille University, Marseille, France
(Dated: 16 April 2013)
The dynamical analysis of large biological regulatory networks requires the development of scalable methods
for mathematical modeling. Following the approach initially introduced by R. Thomas, we formalize the
interactions between the components of a network in terms of discrete variables, functions and parameters.
Model simulations result in directed graphs, called state transition graphs. We are particularly interested
in reachability properties and asymptotic behaviors, which correspond to terminal strongly connected com-
ponents (or ”attractors”) in the state transition graph. A well-known problem is the exponential increase
of the size of state transition graphs with the number of network components, in particular when using the
biologically realistic asynchronous updating assumption. To address this problem, we have developed several
complementary methods enabling the analysis of the behavior of large and complex logical models: (i) the
definition of transition priority classes to simplify the dynamics; (ii) a model reduction method preserving
essentialdynamicalproperties,(iii)anovelalgorithmtocompactstatetransitiongraphsanddirectlygenerate
compressedrepresentations,emphasizingrelevanttransientandasymptoticdynamicalproperties. Thepower
of an approach combining these different methods is demonstrated by applying them to a recent multilevel
logicalmodelforthenetworkcontrollingCD4+Thelpercellresponsetoantigenpresentationandtoadozen
cytokines. This model accounts for the differentiation of canonical Th1 and Th2 lymphocytes, as well as of
inflammatory Th17 and regulatory T cells, along with many hybrid subtypes. All these methods have been
implementedintothesoftwareGINsim,whichenablesthedefinition,theanalysisandthesimulationoflogical
regulatory graphs.
Keywords: logical modeling, discrete attractors, T-helper cell differentiation, state transition graph, model
reduction, model-checking
The dynamical analysis of comprehensive biologi- has proved effective in its application to a variety of
cal regulatory networks requires the development regulatory and signaling networks, from yeast cell cycle
of scalable mathematical modeling methods. In control15 to T lymphocyte differentiation27.
this context, discrete (Boolean or multi-valued)
The logical modeling approach has been implemented
logical modeling is increasingly used to handle
into the software GINsim, which enables the definition
and analyze large molecular networks18. This ar-
of logical regulatory graphs and provides a number of
ticle focuses on the presentation of several ap-
original functionalities. These include the construction
proaches to cope with the inherent exponential
of synchronous or asynchronous state transition graphs
growth of the discrete state space as the size of
that represent model dynamical behaviors, along with
the regulatory networks considered increases.
algorithms enabling the determination of all logical sta-
ble states and the analysis of the roles of regulatory
circuits9,26. However, when focusing on transient as-
pects of the dynamics or on the reachability of the at-
I. INTRODUCTION
tractorsfromspecificinitialconditions,wearefacingthe
recurrentcombinatorialexplosioninherentinthesemod-
To model biological regulatory networks, we rely on
els: the size of the state space grows exponentially with
the logical approach initially introduced by R. Thomas,
the number of regulatory components involved in the
where the interactions between the components of a net-
model. This problem is particularly acute in the case of
work are formalized in terms of discrete variables, func-
asynchronous, non-deterministic updating mode, which
tions and parameters10,37,38. This modeling approach
is usually more biologically realistic than the simpler,
deterministic synchronous updating mode17,37. Here, we
presentanoverviewofthreemaincomplementarystrate-
gies to cope with this combinatorial explosion.
a)Electronicmail: [email protected]
b)Electronicmail: [email protected] Thefirststrategyconsistsinreducingthemodelbefore
c)Electronicmail: thieff[email protected] performing simulations or other kinds of analysis29.
2
The second strategy simplifies the state transition Logicalfunctions
graphsbyforcingchoicesbetweenalternativetransitions; (cid:0)☎ (cid:0)✁✂ KCI(s)=2 ¬Cro ∨ CII
thiscanbeachievedbydefiningpriority(a/synchronous) KCro(s)=3 ¬CI:2 ∧ ¬Cro:3
transition classes, which are similar to time-scale based KCro(s)=2 ¬CI:2 ∧ Cro:3
assumptions often used to simplify the dynamical analy- KCII(s)=1 ¬CI:2 ∧ ¬Cro:3 ∧ N
sis of ODE models14.
The third, novel strategy consists in compressing the ✄ (cid:0)☎☎ KN(s)=1 ¬CI ∧ ¬Cro:2
state transition graph into a novel graph representation,
FIG.1. Logical regulatory graph of the bacteriophage
called hierarchical transition graph, which keeps track lambdaswitch35. Left: theinteractiongraph,withthefour
of attractors and their basins of attraction, as well as componentsCI,Cro,CIIandN.Right: thelogicalfunctions,
of transient oscillatory properties; here, we further pro- s denoting the vector (s ,s ,s ,s ) of the component
CI Cro CII N
pose an algorithm for the construction of this hierarchi- levels. For legibility, we rewrite each rule in terms of logical
cal graph. We also show that this method can be used variables, e.g. CI denotes an interaction going out CI with
in combination with the two aforementioned approaches a threshold 1, and it is true whenever sCI ≥ 1, while CI:2
denotes an interaction going out CI with a threshold 2; it is
to get insights into the dynamics of complex logical reg-
truewhenevers ≥2. Here,foreachcomponent,weprovide
ulatory graphs. CI
the rule(s) leading to a non-zero value of the logical function
In addition, model-checking approaches rely on sym-
(meaning that when none of these conditions is fulfilled, the
bolicrepresentationsofthedynamics, exploringonlythe
value is 0). For instance, the rule for K (s)=2 is satisfied
CI
necessarystatespacerequiredfortheverificationofprop- for 30 states (those such that s = 0 or s = 1); for all
Cro CII
erties expressed as temporal logic formulas. other states, CI’s target value is 0. Note that values 1 of CI
SectionIIintroducesthebasicsofthemultilevellogical and Cro are always transient for this set of rules.
formalism and provides an overview of selected methods
enablingtheanalysisofthedynamicsoflargelogicalreg-
ulatory networks. • G = {gi}i=1...n is the set of n regulatory compo-
The definition of hierarchical transition graphs is at nents. Each component gi is associated with a dis-
the core of Section III, referring to the relatively sim- crete variable si in Di ={0,...,maxi}. A state is
pleexampleofthebacteriophagelambdacoreregulatory thus defined as a vector s∈S =Πgj∈G Dj.
network.
• K ={K } is the set of logical functions;
Section IV takes advantage of a recent comprehensive i i=1...n
K :S →D defines for each state, the target level
modeloftheregulatorynetworkcontrollingT-helpercell i i
of g .
differentiationinresponsetoantigenpresentationandto i
a dozen cytokines27 to illustrate the power of the com- The arcs are deduced from the functions in K; there is
pression of state transition graphs into hierarchical tran- a regulatory interaction from g to g iff there are two
i j
sition graphs, as well as the insights gained into the cor- states s and s, differing only by the value of g , that lead
(cid:48) i
responding logical dynamical behavior. to different values of K :
i
Finally, Section V proposes some global conclusions
and discusses current challenges and further prospects. ∃s,s(cid:48) ∈S s(cid:48)k =sk∀k (cid:54)=i, andsi (cid:54)=s(cid:48)i,s.t.Kj(s)(cid:54)=Kj(s(cid:48)).
Figure 1 illustrates this definition with a logical regu-
latory graph for the bacteriophage lambda switch.
II. LOGICAL MODELING AND ANALYSIS OF
The dynamics of logical models are represented in the
REGULATORY GRAPHS
formofstatetransitiongraphsasdefinedinthefollowing
subsection.
Thissectionintroducesthebasicsofthelogicalformal-
ism and presents a short overview of existing methods
that enhance the dynamical analysis of logical models.
B. State Transition Graphs
Definition 2. Given a logical regulatory graph R =
A. Logical Regulatory Graphs
(G, K), its(full)StateTransitionGraph(STG),denoted
by E =(S, T), is a directed graph with:
A logical model is defined by an interaction graph
where the nodes denote regulatory components (genes, • S the state space of R: S =Π D ,
proteins, etc.) and the arcs denote regulatory interac- gj∈G j
tions. Moreover, a discrete variable is associated to each • T : S2 →{0,1}thetransitionfunction: thereisan
component,accountingforitslevelofactivity(orexpres- arcconnectingastatestoitssuccessors(cid:48) whenever
sion). Logical functions define the dynamical evolution T(s,s(cid:48))=1.
of the model.
The transition function is defined according to an up-
Definition 1. A LogicalRegulatoryGraph R=(G, K) dating policy, which indicates the components to be up-
is a graph, where dated in each transition. Here, for sake of brevity, we
3
only consider the asynchronous updating policy (Defini- Hence, A can be reached from any state in B or in
∗ A∗
tion 3). All results could be extended to other updating B ; no other attractor can be reached from any state
A∗
policies (including mixed (a)synchronous priority classes in B .
A∗
as presented in Section IIC2).
Definition 3. Given a logical regulatory graph R =
C. Coping with large dynamics
(G, K), the transition function defined according to the
asynchronous updating policy specifies, for each state s,
Given a logical regulatory graph, the associated state
its successor states (as many as the components called to space has (cid:81) |D | elements (i.e. 2G, in the case of
update in s): ∀(s,s(cid:48))∈S2, Booleanvariagib∈lGes),mi eaningthatitssiz|e|growsexponen-
1, if ∃g ∈G s.t. K (s)(cid:54)=s , tially with the number of regulatory components. Most
T(s,s(cid:48))=0,so(cid:48)ith=erswiii+se.|KKii((ss))−−ssiii| and ∀igj (cid:54)=gi,s(cid:48)j =sj pthriospperrtoibelsemarebythluessseNnPin-gcotmheplseitzee,obfuttheonseeacracnh smpiatcige.ate
Here,webrieflyreviewstrategiestoeasetheanalysisof
largedynamics. Afirstapproachconsistsinreducingthe
Note that updated are performed stepwise and thus
model, whileensuringthepreservationofkeyproperties.
transitions connect neighbouring states (i.e. their Ham-
Anotherstrategylessensthenumberoftransitionsofthe
ming distance is 1).
STG (hence simplifying the dynamics) assigning priori-
Weareofteninterestedinasub-graphofthefullSTG,
tiestoupdatingcalls,relyingonbiologicallywell-founded
which is generated considering a (set of) initial state(s).
assumptions. Othermethodsenablethereductionofthe
Then, the property of interest relates to the attractors
size of a STG, either by compacting it without losing
reachable from this (set of) initial condition(s). Fig-
any information, by applying appropriate reductions, or
ure 2A displays the STG obtained for the phage lambda
byconsideringalternaterepresentations. Finally, weend
model,startingfromastatewhereallvariablesaresetto
this section with a short discussion on model-checking
zero. Attractors, which denote asymptotical behaviors,
applied to multilevel logical models.
are defined in a STG as its terminal strongly connected
components(SCCs). Recallthatstronglyconnectedcom-
ponents are defined as the maximal strongly connected
1. Model reduction
subgraphs (i.e. there is a path from each node to every
other node)5,13.
GivenE =(S,T)aSTG,weintroducefurthernotation A first strategy to reduce the complexity of a model is
below. to reduce its size, by removing some components. This
is often done manually by the modeler, defining direct
• Sccisthesetofthestronglyconnectedcomponents interactions even when it is known that the regulatory
(SCCs) of E ; effects involve intermediate components. Obviously, by
lesseningthenumberofcomponents,suchreductionslead
• ∀s,s ∈ S, s (cid:32) s means that there exists a path
(cid:48) (cid:48) to smaller state spaces and hence simplified dynamics.
fromstos (weconsiderthatasequenceofaunique
(cid:48) We have proposed to automate such model reductions
state forms a path of length 0, hence s(cid:32)s);
and characterized their impact on the dynamics29. Basi-
• ∀s ∈ S, ∀C ∈ Scc, s (cid:32) C means that there exists cally,thereductionofacomponentamountstoattribute
a path from s to any state s ∈C; its regulatory role to its own regulators and to modify
(cid:48)
the logical function driving the behavior of its targets
• S is the set of the trivial SCCs (i.e. reduced to a accordingly. The reduction of a self-regulated compo-
single node) S ={{s}∈Scc,s∈S}; nent is forbidden for it would not fit this rationale and
wouldchangethedynamicalpropertiesofthemodel. In-
• C is the set of the complex SCCs: C = {C ∈
deed, we could prove that, provided this restriction on
Scc,#C ≥2} (or C =Scc\S);
self-regulated components, the stable states and elemen-
tary terminal cycles of a reduced model exactly corre-
• The sign denotes terminal elements of Scc that
∗
spond to those of the original model. Moreover, the re-
will be referred to as attractors: C is terminal iff
∗
duced model displays at least as many complex attrac-
∀s ∈ C , ∀s ∈ S, T(s,s) = 1 ⇒ s ∈ C . The
∗ (cid:48) (cid:48) (cid:48) ∗
tors, some corresponding to complex attractors of the
non-terminal components are transient. In addi-
tion,C (resp. S )denotesthesetofthecomplex original model, while others correspond to original tran-
∗ ∗
sient oscillatory behaviors. In short, the main property
attractors (resp. the set of stable states).
of the proposed reduction is that it may suppress some
Definition 4. Let A∗ be an attractor, we define BA∗ the transitions or paths, but never generates new ones.
basin of attraction of A∗: BA∗ = {s ∈ S,s (cid:32) A∗}. We Considering signaling networks including non-
further define BA∗, the strict basin of attraction of A∗, regulated input components, which usually account for
external stimuli, Saadatpour et al. recently proposed to
BA∗ ={s∈BA∗ s.t. ∀X∗ ∈(C∗∪S∗)\{A∗}, s(cid:54)∈BX∗}. reduceinputcascadesthatstabilizeunderconstantinput
4
conditions32. This reduction has obviously no impact consideredas“uncontrolled”variables,whichareallowed
on the number and nature of the attractors, although tofreelyvaryateachtimestep. Anaturalreductioncon-
it might change their reachability. Another type of sists in projecting the state space on the set of internal
components that can be ignored are (pseudo-)output componentsandlabellingeachtransitionwiththevalues
components that have no outgoing interactions or that of the input components that enable that transition (for
only regulate (pseudo-)output components28. Indeed, more details, see 24).
such output cascades have no impact, neither on the Other strategies, mainly developed by the formal ver-
number and nature of the attractors, nor on their ification community, reduce the state space yet ensur-
reachability. ing that truth values of (linear) temporal logic formulas
are preserved. This is the case of partial-order reduction
methods that basically consist in identifying, for each
2. Priority classes state, a subset of transitions to explore (hence not ex-
ploring all the successors). Alternative (rather similar)
Asynchronous state transitions graphs can be some- definitions of these sets have been proposed, called stub-
times simplified by reducing the number of transitions, born, ample or persistent sets12,19.
using relatively simple temporal assumptions. Indeed, Relying on the Petri net representation of logical reg-
in all states, the asynchronous scheme defines as many ulatory graphs8, and using Petri net tools (e.g. TINA6),
transitions as the number of components called to up- wehaverecentlyappliedsuchapartial-orderreductionto
date, thus potentially generating spurious trajectories. checkreachabilitypropertiesonalargelogicalmodel(en-
A number of these can be ignored by defining priority compassing 72 regulatory components). For this specific
classes ranking updating calls14. When two calls with model,duetothestructureofitsdynamics,partial-order
distinct priority ranks are enabled in a state, the one reduction proved to be poorly effective. However, there
with the lowest rank is discarded. Updatings belonging is certainly room for improvements of these methods16,
to the same class can be treated synchronously or asyn- and further work might identify a class of logical graphs
chronously. In GINsim, it is thus possible to partition more amenable to this kind of reductions.
component updatings into distinct classes that imple-
mentsuchapriorityscheme9,26. Needlesstosay,priority
classesshouldbebiologicallywell-foundedtoensurethat 4. Model-checking
discarded trajectories are indeed irrelevant.
Duringtherecentyears, formalverificationtechniques
based on model-checking have been successfully applied
3. Lessening the size of the state transition graph totheanalysisofmolecularnetworkmodels4,7,24,25. This
approach is directly applicable to the verification of log-
Severalstudieshaveaddressedtheproblemofthecom- ical regulatory graphs, which constitute a class of finite
binatorial explosion of the state space of asynchronous state systems. In general, experimentally observed bio-
transition systems. logical behaviors can be expressed in terms of temporal
Given a state transition graph (STG), an informative logic statements, and model-checking algorithms used to
view of the dynamics is provided by the graph of its automaticallyverifyifamodelsatisfiesthesestatements.
stronglyconnectedcomponents(SCCs),whereeachnode Whenusingexplicitrepresentationsofstatesandtran-
accounts for one SCC (possibly keeping the information sitions, model-checking may use partial-order reduction
of the states it encompasses). The resulting graph is to lessen the size of the search space. However, symbolic
a directed acyclic graph, which is often much smaller model-checkingreliesonimplicitrepresentations, scaling
than the original STG, yet keeping all the reachability better for large models. The choice of the temporal logic
information (see Figure 2A-B). Tarjan defined an effi- depends on the type of property to be checked12. Here,
cientalgorithmtocomputetheSCCsofadirectedgraph we are mainly interested in attractor reachability from a
(linear in the number of nodes and arcs)34. Tournier (setof)initialcondition(s)aswellasintheconditionsen-
and Chaves39 have already applied SCC decomposition abling such trajectories. This supposes a previous char-
to STGs. However, SCC compaction remains limited in acterization of the attractors, among which the stable
the case of networks with long or numerous regulatory states can be efficiently identified beforehand9.
cascades,whichgiverisetomultiplelinear(althoughpo- GINsim includes an export converting logical mod-
tentially branching) pathways in the resulting STGs. els into NuSMV symbolic descriptions24. NuSMV is a
AnotherapproachthatkeepsthewholeSTGstructure symbolic model-checking tool capable of verifying finite
applies to models that encompass a significant number state machines against a set of property requirements,
of input components. These account for external stimuli expressedastemporallogicformulas11. Thisexportsup-
(e.g. environmental cues) and the corresponding vari- ports the definition of priority classes, and takes advan-
ablesaregenerallymaintainedconstant. Inthiscase,the tage of the reduction over input components evoked in
STGs corresponding to different combinations of input Section IIC3, these being specified either as constant or
values are disconnected. Input components may also be as freely varying variables. Noteworthy, in the case of
5
varying inputs, the notion of stable states needs to be Property 1. 1. A path connecting any HTG compo-
extended: astatemaybestableforgivenvaluesofinput nent to a non-irreversible component implies the
components,andnotforothers28. Formodelswithinput existence of a path in the corresponding STG:
components, it is thus possible to analyze switches be-
tweencellulartypes(i.e. stablestates)andverifythecor- ∀C ∈C ∪I, C ∈C,
(cid:48)
responding input component variations (see Section IV
and Figure 6). C (cid:32)H C(cid:48) =⇒s(cid:32)s(cid:48) ∀s∈C, ∀s(cid:48) ∈C(cid:48).
2. A path between two states implies the existence of a
III. HIERARCHICAL TRANSITION GRAPH path between the HTG components they belong to:
This section deals with the definition and properties ∀s,s(cid:48) ∈S,s(cid:32)s(cid:48) =⇒C (cid:32)H C(cid:48),withs∈C,s(cid:48) ∈C(cid:48).
of a novel, compact hierarchical graph, where a set of
states is shrunken into a single node, whenever it forms
astronglyconnectedcomponent(SCC),ora(setof)lin-
ear chain(s) leading to the same set of SCCs and attrac- Proof. 1. Let C (cid:32)H C(cid:48), with C ∈ C ∪ I and C(cid:48) ∈
tors. ComparedtotheSCCgraphmentionedabove,this C ∪S . Then C ∈ σ(C) : ∀s ∈ C, s (cid:32) C and
∗ (cid:48) (cid:48)
graphgenerallycorrespondstoafurtherreductionofthe the first item of Property 1 is proved by definition.
statetransitiongraph(STG).Furthermore,theresulting
grouping of states greatly eases the interpretation of the 2. Let s, s(cid:48) ∈ S s.t. s (cid:32) s(cid:48), and denote C and C(cid:48)
structureofthedynamicsintermsofbasinsofattraction. the components of the HTG such that s ∈ C and
s ∈C .
(cid:48) (cid:48)
• If C =C , the statement is obviously true.
A. Definitions (cid:48)
• If C (cid:54)= C , let (s = s ,s ,...s = s) be the
(cid:48) 1 2 k (cid:48)
Let us first define the application σ that associates to pathfromstos intheSTG:∀i=1,...k−1,
(cid:48)
each SCC C, the set of SCCs, complex or terminal, that s ∈ S and T(s ,s ) = 1; if ∀i = 1,...k, C
i i i+1 i
are reachable from C, including C itself if it is complex denotes the component of Scc such that s ∈
i
or terminal: C , we have T(C ,C ) = 1 or C = C .
i i i+1 i i+1
σ(X)=(cid:8)C ∈C ∪S s.t. X =C or ∀ s∈X, s(cid:32)C(cid:9). Hence, following the path s (cid:32) s, we obtain
∗ (cid:48)
Furthermore, we define I ⊂ 2S, the set of irreversible that C (cid:32)H C(cid:48).
transient componentsinwhichtrivialnon-terminalSCCs
(elements of S \S ) that have the same σ-image are
∗
grouped together:
Remark 1. Property 1 does not ensure equivalence of
I =(cid:8)I ∈2S s.t. ∀s∈I, {s}∈S \S and path existence in STG and related HTG. Indeed, in item
∗
s, s ∈I ⇒σ({s})=σ({s})}. 1, we have the restriction that C(cid:48) (cid:54)∈ I: when C (cid:32)H C(cid:48),
(cid:48) (cid:48)
with C ∈I, given s∈C, we cannot ensure the existence
(cid:48)
Definition 5. Given a STG E =(S,T), we denote H= of a path in the STG from s to a state s(cid:48) ∈C(cid:48).
(C∪I∪S∗,T)itscorrespondingHierarchicalTransition In Figure 2 we have such a situation, where C (cid:32)H C(cid:48),
Graph(HTG),whereT : (C∪I∪S∗)2 →{0,1}defines with C(cid:48) ∈ I, and ∃s ∈ C s.t. there is no path in the
the arcs of H, STG from s to any state in C . Indeed, considering Fig-
(cid:48)
ure 2, panel C, the irreversible component i#7 contains
T(C,C )=1⇐⇒∃s∈C, ∃s ∈C s.t. T(s,s)=1.
(cid:48) (cid:48) (cid:48) (cid:48) the state 1000 (see panel B), and the arc from i#7 to
i#3 indicates that there exists a state s ∈ i#7 and a
EachcomplexSCCoftheSTGiscontractedtoasingle state s ∈ i#3 such that s (cid:32) s (e.g. T(1011,2011)=1,
(cid:48) (cid:48)
nodeintheHTG.Similarly,asingleHTGnodeaccounts
in panel A or B). However, there is no path from state
for all trivial SCCs sharing the same σ-image. Figure 2
1000 to any state of i#3 as illustrated in panel D.
provides an illustration of the HTG construction for the
lambda phage model. Another typical situation for which we have C (cid:32)H C(cid:48),
C ∈ I, and no path in the STG from s ∈ C to s ∈ C ,
(cid:48) (cid:48) (cid:48)
may arise when a hierarchical (irreversible) component
contains disconnected states.
B. Properties
We propose an algorithm to generate HTGs of logical
FortwocomponentsC, C(cid:48) ∈C ∪I∪S∗, thenotation regulatory graphs, given a (set of) initial condition(s).
C (cid:32)H C(cid:48) indicates the existence of a path from C to C(cid:48) Described in the supplementary file43, this algorithm is
in the HTG. The following property relates paths in the basedonTarjan’smethod34 andcompactsaSTGon-the-
STG to paths in the corresponding HTG. fly.
6
A B
0310 1310 2310 ct#2 ct#31
0311 1311 2311
0300 1300 2300
ca#2
0301 1301 2301
0210 1210 2210 i-1010 i-2010
0211 1211 2211
0200 1200 2200 i-0011 i-1011 i-2011
i-0000 i-1000 ss-2000
0201 1201 2201
i-0001 i-1001 i-2001
1110 2110 C D
0111 1111 2111 CI- CII-
0100 1100 2100 ct#2 Cro+ N- ct#31 ct#2 ct#31
0101 1101 2101 CI- Cro+CNI-I-
1010 2010 ca#2 Cro- Cro+ Cro- ca#2
i#7 i#7
0011 1011 2011
0000 1000 2000 CI+ CI+
0001 1001 2001 ss-2000 CNI-I- i#3 ss-2000 i#3
FIG. 2. Lambda phage model: different views of the dynamics. All the graphs, except that of panel (D), have been
generated starting from the initial state (CI, Cro, CII, N) = 0000. (A) State transition graph (STG) with the unique stable
state indicated in a rectangular node; states sharing the same color belong to the same strongly connected component (SCC).
(B) Corresponding graph of the SCCs (same coloring as in panel A). (C) Corresponding hierarchical transition graph (HTG);
arc labels refer to transitions in the underlying STG (i.e. updates of regulatory components). (D) HTG obtained using 1000
as initial state, which belongs to component i#7 in panel C; note the absence of a path from this state to the component i#3
(cf. Remark 1). In the SCC and HTG graphs, node labels indicate the nature of the components: i irreversible; ct complex
transient SCC; ss stable state; ca complex attractor; followed by the numbers of states in the component. For components
reduced to a single state, the value of this state is displayed.
C. Basins of attraction IV. APPLICATION TO TH CELL DIFFERENTIATION
T-helperlymphocytesplayakeyroleintheregulation
Aclassicalwaytostudythedynamicsistofocusonat-
of the immune response in vertebrate. Various T-helper
tractors and their basins of attraction (cf. Definition 4).
subtypes (Th1, Th2, Th17, Treg) have been identified
When using the synchronous dynamics, their computa-
over the years, characterized by the expression of spe-
tion is facilitated by the fact that all states have at most
cific transcription factors and cytokines, which have a
one successor (for more details, see 42). But in the case
critical influence on the selection of specific immune re-
of concurrent behavior, it is computationally much more
sponses, driving pro-inflammatory or allergic responses,
costly(see1,39,and41forthefullyasynchronouscase).
promotingalternativeantibodyclasses,oryetpreventing
Byconstruction,theHTGnodesgrouptogetherstates (auto)immunity by inhibiting the activation and prolif-
of the STG and thus allow to easily recover the basins of eration of other cells.
attractions. Indeed, given A∗ an attractor in the STG, Several modeling studies have been proposed to shed
and C ∈C ∪I∪S∗ a node of the HTG, the states of C light on the regulatory network controlling T-helper cell
are in BA∗, the basin of attraction of A∗, iff A∗ ∈ σ(C). activationanddifferentiation(seee.g. 20,21,23,33,and
The states of C are in BA∗, the strict basin of attraction 40 and references therein).
of A , iff σ(C) ∩ (C ∪ S ) = {A }. Hence, for all
∗ ∗ ∗ ∗ To gain insight into the heterogeneity and the plas-
attractors, it is much easier to identify their basins of
ticity of late T-helper lineages, we have recently built
attraction on the basis of the HTG.
an integrated logical model of the core regulatory net-
Irreversible decisions are taken at the intersections of work and main signaling pathways controlling Th cell
the basins of attraction. Given two consecutive nodes differentiation27 (Figure 3). Encompassing 65 compo-
X and X in the HTG (T(X ,X ) = 1), crucial deci- nents(including13inputs,correspondingtoantigenpre-
1 2 1 2
sionscanbeassociatedwiththearclinkingthesenodesif sentationandadozendifferentcytokines),thismultilevel
σ(X )∩(C ∪S )(cid:54)=σ(X )∩(C ∪S )(i.e. thesystem logical model proved to be too complex to be straight-
1 ∗ ∗ 2 ∗ ∗
enters a more restricted basin of attraction). Consider forwardlysimulated. Thissituationmotivatedthedevel-
twoattractorsA andA suchthatA ∈σ(X )∩σ(X ), opment of the reduction method mentioned above27,29.
∗1 ∗2 ∗1 1 2
A ∈ σ(X ) and A (cid:54)∈ σ(X ). We say that the transi- In the case of our T-helper model, we were led to hide
∗2 1 ∗2 2
tion (X ,X ) belongs to the boundary of B : removing 31 components (shown in grey in Figure 3) and thus ob-
1 2 A∗
all such transitions would isolate B from2the rest of tained a more compact model encompassing 34 nodes
A∗
theHTG.InFigure2(C),transitions(2ct#31, ct#2)and (including the same 13 inputs).
(ct#31, ca#2) constitute the boundary of Bss-2000. Using the resulting reduced T-helper model, we have
7
IFNB_e IFNG_e IL27_e IL6_e IL21_e IL23_e IL10_e TGFB_e IL12_e IL4_e IL15_e IL2_e APC
proliferation
CGC
IFNGR2 IL15RA
GP130 IL2RA
IL12RB1
IFNGR1 IL27RA IL6RA IL10RA IL12RB2 IL4RA IL2RB
IL10RB
IFNBR IFNGR IL27R IL6R IL21R IL23R IL10R TGFBR IL12R IL4R IL15R IL2R CD28 TCR
NFAT
STAT1 STAT3 STAT4 STAT6 STAT5
IKB
IFNG IL21 IL23 IL10 TGFB IL4 IL2 IL17 NFKB
SMAD3
IRF1
RUNX3 TBET GATA3 RORGT FOXP3
FIG. 3. Th differentiation regulatory graph. The top nodes correspond to inputs (APC and external input cytokines),
while nodes placed at the bottom correspond to key transcription factors. Nodes considered for reduction are emphasized in
grey. Greenarrowsdenoteactivations,redbluntarcsdenoteinhibitionswhilethebluearcfromNFKBtoIL17denotesadual
interaction (see 27 for details).
performed a series of simulations to assess the effects otherexpressingIL2. Fromthearclabels, itfollowsthat
of heterogeneous environments on Th cell differentia- the selection between these two stable states depends on
tion. This led us to identify stable states correspond- the concurrent activation of RORGT and FOXP3.
ing to canonical Th1, Th2, Th17 and Treg subtypes, but Figure 5 (top) displays the HTG obtained when sim-
also to hybrid cell types co-expressing combinations of ulating a naive T-helper cell stimulated in mixed pro-
Th1, Th2, Treg and Th17 markers in an environment- Th2/pro-Th17 environment, i.e. in the presence of IL4,
dependent fashion. IL6, TGFB, and in the absence of IL2. The resulting
Here, we apply the HTG construction to this reduced HTG merely comprises 13 nodes, to be contrasted with
model in order to demonstrate how the dynamics can be the 1146 states of the corresponding STG. Furthermore,
compressed in a meaningful way, emphasizing the struc- the HTG structure emphasizes the progressive commit-
ture of the underlying STG, as well as crucial decision ment of cells when following pathways from the root
points along dynamical pathways. In this respect, we to the leaves (stable states). The states encompassed
have selected a limited number of simulations leading to by other nodes belong to two or more basins of attrac-
STGs of increasing complexity. tion. Note that the system can reach four stable states,
Figure 4 displays the HTGs obtained when simulating more precisely two pairs of activated versus anergic Th2
a naive T-helper cell stimulated by an antigen present- RORGT+ subtypes. Within each of these pairs, the sta-
ing cells in the presence of IL2 alone, or in the presence blestatesdifferbytheexpressionlevelofFOXP3. Thela-
of pro-Th1, Th2, Treg or Th17 cytokines. In all cases bels associated with the arcs clearly emphasize the tran-
but the last one, we obtain a unique stable state corre- sitions implementing differentiation decisions. As illus-
sponding to the expected cellular state (activated Th0, trated in Figure 5 (bottom), the use of priorities signif-
Th1, Th2, or Treg). In each of these HTGs, all other icantly decreases the size of the dynamics; selecting up-
states reachable from the initial conditions are grouped datesofILR2,NFATandanyoftheSTATfactorsagainst
together into a single irreversible transient component, other component updates led to an HTG of 5 nodes (en-
encompassing between 25 and 255 states. The label as- compassing31states)insteadof13nodes(encompassing
sociated with each arc denotes the ultimate elementary 1146states),wherethetwoanergiccellulartypesarethe
transitionsgoingouttheHTGnode. Incontrast, inpro- only reachable stable states.
Th17 conditions, the system can reach two different sta- A thorough discussion of the biological significance of
ble states expressing Th17 transcription factor RORGT, theseobservationswouldgobeyondthescopeofthisarti-
IL10, IL21 and IL23, one expressing also FOXP3, the cle. However, these examples demonstrate the compres-
8
APC + IL2 Pro Th1 (IFNG) Pro TH2 (IL4, IL6)
i#25 i#79 i#255
IL2-
IL2+ IFNG+ IL4+
Proliferation+ IL2- IL10+
IL21+
IL23+
ss-1001000000000210100000000100020100 ss-1011000000000211000000100110020100 ss-1001110000000210011110010101021100
Pro Th17 (TGFB, IL6)
i#66
Pro Treg (TGFB) RORGT+ FOXP3+
i#37 i#39 i#91
ITLG2F-b+ IIILLL212+01++ RORGT+ TRGOFRBG+T+ IIILLL212-01++
IL23+ IL23+
ss-1001000000001210000001001100020100 ss-1001010000001210101110000101020110 ss-1001010000001210001111001101020110
FIG. 4. HTG representation of asynchronous simulations of naive Th cells in simple polarizing environments.
TheseHTGscorrespondtothesimulationofTh0cellsinthepresenceofAPCsignalling+IL2alone,withINFG(proTh1),or
withIL4+IL6(proTh2),orwithTGFB(proTreg),orwithTGFBandIL6(proTh17),fromtoptobottomandlefttoright.
Allothernodesaresettozeroattheinitialstate. Inthenotationofthelogicalstablestates(prefixedby”ss-”),thenodeorder
considered starts with APC, followed by the external input cytokines IFNB, IL2, IL4, IL6, IL10, IL12, IL15, IL21, IL23, IL27,
TGFB, followed the cytokines produced by the Th cell considered IL2, IL4, IL10, IL21, IL23, TGFB, then the transcription
andsignaltransductionfactorsTBET,TBET,GATA3,FOXP3,NFAT,STAT1,STAT3,STAT4,STAT5,STAT6,followedby
the proliferation node, followed by RORGT and IL17. Arc labels indicate transitions (regulatory component updates) driving
the system out of an HTG node toward another one.
sion and clarification of asynchronous simulations that V. CONCLUSION AND PROSPECTS
can be achieved using the HTG representation.
Hierarchicaltransitiongraphs(HTGs)emphasizerele-
vanttransientandasymptoticdynamicalproperties. We
have defined a novel algorithm enabling the compaction
of state transition graphs and the generation of HTGs
on-the-fly. This approach has been implemented into a
Finally, Figure 6 displays the reachability analysis be- developmentversionofthesoftwareGINsim,availableas
tween cellular types through the use of model-checking. a pre-release30.
For this, we considered the environmental conditions de- We have applied this approach to a comprehensive
fined by specific valuations of the 13 input components model for T-helper cell differentiation. Although this
(see Figure 5 of the original study, i.e. in reference 27) model still needs to be further refined and tested, the
and the stable expression patterns (see Figure 6 in ref- analysispresentedhereclearlydemonstratestheassetsof
erence 27). These stable states correspond to cellular theHTGrepresentation,whichleadstosignificantgraph
subtypes, which are stable under specific environmental compression and clearly emphasizes the organization of
conditions. So,foreachinputcombination,weverifythe the state space into attractors and basins of attraction.
existence of a direct path between each possible pair of Interestingly, applying our algorithm for HTG con-
cellular types. More precisely, we check whether there struction onto a HTG produces a further compacted
is a fixed valuation of inputs such that there is a direct graph-based representation of the dynamics, where the
pathbetweentwocellularsubtypes,C1andC2(without nodes correspond to basins of attraction.
goingthroughothercellularsubtypes)andC2beingsta- Should a given dynamics be too large and complex to
ble. Using this approach, we could reproduce the results be effectively compacted using the HTG representation,
obtained at the cost of extensive simulations in the orig- wecanrelyoncomplementarymethodspresentedinthis
inal study (Figure 7 in reference 27). Three main groups manuscript. These methods aim at reducing the size of
aredefinedovertheThcellsubtypes(Th0,Th1andTh2, the search space, including the model reduction method
see Figure 8 in reference 27). We could also verify that that preserves key dynamical properties and the defini-
Th1 and Th2 subtypes can never switch back to a Th0 tion of transition priority classes, relying on biologically
one, even when inputs are allowed to vary freely. How- well-founded assumptions. Moreover, we have presented
ever, in such case, switches between all cellular subtypes model-checking techniques to analyse reachability prop-
are possible within each group. erties. Used jointly, these methods enable the dynamical
9
i#112(1)
RORGT+ IL2R-
IL2- STAT5+ FOXP3+
STAT5+
FOXP3+
STAT6+
i#56 i#24 i#143 i#54
STAT5+
IL2R- RORGT+ FOXP3+ IILL22-R- STAT5+
IL2- RORGT+
STAT5+ FOXP3+ i#595
STAT6+ i#35
i#112(2) IL2R-
IL2- IL4-
STAT5+ i#11 RORGT+ IIILLL421RR0++- RORGT+ IL2R- IL2- IIILLL122013+++
IL2R- IL21+ IL2RA+ TGFB+
IL2RA+ IL23+ GATA3+ GATA3+
GATA3+ RORGT+ RORGT+
ss-1000110000001010011110010101021110
ss-1000110000001010001111011101021110
Activated GATA3+ RORGT+
ss-1000110000001010000000010101011010 IL4+ IL10+ IL21+ IL23+ Activated GATA3+ RORGT+ FOXP3+
Anergic GATA3+ RORGT+ ss-1000110000001010000000011101011010 IL10+ IL21+ IL23+ TGFB+
Anergic GATA3+ RORGT+ FOXP3+
i#19
FOXP3+
RORGT+
i#7
i#3
RORGT+ IL2RA+
IL2RA+ GATA3+
GATA3+ RORGT+
ss-1000110000001010000000010101011010 ss-1000110000001010000000011101011010
Anergic GATA3+ RORGT+ Anergic GATA3+ RORGT+ FOXP3+
FIG. 5. Compressed representation (HTG) of the asynchronous simulation of Th0 in the presence of APC
signaling + IL4 + IL6 + TGFB (combination of pro Th2 and Th17 cytokines, in the absence of IL2). Four
stable states can be reached: two pairs of activated versus anergic Th2 RORGT+ cells, differing by the expression of FOXP3.
The bottom part shows the HTG obtained for the same initial conditions but using two asynchronous priority classes. In
this configuration, transitions involving IL2R, NFAT, or any of the STAT factor are selected against those involving any
other component. In contrast with the results obtained without prioritization, only two stable states can be reached, both
corresponding to anergic Th2 RORGT+ cells, which differ by the expression of FOXP3. The labels associated with the arcs
emphasize the crucial transitions underlying the choice of one or the other differentiation state.
analysis of logical models of unprecedented sizes. wellsuitedtothedynamicalanalysisofcomplexnetworks
ItisworthnotingthattheHTGstructurecouldbecon- endowed with differentiation properties (i.e. present-
sideredinthecontextofotherformalapproaches relying ing multiple alternative stable states, which can be all
on state transition graphs, including Petri nets (see e.g. reachedfromgiveninitialconditions),astheycapturethe
8andreferencestherein)andpiecewise-lineardifferential general organization of the corresponding STGs. Hence,
equation(PLDE)models3,17. Model-checkingtechniques based on the analysis of HTG structures, we should be
also apply to these models, once their dynamics can be abletoidentifythecircuitsatthecoreofcellcommitment
representedbyKripkestructures12. Ourmodelreduction andtherebyfocusonthegenesresponsibleforirreversible
could be applied to PLDE models, but its impact on the decisions.
dynamics still needs to be clarified.
HTG construction could be optimized and improved, Analternativestrategytoanalyzelargeregulatorynet-
e.g. using parallel algorithms. Although depth- works takes advantage of their modularity. Recently,
first search algorithms are known to be difficult to we have defined a compositional framework that relies
parallelize31, different methods have been proposed to on process algebra to incrementally compose, abstract
tackle this problem2. and minimize (using the safety equivalence) logical reg-
Further analysis relying on HTG structures should al- ulatory modules, enabling impressive reductions of the
lowtheassessmentoffinerproperties. Forinstance,some dynamics22. However, as proper methods to decompose
well-established rules (Thomas’ rules36) assert that dif- large networks into functional modules are still lacking,
ferentiation (resp. homeostasis) phenomena lean on the wehavefocussedonregulatorynetworksthatencompass
action of a positive (resp. negative) circuit in the reg- several identical modules connected by inter-cellular sig-
ulatory graphs. In practice, circuit functionality analy- naling. Inthisrespect,onecouldtakeadvantageofHTG
sis often points to combinations of intertwined circuits, structures to come up with relevant model decomposi-
which are difficult to analyze. HTGs appear particularly tion.