Faulon J.L., Bender A. Handbook of Chemoinformatics Algorithms

Подождите немного. Документ загружается.

328 Handbook of Chemoinformatics Algorithms

11.5 REACTION NETWORK GENERATION ALGORITHM

The network generator algorithm presented next is based on the Dugundji–Ugi theory

[15] and the computer programs RAIN developed by Fontain and Reitsam [24] and

NetGen developed by Broadbelt et al. [4]. The generator is limited to monomolecular

and bimolecular reactions but could easily be extended to more complex reactions.

Following the Dugundji–Ugi methodology, reactants are represented by be-matrices.

The reactions are represented in the form of elementary transformations.As any reac-

tion, an elementary transformation is composed of two configurations containing

atoms participating in a reaction before and after the reaction has taken place. Unlike

full reactions, the configurations with elementary transformations do not take into

account the full structure of the species involved, but focuses on the reaction center

and its immediate surrounding. In other words, elementary transformations focus on

the electronic transformations that atoms undergo when reacting. One notes that the

codifications (RC numbers, be- and r-matrices, and signatures) presented in the Sec-

tion 11.3 can all be restricted to reaction centers, and thus can all be used to code

elementary transformations. Examples of elementary transformations using be- and

r-matrices are given in Figure 11.7 for thermal cracking. Fontain and Reitsan [24]

compiled the complete set of elementary transformations that carbon atoms can take

in organic reactions. This set is composed of 324 transformations; however, the set

can be greatly reduced depending on the studied process. For instance, it is known

[25] that hydrocarbon thermal cracking reactions can be generated from the five

elementary transformations listed in Figure 11.7. As another example, Susnow et al.

[7] proposed 17 transformations to model methane combustion. As a final example,

Hatzimanikatis et al. [26] have computed elementary transformation for metabolic

reactions. r-matrices were generated for all reactions stored in the metabolic KEGG

database (http://www.genome.jp/kegg/) resulting in about 250 unique elementary

transformations.

The input of the deterministic reaction network generator (Algorithm 11.1) is a

list of reactant species, a list of constraints, and a list of elementary transformations.

Constraints are user defined; examples of constraints are maximum number of species

and reactions generated, maximum species size, maximum number of reactants per

reaction, maximum number of lone electrons, and so on. The output of the algorithm

is a network composed of all possible species and reactions that can be created from

the input. The algorithm described below is illustrated for ethane thermal cracking in

Figure 11.8.

The algorithm starts by computing all the possible species (L

) that can be derived

by applying the elementary transformations (L

) to the initial species (L

) while

respecting the constraints. To prohibit duplication of species, L

is composed of

species that are not already present in L

.When applying elementary transformations,

one has to make the distinction between monomolecular and bimolecular reactions.

For each monomolecular elementary transformation, L

is composed of all the pos-

sible products that can be generated by applying the elementary transformation in all

possible ways to the species of L

. For each bimolecular reaction, L

is composed

of the products derived by applying the elementary transformations to all possible

pairs of species in L

. The list of reactions L

is updated each time an elementary

Reaction Network Generation 329

Homolysis

Step 1

–CH

={CH

–CH

}

•

}

•

Step 2

–CH

•

–CH

={CH

•

–CH

}

={CH

–CH

C=CH

•

–CH

C=CH

•

Step 3

•

–CH

•

+CH

•

–CH

•

–CH

+ CH

•

–CH

•

–CH

FIGURE 11.8 The first three steps of ethane thermal cracking using the elementary transfor-

mations listed in Figure 11.7. These matrices comprise only atoms for which the configuration

is changed. The lists of species (L

, L

, and L

) are given on the right side of the

figure.

transformation applies to a species in L

. Each reaction is represented by n-tuples

composed of the elementary transformation (from L

), the reacting species (from

), the product species (from L

or L

), and the rate constant of the reaction.

Once the list L

has been computed, the algorithm proceeds and computes recur-

sively L

, L

, ..., L

. Note that in order to compute L

(i > 1) one has to consider

monomolecular reactions only from the set L

si−1

since all monomolecular reactions

from the sets L

si−2

, L

si−3

, ...have already been computed in the previous steps. For

the same reason, the products of bimolecular reactions in L

are generated for every

possible pair of species so far generated having at least one element in L

si−1

. To avoid

redundancies of species in the lists L

, L

, ..., L

at any step i > 1, a new species

is added to L

only if the species is not in L

si−1

si−2

+···+L

The process halts at any step i, where the corresponding lists of species L

and

reaction L

are empty. Since bimolecular reactions can potentially create prod-

ucts of infinite molecular weight, in order to keep a finite network size one has

to set a limit for the maximum species size; let n be this limit. Note that with the

exception of polymerization reactions it is reasonable to assume that all species

in a given network will be limited in size. In turn, note that if the species have

a limited size, then the network generation algorithm is guaranteed to converge.

330 Handbook of Chemoinformatics Algorithms

Indeed, the maximum number of species, N, that can be formed is bounded by

the total number of molecular structures having a number of atoms raging from

1 to the specified value of n. This latter number is equal to the sum of the num-

bers of isomers containing 1, 2, ..., n atoms and scales exponentially with n even

for simple compounds such as hydrocarbons. Finally, note that the number of reac-

tions is also finite and is bounded by N

, which is an upper bound on the number

of ways of selecting four species (i.e., two reactants and two products) among

N species.

The network generator Algorithm 11.1 uses specific data structures to describe

species and reactions. A molecular species s is represented by a molecular graph

G(s). This graph can in turn be represented by any codification system presented in

Section 11.3, be-matrix for instance.An elementary transformation (et) is represented

by two molecular graphs, B(et) the et of the surrounding atoms before the reaction has

taken place, and E(et) the et after the reaction has taken place. These molecular graphs

can be represented by any of the codification systems presented in Section 11.3. A

reaction, r,isann-tuple composed of an elementary transformation, a list of reactant

species L

(r), a list of products L

(r), and a rate constant k(r). L

(r) is composed of

either one or two species for monomolecular and bimolecular reactions, respectively.

Additionally, the algorithm maintains several lists of species and reactions: the list of

reactants (L

), the list of species (L

) and reactions (L

) created at the current step

i, and the list of species generated at the previous step (L

si−1

). As already mentioned,

species are added to lists when not already present; since species are represented by

molecular graphs, the addition of species requires to check for graph isomorphism

between the species to be inserted and the species already in the list. Molecular

graph isomorphism algorithms can be found in Chapter 2. In Algorithm 11.1, the

function “species-constraints” verifies the user constraints and the function “rate-

constant” computes the rate constant of a given reaction using techniques described

in Section 11.4.

ALGORITHM 11.1 PSEUDOCODE FOR DETERMINISTIC

NETWORK GENERATOR

deterministic-network-generator(L

, L

)

input: - L

ys0

: list of initial species

-L

: list of elementary transformations

output: - L

: list of all species in network

-L

: list of all reactions in network

= L

, L

= Ø

si−1

= L

# L

si−1

is the list of species at previous

step i-1

, L

) = generate-species-reactions (L

si−1

, L

)

remove from L

all species present in L

remove from L

all reactions present in L

= L

, L

= L

si−1

= L

Reaction Network Generation 331

until (L

= Ø and L

= Ø)

generate-species-reactions(L

, L

)

input: - L

, L

: list of species

-L

: list of elementary transformations

output: - L

: list of species created at step i

-L

: list of reactions created at step i

For all tranformations et in L

For all species s

in L

and s

in L

# s

= Ø for monomolecular reactions

= generate-product(G(s

) +G(s

), et)

, L

) = update-species-reactions (et, s

, s

, L

)

done

return (L

, L

)

generate-product(Gb,et)

input: - G

: molecular graph of the reactant(s)

- et: elementary transformation

output: - L

: list of products

= Ø

B(et), E(et) = Graphs of the reactants and products

of et

For all subgraph b of G

s.t. b is isomorphic

to B(et)do

e = G

−b + E(et) # b is replaced in G

by E(et)

if species-constraints(e) and e is not in L

then L

= L

done

return L

update-species-reactions(et,s

, s

, L

)

input: - L

: list of graphs returned by

generate-product

- et: elementary transformation

- s

, s

: reactant species

output: - L

: list of species created at step i

-L

: list of reactions created at step i

= Ø, L

= Ø

For all graphs G

in L

compute L

the list of connected components of G

remove from L

all element already present

in L

= L

k = rate-constant(et, s

, s

, L

)

= L

+(et, s

, s

, L

, k)

done

return (L

, L

)

332 Handbook of Chemoinformatics Algorithms

Deterministic algorithms similar to the one described above have been applied to

processes such as pyrolysis, metabolism, and signal transduction. Broadbelt et al. [4]

developed a network generator named NetGen and applied it to hydrocarbon pyrol-

ysis (ethane and cyclohexane). Another pyrolysis network generation process can be

found in Ref. [5] for the thermal cracking of butane. With both applications, the ele-

mentary transformations and their associated be- and r-matrices are given in Figure

11.7. All studies generated networks with reactions rates. In Ref. [4], the rates were

estimated using linear free energy, and quantitative structure–reactivity relationships

as in Faulon and Sault [5]. Finally, in both cases, the numbers of species and reactions

were found to scale exponentially with the number of atoms of the initial species.

Beyond pyrolysis, it has been proposed to use techniques similar to NetGen to gen-

erate reaction networks to predict toxicity of complex reaction mixtures [27] focusing

on degradation of chemicals by cytochrome P450. The deterministic generator pro-

posed by Faulon and Sault has also been used to generate and analyze complex reaction

networks in interstellar chemistry [28].

Broadbelt et al. adapted NetGen with a new software named BNICE to explore

the diversity of complex metabolic networks [26]. With BNICE, elementary trans-

formations were computed for the KEGG database resulting in about 250 unique

elementary transformations. The transformations were then appliedto the biosynthesis

pathways of aromatic amino acids, phenylalanine, tyrosine, and tryptophan. The three

amino acids are synthesized from chrorismate and the cofactors and cosubstrates—

glutamate, gutamine, serine, NAD

/NAD, and 5-phospho-α-D-ribose-1-diphsopahe

(PRPP). The native pathways from chorismate to phenylalanine and tyrosine com-

prise three reactions between less than 10 compounds, and the native pathway leading

to tryptophan is composed of five reactions between 19 compounds. The 250 elemen-

tary transformations were applied to the initial substrate, cosubstrate, and cofactors,

producing 246 compounds with the phenylalanine pathway and 289 and 58 for the

tyrosine and tryptophan pathways. These compounds fall into three categories: the

compounds that are part of the native pathways (i.e., compounds in KEGG), com-

pounds found in the CAS database (http://www.cas.org/) but not in KEGG, and novel

compounds. The main outcome of the study was the in silico discovery of novel alter-

native biosynthesis pathways to aromatic amino acid biosynthesis, which remain

to be experimentally verified through enzyme and pathway engineering. BNICE

was also applied to polyketide biosynthesis [29]. While about 10,000 polyketides

structures have been discovered experimentally, BNICE raised this number over

7 millions.

To generate reaction networks for product biocatalysis and biodegradation, a semi-

automated system (UM-PPS) with a database (UM-BBD) has been developed at the

University of Minnesota [30] (http://umbbd.msi.umn.edu/predict/). In this system,

reaction rules (i.e., elementary transformation) are applied to an initial compound

entered by the user. Because several reaction rules applied to an initial compound

may result in different products, the user has to select the next product and the next

rule to apply. The process is iterated until no more rules can be applied, the user

selecting next products and rules at each step. At the time of writing (March 2009) the

database contained 259 biotransformation rules. Rules generally transform one func-

tional group into another, for instance, a cyano group can be hydrated with one water

Reaction Network Generation 333

molecule to form an amide or hydrolyzed by two water molecules to form carboxylic

acid and ammonia. Functional groups are searched through SMARTS matching and

transformations are carried using an algorithm whose outcome is identical to the

generate-product routine of Algorithm 11.1.

Deterministic network generators have also been used to create models char-

acterizing the dynamics of signal transduction networks. In particular, BioNetGen

(http://bionetgen.org) uses an algorithm similar toAlgorithm 11.1 with reaction rules

instead of elementary transformations [31].As seen in Section 11.3, elementary trans-

formations depict graph modifications of reactant species undergoing a reaction;

similarly, reaction rules code for modifications between reacting biological species

(ligand, protein, and complexes). BioNetGen and reaction rules have been used for

many biological applications and further information is given in Chapter 15.

As already mentioned, one of the main drawbacks of deterministic network gener-

ation is computational complexity due to the fact that an exponentially large number

of species may be generated. Techniques that overcome this drawback are discussed

in the next section.

11.6 REACTION NETWORK SAMPLING ALGORITHM

As discussed in Algorithm 11.1, the only reduction technique applicable to gen-

eral systems is the detailed reduction method. The caveat of the method, however,

is that the entire network must be known prior to reduction. Although the size of

the network generated by the deterministic algorithm is finite due to the limited

size of the species, as already mentioned the number of species in the network

grows exponentially with n, the number of atoms of the largest species. To over-

come the combinatorial explosion issue, network generation and reduction cannot be

processed in sequence if the goal is to derive a computationally tractable technique.

In this section, three methods are outlined where network generation and reduction

are performed simultaneously, concentration-sampling-network-generator (CSNG),

MC-sampling-network-generator (MCNG), and single-molecule simulator (SMG).

These algorithms are stochastic in nature and are efficient. The idea of the CSNG

algorithm was first published by Susnow et al. [7]. The MCNG algorithm was first

published by Faulon and Sault [5]. The SMG algorithm simulates the dynamics of

the species of an initial molecular graph comprising all reactants, without actually

generating a network. The idea of this algorithm was first published by Morton-Firth

and Bray [32], for predicting temporal fluctuations in an intracellular signaling path-

way. All the three algorithms require on-the-fly assignment of the rate constants of

the reactions generated; these computations are detailed in Section 11.4.

11.6.1 CONCENTRATION-SAMPLING NETWORK-GENERATOR ALGORITHM

With the concentration-sampling algorithm, the reduction of the network is based on

the species concentrations. The main assumption of this method is that if a species

created at any given step i has low concentration values over the reaction time, that

species will generate products with concentrations at the most equal to twice that low

value (in the particular case the species dissociates into two identical new species).

334 Handbook of Chemoinformatics Algorithms

Therefore, removing low-concentration species from the list L

should have a neg-

ligible effect on the final product distribution. Hence, in the present case, the total

number of network species is not limited, but the number of species generated at each

step is bounded by a predefined value M

. Note that this assumption is valid only if

the species are consumed during the reactions and do not act as catalysts. Indeed, if

a species is a catalyst even with low concentration, this species could have a rela-

tively large impact on the final product distribution. Hence, catalytic species must be

identified prior to using this scheme.

Using the above concentration assumption, the algorithm works as follows.At each

step i, the deterministic routine generate-species-reaction is run to compute the list of

new species. At this point the network is composed of all species in L

+···+L

and associated reactions. Although the network may not be complete, since the reac-

tion rates are calculated on-the-fly, the time evolution of the species concentrations

can be computed by solving the system of differential equations associated with the

partial network. The algorithm presented next uses the SSA MC–Gillespie algorithm

[22] (MC–Gillespie, cf. Section 11.4 for further details) to solve the system. The

MC–Gillespie integration (routine MC–Gillespie-step) is called M

times. That rou-

tine updates species concentration according to Equations 11.13 and 11.14. Precisely,

if r is the selected reaction, then the numbers of particles of all species in the list of

reactants for r (L

(r)) are decremented by 1 and the numbers of particles of all species

in the list of products for r (L

(r)) are incremented by 1. During the MC–Gillespie

integration, the CSNG algorithm retains for each species present in the network the

maximum concentration (i.e., maximum number of particles) reached over the time

period the system is integrated. This operation requires to maintain two lists: L

[s]

the list of species concentrations calculated for a given time step t, and L

[smax]

, the

list of species maximum concentrations calculated for all M

integration time steps.

Species are then sorted by decreasing concentration in L

[smax]

, the first M

elements

of the sorted list are kept in L

, while the others are eliminated. The numbers M

and

are user input. The algorithm is given in Algorithm 11.2, and the routine generate-

species-reaction is given in Algorithm 11.1. The algorithm 11.2 takes as input the

list of reactants (L

), reactant concentration (L

[s0]

), and elementary transformations

); it returns the list of species (L

) and reactions (L

) of the generated network.

ALGORITHM 11.2 PSEUDOCODE FOR CSNG

concentration-sampling-network-generator(L

, L

[s0]

, L

)

input: - L

: list of initial species (reactants)

-L

[s0]

: list of initial species concentrations

-L

: list of elementary transformations

output: - L

: list of all species in network

-L

: list of all reactions in network

= L

, L

=∅

si−1

= L

# L

si−1

is the list of species at previous

step i-1

do forever

, L

) = generate-species-reactions (L

si−1

, L

)

Reaction Network Generation 335

remove from L

all species present in L

remove from L

all reactions present in L

if (L

=∅ and L

=∅) then end

= L

, L

= L

, L

) = reduce-mechanism-concentration(L

, L

[s0]

, L

)

si−1

= L

done

reduce-mechanism-concentration(L

, L

[s0]

, L

)

input: - L

: list of species

-L

[s0]

: list of initial species (reactants)

concentrations

-L

: list of species created at step i

output: - L

: reduced list of species

-L

: reduced list of species created at step i

[s]

= L

[s0]

[smax]

= list of species maximum concentration

over time

initialized to Ø

t = 0

For M

steps do

[s]

, t) = MC-Gillespie-step(L

, L

[s]

, L

, t)

[smax]

= MAX(L

[smax]

, L

[s]

)

done

While (|L

| > M

) do

findsinL

having the lowest value in L

[smax]

= L

−s, L

= L

−s

done

return (L

, L

)

MC-Gillespie-step(L

, L

[s]

, L

, t)

input: - L

: list of species

-L

[s]

: list of species concentration

-L

: list of reactions

- t: time

output: - L

[s]

: updated list of species concentration

- τ: time after event occurs

τ = time of next event using eq. (11.13)

r = selected reaction in L

occurring at time

t +τ using eq.(11.14)

t = t +τ

(r) and L

(r) are the lists of reactants

and products for r

For all s in L

(r) do [s]=[s]−1 done

For all s in L

(r) do [s]=[s]+1 done

return (t, L

[s]

)

336 Handbook of Chemoinformatics Algorithms

The above algorithm was used by Faulon and Sault for thermal cracking of butane

[5]. The results given by the CSNG algorithm are identical to those of the deterministic

DNG algorithm when limiting the number of particles created at each step to no

more than 8 (for a full deterministic network comprising about 100 species and 1000

reactions).

A different implementation of the CSNG algorithm was developed by Susnow

et al. [7]. In this implementation, species concentrations are computed solving ODEs

rather than SSAs, and the selection of the species to keep when reducing the network is

based not on concentration but on rate of formation. In other words, inAlgorithm 11.2,

new species are kept in L

when their rate of formation is greater than a user-defined

threshold R

min

. Susnow et al. applied their method for the pyrolysis of ethane and

butane. The network generation process was iterated until the conversion of the initial

reactants was below a user-defined threshold. By changing the R

min

value, Susnow

et al. were able to reduce the number of reacted species and differential equations by

an order of magnitude while maintaining suitable complete mechanism.

11.6.2 MC-SAMPLING-NETWORK-GENERATOR ALGORITHM

The idea of the MC-sampling algorithm is to perform at the same time both the

MC integration and the network generation. The advantage of this technique is that

there are no assumptions regarding catalyst species. As in the previous case, one

starts with an initial reactant concentration (L

[s0]

) given in the form of numbers

of particles. At each step, the set of new species is computed using the generate-

species-reactions routine of Algorithm 11.1, but in the present case, these species are

generated by applying the ets only for species having nonzero concentration (i.e., set

s∗

in Algorithm 11.3). The concentrations of the new species are set to zero and

updated using the MC–Gillespie-step integration routine given in Algorithm 11.2.

The process is iterated until the number of steps exceeds a predefined M

value.

ALGORITHM 11.3 PSEUDOCODE FOR MCNG

MC-sampling-network-generator(L

, L

[s0]

, L

)

input: - L

: list of initial species (reactants)

-L

[s0]

: list of initial species concentrations

-L

: list of elementary transformations

output: - L

: list of all species in network

-L

: list of all reactions in network

= L

, L

=∅, L

[s]

= L

[s0]

t = 0

For M

steps do

∗

= species in L

with non-zero concentration

∗

= species in L

with non-zero concentration

, L

) = generate −species −reactions

∗

, L

∗

, L

)

remove from L

all species present in L

remove from L

all reactions present in L

Reaction Network Generation 337

= L

, L

= L

For all species s in L

do [s]=0 done

[s]

, t) = MC-Gillespie-step(L

, L

[s]

, L

, t)

done

Faulon and Sault applied the above scheme to generate the reaction network of

butane pyrolysis [5]. When comparing the results with the DNG algorithm, good

agreement was found as long as the number of initial particles M

was around

or above 2000; furthermore, the same threshold was found for all alkanes up to

octane.

MCNG has been implemented and used to model signaling and protein inter-

action networks. Faeder et al. [31] implemented MCMG within BioNetGen for

signaling cascades such as those found with Toll-like receptors. Lok and Brent

[11] developed Moleculizer (http://sourceforge.net/projects/moleculizer/), a software

implementing a procedure very much like the MCNG algorithm [33]. Moleculizer

has been used to model the receptor–G-protein complex and the MAP kinase cascade

and scaffold complex.

11.6.3 SMS ALGORITHM

This approach, which was first proposed by Monthon-Firth and Bray [32], led to the

development of a computer code named StochSim (http://www.ebi.ac.uk/∼lenov/

stochsim.html). So far it appears that this algorithm has been developed and used

exclusively to model biological network dynamics. Below, a general algorithm is

proposed, which is not necessarily optimized, but provides some ideas on how the

SMS approach could be applied to chemicals. In the StochSim method, the dynamics

of the reaction network is simulated without actually compiling species and species

concentrations. Instead, the algorithm considers each molecule as a separate entity, all

molecules have thus the same concentration (e.g., one particle), and reactions occurs

one molecule at a time. Starting with an initial simulation graph (G) comprising all

reactants, each reactant being duplicated a number of times equal to its initial concen-

tration, the main task of the algorithm is to apply elementary transformations following

an SSA procedure. Precisely, at each time step, one compiles all the reactions (L

) that

can occur according to the list of elementary transformations (L

) and the simulation

graph G. In their original paper, Monthon-Firth and Bray [32] provide equations where

the time increment (t) is fixed; since there is no guarantee that reactions can occur in

a fixed time interval, the dynamics may encompass null events (no reactions occur

between t and t + τ). In the implementation given below, Equations 11.13 and 11.14

are used to select a time delay t and a reaction r to apply to G. Let us note than when

solving these equations, all concentrations are equal to 1, and that Equation 11.12 does

not apply as each molecule is considered to be a separate entity. In Algorithm 11.4

given below, the reaction selected by equation 11.14 is fired at any subgraph (b)ofG

matching (e.g., isomorphic to) the et of the reactants (B(et)) of an elementary transfor-

mation. The subgraphs are compiled for all elementary transformations by the routine

generate-reactions.