A Qualitative Framework Dedicated to Toxicology

Benjamin Miraglio

, Gilles Bernot

, Jean-Paul Comet

and Christine Risso-de Faverney

Universit

e C

ote d’Azur, CNRS, I3S, Sophia Antipolis, France

Universit

e C

ote d’Azur, CNRS, ECOMERS, Nice, France

Keywords:

Computational Toxicology, Discrete Dynamic Systems, Rule-Based Modelling, Temporal Logic.

Abstract:

Emerging constraints have led the toxicology community to complete the classical paradigm of toxicology

with the study of molecular events underlying the toxicity of a chemical substance. This evolution motivates

the emergence of new modelling approaches for toxicology. In this article, we introduce a qualitative rule-

based formalism dedicated to the domain of toxicology. This new formalism departs from other rule-based

formalisms such as BioChAM because it directly encodes possible alterations of equilibrium, instead of mak-

ing equilibriums emerge from the dynamics of the model. Using a simple example of the energy metabolism,

we show that this formalism is able to describe both the normal evolution of a biological system and its possible

toxic disruptions.

1 INTRODUCTION

Toxicology can be deﬁned as the study of adverse

effects caused by exogenous chemical substances to

biological systems. The classical paradigm assumes

that the more an organism is exposed (in dose and/or

time) to a compound, the greater the compound ef-

fects will be. In these conditions, any chemical sub-

stance can therefore cause harmful effects to an or-

ganism if this organism is exposed during a long

enough time to a high enough dose of chemical.

That concept serves as a basis to the dose-response

relationship, which enables toxicologists to establish

a causality between the exposure to a chemical and its

induced observed effects. It also allows toxicologists

to determine the threshold of toxicity, namely the low-

est exposure (in dose and/or time) where an induced

effect is observed.

Many experiments carried out recently have

pointed out the limitations of this paradigm. Indeed,

toxicity assessment is quite complex and, besides the

dose and time of exposure, a lot of other factors can

affect the results of toxicity tests. In particular, these

factors include temperature, food, light, the route of

exposure and the chemical interactions of the tested

substance with other chemical compounds. Other fac-

tors related to the test subject itself, including age,

sex, genetics, health status, hormonal status or win-

dow of exposure may also greatly inﬂuence the vul-

nerability of an organism to a chemical substance.

To answer these limitations, an increasing trend in

toxicology is to focus on the causal sequence of key

events occuring during the toxic response and lead-

ing to an observable effect. These sequences, called

pathways of toxicity, lay the basis of the mechanistic

toxicology and include events from molecular, cellu-

lar and even organ scales.

As mechanistic toxicology allows a better under-

standing of molecular mechanisms leading to adverse

effects, it can cope with many difﬁculties mentioned

earlier, such as the extrapolation of toxicity ﬁndings

obtained from laboratory animals to humans or the

consideration of additional factors in toxicity assess-

ments. Moreover, as distinct pathways of toxicity can

share the same key events, knowledge obtained when

studying one chemical could be reused when assess-

ing other chemicals.

Concurrently, as the potential toxicity of chemical

exposure became an area of great concern to both the

public (Colborn et al., 1996; Kepner, 2004) and the

regulatory authorities (Backhaus et al., 2010), the pro-

duction of chemical compounds is increasingly reg-

ulated worldwide. Manufacturers must now conduct

extensive studies to demonstrate the innocuity of their

products, considerably increasing the cost of develop-

ment of such products.

This context favours the emergence of different

mathematical modelling approaches, and so far, most

of these approaches are quantitative and aim at either

inferring the toxic threshold of a chemical substance

Miraglio B., Bernot G., Comet J. and Risso-de Faverney C.

A Qualitative Framework Dedicated to Toxicology.

DOI: 10.5220/0006168200930103

In Proceedings of the 10th International Joint Conference on Biomedical Engineering Systems and Technologies (BIOSTEC 2017), pages 93-103

ISBN: 978-989-758-214-1

or conﬁrming its speciﬁc pathway of toxicity. These

objectives require a lot of toxicological data. This

can be restrictive given the current acquisition cost of

new biological data. There is thus an incentive to de-

velop methods that do not focus on toxic thresholds

but instead, aim at describing pathways of toxicity

in a qualitative manner. Such an approach would fo-

cus on equilibrium shifts and would therefore require

comparatively less toxicological data.

In this article, we present a new qualitative formal-

ism allowing to enumerate all the conceivable path-

ways of toxicity linked to a compound present in a

given biological system. As a lot of these modelled

pathways are biologically improbable, it is then pos-

sible to encode into temporal logic basic toxicolog-

ical knowledge and ﬁlter out the less relevant ones.

The remaining modelled pathways can ﬁnally serve

as a basis to design more informative experiments and

help toxicologists in their search of new pathways of

toxicity.

The next section of this article is dedicated to the

brief description of related work focusing on formal

frameworks dedicated to model biological systems.

As our formalism is presented alongside examples in-

spired from the energy metabolism, section 3 sketches

an overview of the energy metabolism and its key

components. In section 4, we explain how to use the

new formalism to describe the equilibrium changes

of a system. In section 5, we show how to integrate

toxicological knowledge in the system using temporal

logic. Finally, this formalism is applied to a simpliﬁed

model of the energy metabolism in section 6.

2 RELATED WORK

The rule-based formalism presented throughout this

article was originally inspired from the Boolean se-

mantics of BioChAM (Calzone et al., 2006), an en-

vironment able to model biological systems as net-

works of chemical reactions. However, several speci-

ﬁcities of toxicology make this environment not op-

timal when handling toxicological models, motivat-

ing the development of a more domain-oriented for-

malism. For instance, the notion of abnormal con-

centrations required us to adopt multivalued seman-

tics. Furthermore, the presence of modulations of

reactions, crucial in toxicology, is difﬁcult to han-

dle with BioChAM. This motivated us to directly de-

scribe equilibrium changes in a biological system, ab-

stracting any quantitative computing steps.

Pathway logic (Talcott, 2008) is another frame-

work dedicated to the description of biological sys-

tems while highlighting structural aspects of the cell.

Unfortunately, this framework is also not adapted to

deal with speciﬁcities of toxicology such as abnormal

concentrations or modulations of reactions.

Finally, Bio-PEPA (Ciocchetta and Hillston,

2009) is a modelling approach where biological sys-

tems are formalised as discrete models. These mod-

els include precise biological information such as ki-

netic laws and stoechiometry. Different analysis can

then be performed on these models, ranging from the

construction of continuous time Markov chains to the

translation of models into ordinary differential equa-

tions. However, the discretisation required for con-

tinous time Markov chains and probabilistic model

checking approaches is based on a precise knowl-

edge of kinetic parameters. Unfortunately, this pre-

cise knowledge is often unavailable, leaving room for

a framework where reaction rates are abstracted. Of

course, such a framework should maintain the possi-

bility for toxicologists to easily express concurrency

between equilibrium changes.

3 A QUICK INTRODUCTION TO

METABOLISM

The metabolism can be described as the set of chem-

ical reactions allowing cells to survive. A major dis-

tinction in metabolism is made between catabolism

and anabolism. Catabolism refers to the set of reac-

tions degrading molecules in smaller parts. Catabolic

reactions are mainly oxidations and tend to produce

energy for the cell. On the contrary, anabolism gath-

ers synthesis reactions producing key molecules for

the cell and tend to require energy.

In addition to the catabolism/anabolism separa-

tion, metabolic reactions can also be clustered de-

pending on the type of molecules they handle. Hence,

carbohydrate metabolism refers to pathways manag-

ing glucose and other carbohydrates and the lipid

metabolism regroups reactions that both break down

and synthesise lipids. Obviously, these different

clusters of reactions are heavily interconnected, and

the most notable intermediary molecules are called

metabolic crossroads.

The formalism developed in this article will be il-

lustrated by a simple model of the energy metabolism

of a mammalian hepatic cell in aerobic conditions.

Here, the energy metabolism can be understood as a

set of reactions including pathways belonging to car-

bohydrate and lipid metabolisms as well as additional

chemical reactions occuring in mitochondria. This set

of reactions is summarised in Figure 1 and frequent

references to this ﬁgure will be made during the de-

scription of the reactions. Note that the subcellular lo-

BIOINFORMATICS 2017 - 8th International Conference on Bioinformatics Models, Methods and Algorithms

Figure 1: Representation of a simpliﬁed energy metabolism

model. Arrow captions correspond to rule identiﬁers devel-

oped in Section 6.

cation of the reactions (cytosol or mitochondria) is ab-

stracted in this ﬁgure. Moreover, reactions favoured

in low energy situations are depicted by using plain

arrows while reactions promoted in high energy cases

are shown with dotted arrows.

Carbohydrate Metabolism. Carbohydrates consti-

tute one of the main energy supply for cells and their

metabolism is essentially located in the cell cytosol.

The main catabolic pathway is glycolysis (GlyL),

which converts glucose (Gluc) and several other car-

bohydrates into pyruvate (Pyr) (Pilkis and Granner,

1992). Pyr then migrates to mitochondria where it

is decarboxylated (PyDC) into acetyl-CoA (ACoA),

an important metabolic crossroads. Note that the de-

carboxylation of Pyr into ACoA is performed by the

pyruvate dehydrogenase complex (PDH) (Pettit et al.,

1975). This enzymatic complex will not be included

in the application model developed in Section 6 and

is not represented in ﬁgure 1. However, it is used in

the running example illustrating Section 4.

The anabolic counterpart of glycolysis, the glu-

coneogenesis (GNG), uses energy to create new

molecules of Gluc (Pilkis and Granner, 1992). It

should be noted that Pyr carboxylation (PyC) into ox-

aloacetatic acid (OAA) - the ﬁrst step of gluconeo-

genesis - directly competes with Pyr decarboxylation

(PyDC) to ACoA. The competition outcome depends

on the quantities of available ACoA: if there is a short-

age, Pyr is decarboxylated to ACoA, otherwise, it is

carboxylated in OAA (Owen et al., 2002). The other

reactions composing gluconeogenesis are mainly re-

versed glycolysis reactions.

Lipid Metabolism. Lipids are another important

source of energy for cells. The β-oxidation (βOx) is

the main catabolic pathway and degrades fatty acids

(FA) into ACoA in mitochondria (Schulz, 1991).

Conversely, lipogenesis (LipoG) occurs in the cytosol

and converts ACoA into FA (Hellerstein et al., 1991).

It is interesting to note that both lipid and car-

bohydrate catabolic pathways result in the produc-

tion of ACoA. Indeed, ACoA is actually an important

metabolic crossroads and can be involved in the pro-

duction of either new molecules - notably lipids - or

energy, through its inclusion in the citrate cycle (CC).

The inclusion of ACoA in the citrate cycle results in

the destruction of ACoA, represented in ﬁgure 1 by

the symbol ”∅”.

Citrate Cycle. The citrate cycle (CC), also known

as tricarboxylic acid cycle, is a key component of en-

ergy production (Owen et al., 2002). It starts with

the conjugation of ACoA and OAA to form citrate

(Cit) and is then composed of a sequence of reac-

tions resulting in the regeneration of OAA and the

production of energy (Owen et al., 2002). Hence, the

overall cycle only decreases ACoA levels to produce

energy, therefore, only ACoA is consumed in Figure

1. On top of the production of energy, the citrate cy-

cle also produce reduced compounds. These reduced

compounds are then used to fuel the mitochondrial

respiratory chain (MRC, not represented in Figure 1),

generating even more energy (Chance and Williams,

1956).

Reduction-oxidation Reactions. The reduced

compounds are actually coenzymes such as nicoti-

namide adenine dinucleotide or ﬂavine adenine

dinucleotide that act as electron transporters in the

cell and can have either an oxidised or a reduced

form. For the sake of simplicity, we will summarise

the ratio between reduced and oxidised coenzymes

as the reducing potential of a cell (pRed). A high

(resp. low) reducing potential therefore means high

(resp. low) concentrations in reduced coenzymes

and subsequent low (resp. high) levels of oxidised

coenzymes.

As a matter of fact, the whole metabolism can be

seen as a network of reduction-oxidation (redox) re-

actions and the reducing power of a cell can inform

on the state of the cell metabolism. For instance, both

carbohydrate and lipid catabolic pathways reduce oxi-

dised coenzymes in their reduced form (increasing the

cell reducing potential), while conversely, both an-

abolic pathways oxidise reduced coenzymes in their

oxidised form (decreasing the reducing potential).

A Qualitative Framework Dedicated to Toxicology

For the sake of clarity, modiﬁcations of the reduc-

ing potential are not directly shown in Figure 1 but

can be deduced from the arrows nature. Indeed, re-

actions favoured in low energy states tend to produce

reducing potential (solid arrows in ﬁgure 1) while re-

actions favoured in high energy states tend to decrease

the reducing potential of a cell (dotted arrows in ﬁgure

1).

As previously said, mitochondria are able to create

energy from the reduced coenzymes. This process de-

creases the reducing power, regenerating the oxidised

coenzymes supplies of the cell, allowing the continu-

ation of catabolic reactions.

Another important source of coenzymes oxidation

lays in the cytosolic reduction of Pyr into lactate (Lac)

(Brooks, 1998). If mitochondria functioning falters,

the pyruvate reduction (PyRed) is still able to resup-

ply the cell in oxidised coenzymes, preventing the in-

terruption of catabolic pathways. Since pyruvate re-

duction is reversible, its reversed reaction (lactate de-

hydrogenation or LDH), can also generate reducing

potential while converting Lac in Pyr. The newly cre-

ated Pyr can, in turn, either generate a lot of energy

through its decarboxylation in ACoA or participate to

the gluconeogenesis through carboxylation in OAA

(see Carbohydrate metabolism paragraph).

Disruptions of the Cell Reducing Potential. As

the cell reducing potential is central to a lot of

metabolic reactions, disruptions of the redox equilib-

rium can have great impacts on the organism. For

instance, alcohol dehydrogenase is able to convert

ethanol (Eth, not represented in Figure 1) to acetalde-

hyde while greatly increasing the cell reducing poten-

tial. If this increase is strong enough, it can saturate

the mitochondrial respiratory chain and trigger the re-

duction of Pyr into Lac. The subsequent accumula-

tion of Lac in cells can then lead to troubles such as

metabolic acidosis.

Furthermore, as the presence of Eth triggers Pyr

reduction into Lac, the amount of Pyr available to car-

boxylation in OAA is decreased. The subsequent pro-

duction of Gluc from OAA through the gluconeogen-

esis is then strongly impacted. Since gluconeogenesis

is an important pathway to address hypoglycemia, a

fasting organism can thus see its capacity to recover

from hypoglycemia strongly damaged after its expo-

sition to ethanol (Field et al., 1963).

4 DESCRIBING EQUILIBRIUM

CHANGES

A biological system can be described as a set of bio-

logical entities interacting with each other at different

concentrations. In a given organism, each entity has

a concentration regarded as normal in standard condi-

tions. For instance, the normal blood concentration of

glucose is about 1 g/L in an adult human.

Our domain-oriented formalism allow us to rep-

resent the evolution of the concentration of each en-

tity and to depict abnormal concentrations from which

toxicity can arise. Indeed, we introduce four qualita-

tive abstract levels, which are listed here in increasing

order:

• ε reﬂects a negligible concentration of a given en-

tity, that is to say a concentration too low to trigger

any reaction in the biological system.

• ι conveys an abnormally low concentration, i.e.

a relative lack of this entity that can affect some

mechanisms in the biological system.

• ∆ indicates a normal concentration.

• θ shows an abnormally high concentration,

namely an excess of this entity.

Notation 1. [Concentration levels] We note L the set

{ε, ι,∆, θ} equipped with the total order relation such

that: ε < ι < ∆ < θ. The elements of L are called

concentration levels.

In a given biological system and depending on the

studied issue, not all entities have concentrations re-

garded as abnormally low or high. Therefore, only

the levels ε and ∆ are mandatory for each entity, ι and

θ being optional.

Taking this variation in consideration, the signa-

ture of a biological system allows the deﬁnition of the

set of biological entities considered in the system and,

for each entity, its admissible concentration levels.

Deﬁnition 1. [Signature] A signature is an applica-

tion E : E → P (L) where E is a ﬁnite set and for all

e ∈ E, {ε, ∆} ⊂ E(e). Elements of E are called en-

tities and for each entity e, E (e) is called the set of

admissible levels of e.

For instance, the signature of a ba-

sic energy metabolism model may involve

E = {PDH, Pyr, ACoA, OAA, Cit} can correspond to

the set of ﬁve entities where each entity has its own

set of admissible levels. For example, we may have

E(Pyr) = {ε, ι, ∆, θ}.

After deﬁning the system signature, the state of

the system can be deﬁned as the qualitative level

of each entity present in the system. The previous

example model can be at a state η

where PDH

BIOINFORMATICS 2017 - 8th International Conference on Bioinformatics Models, Methods and Algorithms

is at the level ∆, noted η

(PDH) = ∆ and where

(Pyr) = θ, η

(ACoA) = ι, η

(OAA) = ∆ and

(Cit) = ε. This state can also be written:

= (∆, θ, ι, ∆, ε) (1)

where the order of variable is (PDH, Pyr, ACoA,

OAA, Cit).

Deﬁnition 2. [State] A signature E being given, the

set of states ζ is the set of functions η : E → L such

that for all e ∈ E, η(e) ∈ E (e).

In this formalism, the evolution of the system

is represented by two functions: the incrementa-

tion, noted incr, and the decrementation, noted decr.

These functions apply to one entity at a time and

return the level of this entity just above (resp. be-

low) its current level. Because all entities do not

have the same set of admissible levels, there is

one function deﬁned for each entity. For instance,

if E (OAA) = {ε, ∆, θ} and η

(OAA) = ∆, then

incr

OAA

(η

(OAA)) = θ and decr

OAA

(η

(OAA)) =

ε. It should be noted that the incrementation

(resp. decrementation) function is not deﬁned on

the maximal (resp. minimal) level of the admis-

sible levels. Therefore, in our previous example,

incr

OAA

(η(OAA)) is not deﬁned if η(OAA) = θ.

Besides these functions, the formalism also makes

use of formulas to describe properties about the enti-

ties concentration levels.

Deﬁnition 3. [Formula] The set A of atomic formu-

las on a signature E is the set of expressions of the

form a 6 b where a and b can be any element of E ∪L.

The set F of formulas on a signature E is induc-

tively deﬁned by:

• A ⊂ F .

• if ϕ and ψ are elements of F , then ¬ϕ, ϕ ∧ ψ,

ϕ ∨ ψ, ϕ ⇒ ψ are also elements of F .

Deﬁnition 4. [Satisfaction relation] A state η and

a formula ϕ ∈ F on a signature E being given, the

satisfaction relation η  ϕ is inductively deﬁned by:

• if ϕ is an atom of the form a 6 b, then η  ϕ if and

only if η(a) 6 η(b) where η is the extension of η

to E ∪ L by the identity on L.

• if ϕ is of the form ϕ

∧ϕ

then η  (ϕ

∧ϕ

) if and

only if η  ϕ

and η  ϕ

. We proceed similarly

for the other connectives.

Moreover, ”η  ϕ” is read ”η satisﬁes ϕ”.

We may also use the abbreviation a = b as a short-

cut for (a 6 b) ∧ (b 6 a). We proceed similarly for

a < b, a > b and a > b.

Examples of formulas can be ϕ ≡ (Pyr = θ),

stating an excessive presence of Pyr or ψ ≡ (Cit >

ACoA), stating that the qualitative level of Cit is

strictly superior to the one of ACoA. The state η

previously described in eq. 1, satisﬁes ϕ but not ψ.

To determine the evolution of the system, a set of

rules is then used. A rule can be interpreted as possi-

ble modiﬁcations in the state which can be abstracted

by the following representation:

r : A

+ ··· + A

⇒ A

m+1

+ ··· + A

boost(ϕ) block(ψ)

Beside its identiﬁer r, each rule includes two sets

of entities. The ﬁrst one, for all i in [1, m], constitutes

the set of “reactants” whose level can be reduced by

the application of the rule. The other one, for all i in

[m+1, n], represents the set of “products” whose level

can be increased by the application of the rule. A rule

also includes two modulating conditions boost(ϕ) and

block(ψ) (ϕ and ψ being formulas) representing re-

spectively a possible positive and negative modula-

tion of the rule. The boost(ϕ) (resp. block(ψ)) mod-

ulation takes only effect if ϕ (resp. ψ) is satisﬁed and

its effects are further detailed later on. Of course, if

no modulation is known for a given rule, boost and

block regulations are not displayed in the rule repre-

sentation.

Deﬁnition 5. [Biological action network] A biologi-

cal action network on a signature E, or E -action net-

work, is a set R of rules of the form:

(1) r : A

+ ··· + A

⇒ A

m+1

+ ··· + A

boost(ϕ) block(ψ)

where:

• r is an identiﬁer such that there are not two rules

in N with the same r.

• ∀i = 1 . . . n, A

∈ E.

• {A

. . . A

} ∩ {A

m+1

. . . A

} = ∅.

• ϕ and ψ are elements of F .

For short, we will call such rules E-rules and we

will call state of R a state on the signature of R.

Let us note that a rule representing possible alter-

ations of a biological state, it makes no sense to have

an entity being part of both reactants and products of

a same rule.

Moreover, please notice that a rule can be devoid

of any reactant or product. In the previous deﬁnition,

the index m can be equal to zero (the rule does not

need any reactant) or m can be equal to n (the rule

has no product). A rule without reactant can be con-

sidered as the constitutive production of an entity in a

given model and a rule without product can be inter-

preted as the degradation of an entity. In either cases,

the empty solution is depicted using the symbol.

It is worth mentioning that despite the strong re-

semblance between a rule and a chemical reaction,

a rule must not be interpreted as quanta of reactants

converted into quanta of products but as a possible

evolution of the levels of entities present in the rule.

A Qualitative Framework Dedicated to Toxicology

As a basic example of rule, the condensation of

acetyl-CoA and oxaloacetate to form Cit can be rep-

resented by the following rule:

: ACoA + OAA ⇒ Cit

Since neither positive nor negative modulating

conditions are considered here, only reactants and

products are displayed.

In order to be applicable at a given state, a rule

must meet basic criteria inspired from biology. First,

since the level ε is interpreted as a negligible concen-

tration, a rule is applicable only if all its reactants are

present at least at the level ι. In addition, a rule can-

not be applied if the negative modulating condition

block() applies, namely if the corresponding formula

is satisﬁed.

Deﬁnition 6. [Applicable rule] Let us consider a

state η on a signature E. An E -rule r ∈ R of the form

(1) is said applicable at the state η if and only if:

• ∀i = 1 . . . m, η(A

) 6= ε.

• η 2 ψ.

For instance, let us consider the conversion of Pyr

into acetyl-CoA by the enzyme name Pyr dehydroge-

nase (PDH). If we assume that E(PDH) = {ε, ι, ∆, θ},

the conversion can be written as:

: Pyr ⇒ OAA block(PDH < ∆)

This rule is applicable if and only if the level of

Pyr is strictly greater than ε and the level of PDH is

at least ∆, namely if there is Pyr in the system and a

normal concentration of PDH. Note that the catalysis,

namely the necessary presence of an enzyme to the

proper conduct of a reaction, can be expressed using

the block() condition as in the previous example.

Considering rules as possible alterations of the

equilibrium led us to a difﬁcult choice on the way of

handling the emergence of abnormal levels of prod-

ucts.

On one hand, we can consider that every reactant

of a rule has to be in excess to shift the product(s)

equilibrium(s) to excessive levels. This vision is la-

belled optimistic since we suppose the system less

prone to drift towards excessive states. On the other

hand, with a more pessimistic perspective, only one

reactant in excess is enough to propagate the excess

to the product(s).

Both the optimistic and the pessimistic approaches

need an exception system to take into account par-

ticular biological cases. In the optimistic approach,

the exception is a boost() statement that relaxes the

conditions for a product to reach excessive levels. In

the pessimistic approach, the exception would be a

brake() statement selectively preventing excesses to

propagate.

As the optimistic approach is predominant in tox-

icological data, our formalism implements it thanks

to stricter conditions for a product level to increase.

These constraints can be found in the way of comput-

ing the potential next level of products.

Deﬁnition 7. [Potential next level] Let R be an E -

action network, let η be a state of R, let r ∈ R. We

note η



: E → L the partial function such that η



(e)

is deﬁned if and only if r is applicable and if one of

the following conditions is satisﬁed:

• e ∈ {A

. . . A

} and in this case η



(e) =

decr

(η(e)).

• e ∈ {A

m+1

. . . A

}, η(e) < max(E(e)), and in this

case:

– if η 2 ϕ and η(e) < min

i∈{1...m}

(η(A

)) then



(e) = incr

(η(e)).

– if η  ϕ then η



(e) = incr

(η(e)).

The potential next level of an entity through an

applicable rule refers to the next level of the entity

after the application of that rule. If the entity acts as a

reactant, its potential next level is the one returned by

the decrementation function applied to that entity.

If the entity acts as a product, its potential next

level depends on the boost() statement:

• if the boost() statement is not satisﬁed, a product

level can increase only if all the reactants levels

are strictly greater (this is due to the optimistic

vision explained previously). In this case, the po-

tential next level of a product is thus the one re-

turned by the incrementation function applied to

the product.

• if the boost() statement is satisﬁed, the previous

restriction no longer applies. In such cases, the

potential next level of a product is returned by

the incrementation function applied to it, indepen-

dently of the reactant levels.

Let us note that the potential next level is returned

either by the incrementation or decrementation func-

tion. Therefore, when these functions are not deﬁned,

the potential next level of an entity is also not deﬁned.

Keeping the conversion of Pyr as an example, we

can also specify that an excess of Pyr dehydrogenase

can cause trouble in oxaloacetate levels by adding a

boost() condition to the rule r

: Pyr ⇒ OAA block(PDH < ∆) boost(PDH > ∆)

Here, assuming that the rule is applicable at the

state η

and that η

(OAA) = ∆, the potential next

level of oxaloacetate by this rule can be θ only if

(Pyr) = θ or if η

(PDH) > ∆ (so, η

(PDH) = θ).

Among all the applicable rules at a given state,

only one is applied at a time. When a rule is applied,

one and only one of its entities sees its level chang-

ing to its potential next level. This means that the

BIOINFORMATICS 2017 - 8th International Conference on Bioinformatics Models, Methods and Algorithms

level of an entity has to change in order to consider

that the rule was applied. Importantly, this also means

that a product cannot be updated simultaneously with

a reactant, and conversely. Similar ideas have been

ﬁrstly developed for discrete gene models by Thomas

and Snoussi (Snoussi, 1989; Thomas, 1991). This be-

haviour reﬂects the possibility for an entity to cross

a threshold without all the other entities levels doing

likewise.

In brief, starting from a given state, it is possible

to determine which rules of the system are applicable

at that state. The application of one of these rules then

changes the level of one entity, modifying the system

state. However, it is possible to stay indeﬁnitely at a

same system state thanks to the rule Id (whose appli-

cation does not change the levels of the system entities

and that is always applicable).

It is then possible to establish a transition graph,

mapping all the possible transitions between the states

of a system. An inﬁnite succession of transitions such

that the output state of a transition is the input state

of the next one is here called a path of the transition

graph.

Deﬁnition 8. [Transition graph] The transition

graph of an E -action network R is the labelled graph

whose set of vertices is the set of states ζ and the set

of edges T is the set of transitions of the form η

−→ η

such that one of the following condition is satisﬁed:

• r = Id and η

= η

• r ∈ R and there exists an entity e ∈ E such that



(e) is deﬁned and:

– η

(e) = η



(e)

– ∀ e

∈ E r {e}, η

) = η(e

Remark: The transition graph of an E-action net-

work R canonically deﬁnes a labelled Kripke struc-

ture L = (L, Σ, T ) as follows:

• L(η) = {α ∈ A | η  α}.

• Σ = R ∪ {Id}.

• T can obviously be seen as the set of triplets

(η, r, η

) such that (η

−→ η

) is a transition of T .

A path (π ≡ η

−→ η

−→ . . .

i−1

−−→ η

−→ . . . ) is then

an inﬁnite sequence of labelled transitions such that

the input state of r

is equal to the output state of r

i−1

for all i > 0. The set of paths is called Π

5 INTEGRATING

TOXICOLOGICAL

KNOWLEDGE

As the transition graph of a biological system includes

many biologically improbable paths, it is necessary

to ﬁlter out the irrelevant ones and only characterise

the interesting paths for toxicologists. Temporal logic

and model checking tools have been successfully ap-

plied to biological systems, either using Linear Tem-

poral Logic (Ito et al., 2014) or Computation Tree

Logic (Bernot et al., 2004). Here, since we seek to

ﬁlter paths, we need a logic able to express both state

and transition properties. We thus use the state/event

linear temporal logic (SE-LTL) developed by Chaki

(Chaki et al., 2004).

Since a path can be seen as an inﬁnite alternance

between states and transitions, atomic temporal for-

mulas concern either a state or a transition. For states,

atomic temporal formulas are similar to atomic for-

mulas exposed in Deﬁnition 3. For transitions, atomic

temporal formulas only involve a rule identiﬁer or the

identity operator.

Deﬁnition 9. [Temporal formula] Given an E-

action network R, the set T

of temporal formulas on

R is inductively deﬁned by:

• (A ∪ R ∪ {Id}) ⊂ T

• if ϕ and ψ are formulas of T

, then ¬ϕ, ϕ ∧ ψ,

ϕ ∨ ψ, ϕ ⇒ ψ, Xϕ, Fϕ, Gϕ, ϕUψ are formulas of

Deﬁnition 10. [Temporal formula satisfaction]

Given an E -action network R and a path (π ≡ η

−→

−→ . . . ) ∈ Π

, the satisfaction relation ⊂ Π

×T

is inductively deﬁned on the temporal formulas of T

by :

• π  α where α ∈ A if and only if η

 α,

• π  r where r ∈ R ∪ {Id} if and only if r = r

• π  ϕ ∧ ψ where (ϕ, ψ) ∈ T

if and only if π  ϕ

and π  ψ, other propositional logic connectives

are treated similarly,

• π  Xϕ where ϕ ∈ T

if and only if (η

−→ η

−→

. . . )  ϕ,

• π  Gϕ where ϕ ∈ T

if and only if for all i ∈ N,

(η

−→ η

i+1

−−→ . . . )  ϕ,

• π  Fϕ where ϕ ∈ T

if and only if there exists

i ∈ N, (η

−→ η

i+1

−−→ . . . )  ϕ,

• π  ϕ U ψ where (ϕ, ψ) ∈ T

if and only if there

exists j ∈ N, (η

−→ . . . )  ψ and for all 0 6 i < j,

(η

−→ . . . )  ϕ.

A Qualitative Framework Dedicated to Toxicology

Furthermore, for all r ∈ R of the form r : A

+···+

⇒ A

m+1

+ ··· + A

boost(ϕ) block(ψ), we note

app(r) the temporal formula (

i=1

> ε) ∧ ¬ψ stat-

ing that r is applicable at the current state (see Deﬁni-

tion 6).

In addition, for all e ∈ E, we note ↓ e the temporal

formula stating that the level of the entity e decreases

in the next state:

l ∈ E (e)r{ε}



e = l ∧ X



e = decr

(l)





We proceed similarly for ↑ e.

For instance in our running example, the property

χ characterising paths where an excess of Pyr leads

to a future excess of oxaloacetate can be written as:

G(Pyr > ∆ ⇒ F(OAA > ∆)) and the formula ξ stating

that the rule r

is the ﬁrst applied when Cit is absent

from the system can be written as: G(Cit = ε ⇒ r

In this situation, the path beginning with (η

−→ η

where η

= (∆, θ, ι, ∆, ε) and η

= (∆, θ, ι, θ, ε) satis-

ﬁes χ but not ξ.

Finally, the association of the transition graph of a

system with a set of properties representing the rele-

vant biological pathways is called a constrained net-

work. This constrained network is actually a subset of

paths from the transition graph, with each path in this

subset satisfying all the expressed biological proper-

ties.

Deﬁnition 11. [Constrained network] An E-

constrained network is a couple N = (R, Ax) where

R is an E -action network and Ax is a set of temporal

formulas.

Deﬁnition 12. [Dynamics of a constrained net-

work] Given an E -constrained network N = (R, Ax),

the dynamics of N is the subset Π

of Π

such that

π ∈ Π

belongs to Π

if and only if π  Ax.

Since properties ﬁlter out irrelevant paths from the

transition graph, it is thus possible to use them in con-

junction to formal methods to insure that the ﬁnal con-

strained network respects basic biological and toxico-

logical properties as well as speciﬁc properties related

to the studied issue.

6 APPLICATION TO THE

ENERGY METABOLISM

To illustrate the formalism previously described, we

continue with the simpliﬁed energy metabolism de-

veloped in Section 3. The set of elements introduced

in Figure 1 is thus completed with the cell reducing

potential (pRed) and ethanol (Eth). The resulting sys-

tem signature is E

= {Gluc, Lac, Pyr, ACoA, OAA,

FA, pRed, Eth}.

In order to maximise the amount of possible paths,

let us consider that the set of admissible levels of the

endogenous entities (namely all entities except Eth)

is {ε, ι, ∆, θ}. In parallel, we will consider the ethanol

either absent, present moderately or present in excess

in the cell. Its set of admissible levels is thus {ε, ∆, θ}.

The role of each of these entities was developed

in Section 3. The reducing potential is here desig-

nated as pRed, but as previously described, it is an

abstraction of the balance between oxidised and re-

duced coenzymes. The rule ⇒ pRed thus abstracts

the rule OxidisedCoenzymes ⇒ ReducedCoenzymes.

As the amounts of oxidised and reduced coen-

zymes are interdependent, an increase in pRed means

both an increase in reduced coenzymes and a decrease

in oxidised coenzymes. This means that the rule

⇒ pRed cannot apply when there is an important

lack of oxidised coenzymes, namely when the reduc-

ing potential is in excess. This explains the presence

of block modulations linked to an excess of pRed in

some of the rules presented hereunder (GlyL, PyDC,

LDH, CC, βOx and EthOx).

The following rules summarise the interactions

described in Section 3 and constitute the E

-action

network R

GlyL : Gluc ⇒ Pyr + pRed block(pRed = θ)

PyDC : Pyr ⇒ ACoA+pRed block(pRed = θ∨ACoA > ∆)

PyC : Pyr ⇒ OAA

PyRed : Pyr + pRed ⇒ Lac

LDH : Lac ⇒ Pyr + pRed block(pRed = θ)

CC : ACoA ⇒ pRed block(pRed = θ ∨ OAA = ε)

MRC : pRed ⇒

GNG : OAA + pRed ⇒ Gluc

LipoG : ACoA + pRed ⇒ FA

βOx : FA ⇒ ACoA + OxPw block(pRed = θ)

EthOx : Eth ⇒ pRed block(pRed = θ) boost(Eth > ε)

Rule GlyL abstracts the whole glycolysis with its

transformation of Gluc into Pyr concomitant to the

production of pRed. The different futures of Pyr are

summarised in rules PyDC, PyC and PyRed. It should

be noted that the block modulation of PyDC repre-

sents possible conditions in which Pyr decarboxy-

lation is stopped while LDH is the reverse of rule

PyRed.

The rule CC represents the citrate cycle and

its production of pRed through the consumption of

ACoA. As OAA is an integral part of the cycle, this

rule is stopped by a lack of OAA. The rule MRC rep-

resent the ability for mitochondria to regenerate the

oxidised coenzyme (and therefore to decrease pRed).

Rules GNG and LipoG represent both carbohy-

drate (gluconeogenesis) and lipid (lipogenesis) an-

abolic pathways and their use of pRed. Conversely,

rule βOx abstracts β-oxidation and its production of

BIOINFORMATICS 2017 - 8th International Conference on Bioinformatics Models, Methods and Algorithms

100

Table 1: The set of states present in path π

with rules ap-

plied at each step of the path. In all these states, FA and

pRed levels are normal (∆) while there is a lack of Lac (ι)

and a complete absence of Eth (ε). The level updated be-

tween each state is shown in bold.

State Gluc Pyr OAA ACoA

∆ ι ι ι

∆ ∆ ι ι

ι ∆ ι ι

ι ∆ ∆ ι

∆ ∆ ∆ ι

pRed.

Finally, the detoxiﬁcation of Eth by the cell is de-

picted in EthOx, with the boost modulation represent-

ing the important amount of pRed possibly produced

by the detoxiﬁcation. It should be noted that this

detoxiﬁcation is known to be performed with high pri-

ority by the hepatic cell, but this kind of information

cannot appear in the rule.

Instead, it will be integrated in the model thanks to

temporal formulas as seen in Section 5. For instance,

it is known that Pyr can either be carboxylated (PyC)

or decarboxylated (PyDC), depending on the amount

of available ACoA. This can be summarized in the

temporal formula ϕ

G((ACoA < ∆ ∧ app(PyC) ∧ app(PyDC)) ⇒ ¬PyC)

This property states that when PyC and PyDC are

applicable, and that there is a lack of ACoA, the rule

PyC is not applied. Note that the G operator surround-

ing the formula indicates that the property remains

true at every step of the path. PyC thus never applies

when the previous conditions are satisﬁed.

It is also known that both decarboxylation and car-

boxylation of pyruvate prevail over pyruvate reduc-

tion (PyRed), as written in ϕ

G((app(PyRed) ∧ (app(PyC) ∨ app(PyDC))) ⇒ ¬PyRed)

This property is similar to the previous one and

litterally means that when PyRed and either PyC or

PyDC are applicable, PyRed never applies.

As previously said, these properties can be used to

characterise interesting paths allowed by R

. Let then

be the constrained network associating R

with the

two previous properties (ϕ

and ϕ

), and let us con-

sider the path π

beginning with the following preﬁx

(see Table 1):

GlyL

−−−→ η

GlyL

−−−→ η

PyC

−−→ η

GNG

−−−→ η

PyC

−−→ . . .

This preﬁx starts in state η

. This state can be as-

similated to a fasting state, where there is a normal

amount of energy supplies (Gluc, FA) and a lack of

every other metabolite. The ﬁrst applied rule, GlyL

is the glycolysis, restablishing the normal level of Pyr

(η

). The rule is then triggered again, leading to a de-

Table 2: The set of states present in path π

with rules ap-

plied at each step of the path. In all these states, Lac, ACoA

and FA levels are normal (∆). The level updated between

each state is shown in bold.

State Gluc Pyr OAA pRed Eth

ι ∆ ι ∆ ε

ι ∆ ∆ ∆ ε

∆ ∆ ∆ ∆ ε

Table 3: The set of states present in path π

with rules ap-

plied at each step of the path. In all these states, Lac, ACoA

and FA levels are normal (∆). The level updated between

each state is shown in bold.

State Gluc Pyr OAA pRed Eth

ι ∆ ι ∆ θ

ι ∆ ι ∆ ∆

ι ∆ ι θ ∆

ι ι ι θ ∆

ι ι ι ∆ ∆

crease in Gluc level. The rule PyDC is then applied,

leading to an increase in OAA (η

). Finally, OAA fu-

els the gluconeogenesis (GNG), restablishing normal

glucose levels (η

Although the preﬁx of π

is allowed by the E

action network, it does not satisfy ϕ

. Indeed, the

priming of Pyr carboxylation (PyC) on Pyr decar-

boxylation (PyDC), namely the change from η

, in a fasting situation is very unlikely. As the pre-

ﬁx does not satisfy one of the properties of the con-

strained network, the entire path is ﬁltered out.

Properties are also important to express toxicolog-

ical knowledge. For instance, the priming of ethanol

detoxiﬁcation (EthOx) over all the other rules can be

expressed by ψ

G(app(EthOx) ⇒ EthOx)

Furthermore, the detoxiﬁcation reaction produce

a lot of reducing potential. This means that whenever

an excess of ethanol is detoxiﬁed through EthOx, the

reducing potential of the cell is impacted, hence the

formula ψ

G((Eth = θ ∧ EthOx) ⇒ (↓ Eth ∧ X(EthOx∧ ↑ pRed)))

Finally, if the cell struggles to decrease the pRed

level in presence of Eth, this means that the mitochon-

drial respiratory chain is saturated. In such cases, the

reduction of Pyr (PyRed) is triggered, hence the for-

mula ψ

G((Eth > ε ∧ ¬app(EthOx) ∧ app(PyRed)) ⇒ PyRed)

To illustrate the consequences of such constraints,

let us consider N

, the constrained network associat-

ing R

with the set of properties {ψ

, ψ

}. Let us

also consider the path π

beginning with the follow-

A Qualitative Framework Dedicated to Toxicology

101

ing preﬁx (see Table 2):

PyC

−−→ η

GNG

−−−→ η

PyDC

−−−→ . . .

This path starts in the state η

, where there is no

ethanol, a lack of Gluc and OAA and a normal amount

of every other metabolite. This state corresponds to

an hypoglycemic state where the cell has the ability

to quickly recover its normal glucose level through

its supplies in FA and Pyr. Indeed, π

illustrates one

of the possible paths leading to the regeneration of

glucose, ﬁrst with the application of PyC to restore

OAA levels, then with the application of GNG.

Let us now consider π

, a path similar to π

except

on the beginning state (see Tables 2 and 3):

PyC

−−→ η

GNG

−−−→ η

PyDC

−−−→ . . .

This path starts from η

, which is identical to η

except for the excessive presence of ethanol (see Ta-

ble 3). This path does not satisfy ψ

since ethanol

detoxiﬁcation is not performed as soon as possible.

is thus not retained in N

Finally, let us take a look at path π

(see Table 3):

EthOx

−−−−→ η

EthOx

−−−−→ η

PyRed

−−−−→ η

PyRed

−−−−→ η

EthOx

−−−−→ . . .

Also starting from η

, π

then leads to the detox-

iﬁcation of a part of Eth through EthOx (η

), satis-

fying ψ

. The rule EthOx is then applied again, lead-

ing to the increase in pRed and satisfying ψ

. Rule

PyRed is then applied, decreasing the level of Pyr to

ι. As pRed is still in excess, PyRed is applied a second

time, leading to the recovery of normal levels of pRed

and satisfying ψ

. As π

satisﬁes ψ

, ψ

and ψ

, it

is retained in N

, illustrating the impaired ability for

a cell to regenerate normal glucose levels through the

gluconeogenesis in presence of ethanol.

7 CONCLUSION

In this article is presented a new formal framework

able to handle several speciﬁcities of the toxicology

domain not taken into account so far. This rule-based

modelling framework relies on the direct description

of equilibrium changes happening in a biological sys-

tem. This description does not model the difference

of reaction speed between the model rules, which can

affect the system equilibrium. It is however possi-

ble to integrate biological and toxicological knowl-

edge about rule kinetics through formulas expressed

in temporal logic.

As demonstrated on a simple model of the en-

ergy metabolism, its expressive power allows us to

describe both the equilibrium changes in the biologi-

cal system and knowledge about the prioritisation of

reactions. This knowledge is then used to ﬁlter out

irrelevant paths from the resulting model.

In the future, our formalism will be coupled with

formal methods with the purpose of generating the

comprehensive list of pathways of toxicity present in

a model. Indeed, through the use of biological prop-

erties, it is possible to deﬁne pathological states and

list all the paths leading to these states. The result-

ing paths shall ﬁnally be sorted thanks to additional

toxicological knowledge. Furthermore, ﬁltering the

resulting paths could also highlight gaps in the cur-

rent toxicological knowledge and help toxicologists

in their design of new experiments.

Finally, this formalism will serve as a basis to

develop a software platform dedicated to toxicology.

This platform is currently under development and it is

already possible to run simulations on biological ac-

tion networks. In the future, the platform will also be

able to integrate the temporal formulas and to ﬁlter

out paths from the biological action networks that do

not satisfy the formulas. This will be achieved gen-

erating all the paths allowed by a system biological

action network and directly checking these paths for

their biological relevance thanks to expressed biolog-

ical properties. Finally, by deﬁning states regarded as

pathologic, the platform will then be able to compute

all the paths leading to pathologic states and thus pro-

pose putative pathways of toxicity to toxicologists.

REFERENCES

Backhaus, T., Blanck, H., and Faust, M. (2010). Hazard and

risk assessment of chemical mixtures under REACH:

state of the art, gaps and options for improvement.

Swedish chemicals agency.

Bernot, G., Comet, J.-P., Richard, A., and Guespin, J.

(2004). Application of formal methods to biolog-

ical regulatory networks: extending thomas’ asyn-

chronous logical approach with temporal logic. Jour-

nal of theoretical biology, 229(3):339–347.

Brooks, G. A. (1998). Mammalian fuel utilization during

sustained exercise. Comparative Biochemistry and

Physiology Part B: Biochemistry and Molecular Bi-

ology, 120(1):89–107.

Calzone, L., Fages, F., and Soliman, S. (2006). Biocham:

an environment for modeling biological systems and

formalizing experimental knowledge. Bioinformatics,

22(14):1805–1807.

Chaki, S., Clarke, E. M., Ouaknine, J., Sharygina, N., and

Sinha, N. (2004). State/event-based software model

checking. In International Conference on Integrated

Formal Methods, pages 128–147. Springer.

Chance, B. and Williams, G. (1956). The respiratory chain

and oxidative phosphorylation. Adv Enzymol Relat Ar-

eas Mol Biol, 17:65–134.

Ciocchetta, F. and Hillston, J. (2009). Bio-pepa: A frame-

work for the modelling and analysis of biological sys-

BIOINFORMATICS 2017 - 8th International Conference on Bioinformatics Models, Methods and Algorithms

102

tems. Theoretical Computer Science, 410(33):3065–

3084.

Colborn, T., Dumanoski, D., and Myers, J. P. (1996). Our

stolen future: are we threatening our fertility, intelli-

gence, and survival?: a scientiﬁc detective story. Dut-

ton.

Field, J. B., Williams, H. E., and Mortimore, G. E.

(1963). Studies on the mechanism of ethanol-induced

hypoglycemia. Journal of Clinical Investigation,

42(4):497.

Hellerstein, M., Christiansen, M., Kaempfer, S., Kletke, C.,

Wu, K., Reid, J., Mulligan, K., Hellerstein, N., and

Shackleton, C. (1991). Measurement of de novo hep-

atic lipogenesis in humans using stable isotopes. Jour-

nal of Clinical Investigation, 87(5):1841.

Ito, S., Hagihara, S., and Yonezaki, N. (2014). A qualita-

tive framework for analysing homeostasis in gene net-

works. In proceedings of 5th International Conference

on Bioinformatics Models, Methods and Algorithms

(BIOINFORMATICS2014), pages 5–16.

Kepner, J. (2004). Synergy: The big unknowns of pesticide

exposure. Pesticides and You, 23(4):17–20.

Owen, O. E., Kalhan, S. C., and Hanson, R. W. (2002).

The key role of anaplerosis and cataplerosis for citric

acid cycle function. Journal of Biological Chemistry,

277(34):30409–30412.

Pettit, F. H., Pelley, J. W., and Reed, L. J. (1975). Reg-

ulation of pyruvate dehydrogenase kinase and phos-

phatase by acetyl-coa/coa and nadh/nad ratios. Bio-

chemical and biophysical research communications,

65(2):575–582.

Pilkis, S. J. and Granner, D. (1992). Molecular physiology

of the regulation of hepatic gluconeogenesis and gly-

colysis. Annual review of physiology, 54(1):885–909.

Schulz, H. (1991). Beta oxidation of fatty acids. Biochim-

ica et Biophysica Acta (BBA)-Lipids and Lipid

Metabolism, 1081(2):109–120.

Snoussi, E. H. (1989). Qualitative dynamics of piecewise-

linear differential equations: a discrete mapping ap-

proach. Dynamics and stability of Systems, 4(3-

4):565–583.

Talcott, C. (2008). Pathway logic. In International

School on Formal Methods for the Design of Com-

puter, Communication and Software Systems, pages

21–53. Springer.

Thomas, R. (1991). Regulatory networks seen as asyn-

chronous automata: a logical description. Journal of

theoretical biology, 153(1):1–23.

A Qualitative Framework Dedicated to Toxicology

103