Commit 9cd83751 authored by Mathieu Hemery's avatar Mathieu Hemery
Browse files

Merge branch 'develop' of gitlab.inria.fr:lifeware/biocham into develop

parents ea0a05f7 247a870a
......@@ -3,13 +3,7 @@
[
% Commands
test_rate_independence/0,
test_rate_independence_inputs_sinks/0,
test_rate_independence/2,
rate_independence_reduction_inputs_sinks/0,
rate_independence_reduction/2,
test_rate_independence_invariants/0,
%API
add_merge_reaction/3
test_rate_independence_invariants/0
]
).
......@@ -44,49 +38,7 @@ test_rate_independence :-
format("Time: ~p s~n", [TTime]).
% to generalize to influence models too
% Test rate independence for given inputs and outputs by performing rate-independence preserving reduction
test_rate_independence(Inputs, Outputs):-
biocham_command,
type(Inputs, {object}),
type(Outputs, {object}),
doc('Test graphical sufficient conditions for rate independence of the current model (assuming well-formed kinetics)
for the computation of molecular species in \\argument{Outputs} from species in \\argument{Inputs}.
Warning: this command destroys the current model by applying the \\command{rate_independence_reduction}.'),
rate_independence_reduction(Inputs, Outputs),
test_rate_independence.
test_rate_independence_inputs_sinks :-
biocham_command,
doc('Test graphical sufficient conditions for rate independence of the current model for the computation of graph output species from graph input species (assuming well-formed kinetics). The input species are the species that are not a reaction product or a strict catalyst (i.e. with same stoichiometry). The output species are the species that are not reactant.'),
input_species(Inputs),
sink_species(Outputs),
test_rate_independence(Inputs, Outputs).
% Performs a model reduction which preserves the rate independence property for given inputs and outputs
rate_independence_reduction(Inputs, Outputs):-
biocham_command,
type(Inputs, {object}),
type(Outputs, {object}),
doc('Reduces the current reaction model while preserving rate independence for computing \\argument{Outputs} species from \\argument{Inputs} species.'),
write('inputs: '), writeln(Inputs),
write('outputs: '), writeln(Outputs),
simplify(Inputs, Outputs),
!.
rate_independence_reduction_inputs_sinks:-
biocham_command,
doc('Reduces the current reaction model while preserving the rate independence property for computing the graph sink species from the input species. Warning: the reduction does not preserve the computation, only the rate independence property which may be true or false.'),
input_species(Inputs),
sink_species(Outputs),
write('sources: '), writeln(Inputs),
write('sinks: '), writeln(Outputs),
rate_independence_reduction(Inputs, Outputs).
% API
%%% API
%! is_synthesis_free is det.
%
......@@ -172,271 +124,6 @@ dfs_no_loop(GraphId, Vertex) :-
assertz(seen(Vertex)).
:- dynamic(removableInputs/0).
:- dynamic(restrictionOnTrivialReaction/0).
restrictionOnTrivialReaction.
:- dynamic(noIsolatedReactant/0).
noIsolatedReactant.
:- dynamic(noBetaEmpty/0).
noBetaEmpty. % incorrect otherwise on eliminating d in loop a+b<=>c+d=>e=>c+b
% Elimination of trivial reactions and isolated species with no outgoing fork
simplify(Inputs, Outputs):-
enumerate_molecules(Molecules),
(removableInputs -> Mol=Molecules ; subtract(Molecules, Inputs, Mol)),
subtract(Mol, Outputs, Removable),
reduce(Removable).
reduce(Molecules):-
list_model,
\+ trivial_reaction(Molecules),
\+ trivial_loop(Molecules),
\+ same_reaction,
\+ (member(M, Molecules), eliminate(M)).
reduce(Molecules):-
reduce(Molecules).
% Isolated reactant in one trivial reaction loop
trivial_loop(Molecules):-
member(M, Molecules),
reactions_with_product([M], [Reaction1]),
reactions_with_reactant([M], [Reaction2]),
reaction(Reaction1, _, Reactants, _, [C*M]),
reaction(Reaction2, _, [C*M], _, Reactants),
forall(member(_*X, Reactants), member(X, Molecules)), % not on output
write('--> eliminating '), write(M), write(' isolated in trivial loop: '), write(Reaction1), write(' <--> '), writeln(Reaction2),
delete_reaction(Reaction1),
delete_reaction(Reaction2).
% Elimination of trivial reactions with products=reactants apart from inputs and outputs
trivial_reaction(Molecules):-
item([kind:reaction, id: Id, item: Reaction]),
reaction(Reaction, _Kinetics, Reactants, _Inhibitors, Reactants),
(restrictionOnTrivialReaction -> forall(member((_*X), Reactants), member(X, Molecules)); true),
write('--> removing trivial reaction '), writeln(Reaction),
delete_item(Id).
% Elimination of doublon reactions with same reactants and products
% Todo: the case of equivalent stiochiometry?
same_reaction:-
item([kind:reaction, id: Id, item: Reaction]),
item([kind:reaction, id: Id2, item: Reaction2]),
Id\=Id2,
reaction(Reaction, _Kinetics, Reactants, _Inhibitors, Products),
reaction(Reaction2, _, Reactants, _, Products),
write('--> removing double reaction '), writeln(Reaction2),
delete_item(Id2).
% Isolated reactant in one reaction
% to do? eliminates d in c+d=>e=>c+output and does not conclude to rate independence
% or not to do? does not reduce c+d=>e=>c e=>output and does not conclude to rate independence
eliminate(M):-
\+ noIsolatedReactant,
reactions_with_product([M], []),
reactions_with_reactant([M], [Reaction1]),
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
delete(Reactants1, (_*M), Reac1),
Reac1\=[],
write('--> eliminating '), write(M), write(' isolated reactant in: '), writeln(Reaction1),
delete_reaction(Reaction1),
add_reaction(Kinetics1, Reac1, Inhibitors1, Products1, false).
% Isolated non increasing catalyst in one reaction
eliminate(M):-
reactions_with_reactant([M], [Reaction1]),
reactions_with_product([M], [Reaction1]),
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
member(C*M, Reactants1),
member(D*M, Products1),
C>=D,
write('--> eliminating '), write(M), write(' isolated catalyst in: '), writeln(Reaction1),
delete(Reactants1, (_*M), Reac1),
delete(Products1, (_*M), Prod1),
delete_reaction(Reaction1),
add_reaction(Kinetics1, Reac1, Inhibitors1, Prod1, false).
% Elimination of species M isolated in reactions alpha_i =>2 beta_i + M and M =>1 gamma
% if forall i alpha_i\cap\gamma=\emptyset
% by merging reactions resulting in alpha_i => beta_i + gamma
eliminate(M):-
reactions_with_reactant([M], [Reaction1]),
reaction(Reaction1, _, [C*M], _, _Products1), % gamma empty
reactions_with_product([M], Reactions),
\+ member(Reaction1, Reactions),
forall(
member(Reaction2, Reactions),
(
reaction(Reaction2, _, _Reactants2, _, Products2),
member(C*M, Products2)
%forall(member(D*N, Reactants2), \+ member(D*N, Products1)) % no creation of trivial loop (harmful on output)
)
),
write('--> eliminating '), write(M), write(' reactant only in '), write(Reaction1), write(' product of '), writeln(Reactions),
forall(
member(Reaction2, Reactions),
(
add_merge_reaction(M, Reaction1, Reaction2),
delete_reaction(Reaction2)
)
),
delete_reaction(Reaction1).
% Todo: accept reversible reactions with syntax <=>
% Elimination of M in alpha 2<=>1 M =>3 gamma rewritten to alpha=>gamma
% Elimination of M in alpha 2<=>1 beta+M =>3 gamma rewritten to alpha<=>beta =>gamma
eliminate(M):-
reactions_with_product([M], [Reaction1]),
reactions_with_reactant([M], [ReactionA, ReactionB]),
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
member((C*M), Products1),
(
reaction(ReactionA, _, _, _, Reactants1)
->
Reaction2= ReactionA,
Reaction3= ReactionB
;
reaction(ReactionB, _, _, _, Reactants1),
Reaction2= ReactionB,
Reaction3= ReactionA
),
Reaction1\=Reaction2,
reaction(Reaction2, Kinetics2, Reactants2, Inhibitors2, Products2),
member((C*M), Reactants2),
reaction(Reaction3, Kinetics3, Reactants3, Inhibitors3, Products3),
member((C*M), Reactants3),
delete(Reactants2, (C*M), Beta),
delete(Reactants3, (C*M), Beta),
delete(Products1, (C*M), Beta),
(noBetaEmpty -> Beta=[]; true),
write('--> eliminating '), write(M), write(' in loop '), write(Reaction1),
write(' <--> '),write(Reaction2), write(' --> '), writeln(Reaction3),
(
Beta=[]
-> % merge reactions
add_merge_reaction(M, Reaction1, Reaction3)
; % eliminate M from reactions
add_reaction(Kinetics2, Beta, Inhibitors2, Products2, false),
add_reaction(Kinetics3, Beta, Inhibitors3, Products3, false),
add_reaction(Kinetics1, Reactants1, Inhibitors1, Beta, false)
),
delete_reaction(Reaction2),
delete_reaction(Reaction3),
delete_reaction(Reaction1).
% Elimination of M in alpha 2<=>1 M 4<=>3 gamma rewritten to alpha<=>gamma
% Elimination of M in alpha 2<=>1 beta+M 4<=>3 gamma rewritten to alpha<=>beta <=>gamma
eliminate(M):-
reactions_with_product([M], [Reaction1, Reaction4]),
reactions_with_reactant([M], [ReactionA, ReactionB]),
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
member((C*M), Products1),
(
reaction(ReactionA, _, _, _, Products1)
->
Reaction2= ReactionA,
Reaction3= ReactionB
;
reaction(ReactionB, _, _, _, Products1),
Reaction2= ReactionB,
Reaction3= ReactionA
),
reaction(Reaction2, Kinetics2, Reactants2, Inhibitors2, Products2),
member((C*M), Reactants2),
reaction(Reaction3, Kinetics3, Reactants3, Inhibitors3, Products3),
member((C*M), Reactants3),
reaction(Reaction4, Kinetics4, Reactants4, Inhibitors4, Products4),
member((C*M), Products4),
delete(Reactants2, (C*M), Beta),
delete(Reactants3, (C*M), Beta),
delete(Products1, (C*M), Beta),
delete(Products4, (C*M), Beta),
(noBetaEmpty -> Beta=[]; true),
write('--> eliminating '), write(M), write(' in double loop '), write(Reaction1),
write(' <--> '),write(Reaction2), write(Reaction3), write(' <--> '), writeln(Reaction4),
(
Beta=[]
-> % merge reactions
add_merge_reaction(M, Reaction1, Reaction3),
add_merge_reaction(M, Reaction4, Reaction2)
; % eliminate M from reactions
add_reaction(Kinetics2, Beta, Inhibitors2, Products2, false),
add_reaction(Kinetics3, Beta, Inhibitors3, Products3, false),
add_reaction(Kinetics1, Reactants1, Inhibitors1, Beta, false),
add_reaction(Kinetics4, Reactants4, Inhibitors4, Beta, false)
),
delete_reaction(Reaction2),
delete_reaction(Reaction3),
delete_reaction(Reaction1),
delete_reaction(Reaction4).
:- devdoc(' eliminate M in two reactions by merging them and not removing them ').
add_merge_reaction(M, Reaction1, Reaction2):-
reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
delete(Reactants1, (_*M), Reac1),
delete(Products1, (_*M), Prod1),
delete(Inhibitors1, M, Inhib1),
reaction(Reaction2, Kinetics2, Reactants2, Inhibitors2, Products2),
delete(Reactants2, (_*M), Reac2),
delete(Products2, (_*M), Prod2),
delete(Inhibitors2, M, Inhib2),
kinetics(Reac1, Inhib1, Kinetics1, Kin1), % M will remain in kinetics
kinetics(Reac2, Inhib2, Kinetics2, Kin2),
append(Reac1, Reac2, Reac),
append(Prod1, Prod2, Prod),
append(Inhib1, Inhib2, Inhib),
add_reaction(Kin1*Kin2,Reac,Inhib,Prod,false).
% Dynamic predicate for storing structural information of the graph
:- dynamic(p_invariant/1).
%% :- dynamic(sur_invariant/1).
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment