sufficient condition for rate independence

models.pl

@@ -422,7 +422,7 @@ change_item(Options, Kind, Key, Item) :-
       find_item([noinheritance, kind: Kind, key: Key, id: Id | Options]),
       replace_item(Id, Kind, Key, Item)
     ),
-    error(unknown_item),
+    error(unknown_item(_)),
     add_item([kind: Kind, key: Key, item: Item | Options])
   ).
 
@@ -478,7 +478,7 @@ find_item(Options) :-
   ->
     true
   ;
-    throw(error(unknown_item))
+    throw(error(unknown_item(Options)))
   ).
 
 

rate_independence.pl

@@ -6,7 +6,9 @@
    test_rate_independence_inputs_sinks/0,
    test_rate_independence/2,
    rate_independence_reduction_inputs_sinks/0,
-   rate_independence_reduction/2
+   rate_independence_reduction/2,
+				%API
+   add_merge_reaction/3
   ]
 ).
 
@@ -66,8 +68,7 @@ rate_independence_reduction(Inputs, Outputs):-
   write('inputs: '), writeln(Inputs),
   write('outputs: '), writeln(Outputs),
   simplify(Inputs, Outputs),
-  !,
-  writeln('There is no more reduction preserving rate independence.').
+  !.
 
 rate_independence_reduction_inputs_sinks:-
   biocham_command,
@@ -165,17 +166,16 @@ dfs_no_loop(GraphId, Vertex) :-
   assertz(seen(Vertex)).
 
 
-
 % Elimination of trivial reactions and isolated species with no outgoing fork
-simplify(Inputs, Outputs):-
-  enumerate_molecules(L),
-  subtract(L, Inputs, L1),
-  subtract(L1, Outputs, Molecules),
-  reduce(Molecules).
+simplify(_Inputs, Outputs):-
+  enumerate_molecules(Molecules),
+  subtract(Molecules, Outputs, Removable),
+  reduce(Removable).
 
 reduce(Molecules):-
   list_model,
-  \+ trivial_reaction,
+  \+ trivial_reaction(Molecules),
+  \+ trivial_loop(Molecules),
   \+ same_reaction,
   \+ (member(M, Molecules), eliminate(M)).
 
@@ -183,10 +183,25 @@ reduce(Molecules):-
   reduce(Molecules).
   
 
-% Elimination of trivial reactions with products=reactants
-trivial_reaction:-
+% 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),
+  forall(member((_*X), Reactants), member(X, Molecules)),
 write('--> removing trivial reaction '), writeln(Reaction),
   delete_item(Id).
 
@@ -200,162 +215,190 @@ same_reaction:-
   Id\=Id2,
   reaction(Reaction2, _, Reactants, _, Products),
 write('--> removing double reaction '), writeln(Reaction2),
-  delete_item(Id2).
+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):-
+  % fail, 
+  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).
 
-% Elimination of species M isolated in reactions alpha_i => beta_i + M  and M + gamma => delta 
-% if forall i  alpha_i\cap\gamma=\emptyset and \alpha_i \= \delta
-% by merging reactions resulting in alpha_i + gamma => beta_i + delta
+
+
+% Isolated non increasing catalyst in one reaction
 eliminate(M):-
-  reactions_with_reactant([M], [Reaction1]), 
-  reactions_with_product([M], Reactions), 
-  delete(Reactions, Reaction1, Reactions2),
+  reactions_with_reactant([M], [Reaction1]),
+  reactions_with_product([M], [Reaction1]),
 
-  (
-   Reactions=[Reaction1] % M catalyst not produced elsewhere
-  ->
-   reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
-   member((C*M), Reactants1),
-   member((D*M), Products1), 
-   C=D, % asking now for same stoichiometry (just a catalyst, no synthesis, no degradation)
-write('--> eliminating catalyst '), write(M), write(' isolated in reaction: '), writeln(Reaction1),
-   delete_reaction(Reaction1),
-   delete(Reactants1, (C*M), Reac1),
-   delete(Products1, (D*M), Prod1),
-   (
-    C<D
-   ->
-    E is D-C,
-    add_reaction(Kinetics1, Reac1, Inhibitors1, [(E*M) | Prod1], false)
-   ;
-    C>D
-   ->
-    E is C-D,
-    add_reaction(Kinetics1, [(E*M) | Reac1], Inhibitors1, Prod1, false)
-   ;
-    add_reaction(Kinetics1, Reac1, Inhibitors1, Prod1, false)
-   )
+  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),
   
-  ; % M isolated
-
-   reaction(Reaction1, _, _, _, Products1),
-   forall(
-	  member(Reaction2, Reactions2),
-	  (
-	   \+ (is_reactant(X, Reaction1), is_reactant(X, Reaction2)),
-	   \+ reaction(Reaction2, _, Products1, _, _)
-	  )
-	 ),
-
-write('--> eliminating '), writeln(M), write(' product of '), writeln(Reactions), write(' reactant of '), writeln(Reaction1),
-   delete_reaction(Reaction1),
-   
-   forall(
-	  member(Reaction2, Reactions2),
-	  (
-	   delete_reaction(Reaction2),
-	   add_merge_reactions(M, Reaction1, Reaction2)
-	  )
+  delete(Reactants1, (_*M), Reac1),
+  delete(Products1, (_*M), Prod1),
+  delete_reaction(Reaction1),
+  add_reaction(Kinetics1, Reac1, Inhibitors1, Prod1, false).
+
+
+
+  
+% Warning: Now restricted to GAMMA EMPTY and SAME STOICHIOMETRY 
+% Elimination of species M isolated in reactions alpha_i =>2 beta_i + M  and M (+ gamma) =>1 delta 
+% if forall i  alpha_i\cap\delta=\emptyset and \alpha_i\cap\gamma \= \empty
+% by merging reactions resulting in alpha_i (+ gamma) => beta_i + delta
+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 <=>
 
-% Todo: accept reversible reactions with syntax <=>
 
-% Elimination of M in 2 loops alpha 1<=>3 M (4<)=>2 gamma rewritten alpha (<)=> gamma
-% GENERALIZED TO alpha 1<=>3 beta+M (4<)=>2 gamma in alpha<=>gamma ELIMINATING M ISOLATED AND KEEPING REACTIONS WITH BETA
+% 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_reactant([M], [Reaction1, Reaction2]), % M isolated
-  reactions_with_product([M], [Reaction3 | Reactions]),
+  reactions_with_product([M], [Reaction1]),
+  reactions_with_reactant([M], [ReactionA, ReactionB]), 
+  
+  reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
+  member((C*M), Products1),
+  
   (
-   Reactions=[Reaction4] 
+   reaction(ReactionA, _, _, _, Reactants1)
   ->
-   reaction(Reaction4, Kinetics4, Reactants4, Inhibitors4, Products4),
-   member((S*M), Products4),
-   delete(Products4, (S*M), Beta)
+   Reaction2= ReactionA,
+   Reaction3= ReactionB
   ;
-   Reactions=[],
-   Reaction4=' '
+   reaction(ReactionB, _, _, _, Reactants1),
+   Reaction2= ReactionB,
+   Reaction3= ReactionA
   ),
 
-  
-  reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
-  member((S*M), Reactants1),
-
+  Reaction1\=Reaction2,
+    
   reaction(Reaction2, Kinetics2, Reactants2, Inhibitors2, Products2),
-  member((S*M), Reactants2),
-  
+  member((C*M), Reactants2),
   reaction(Reaction3, Kinetics3, Reactants3, Inhibitors3, Products3),
-  member((S*M), Products3),
+  member((C*M), Reactants3),
 
-  delete(Reactants1, (S*M), Beta),
-  delete(Reactants2, (S*M), Beta),
-  delete(Products3, (S*M), Beta),
+  delete(Reactants2, (C*M), Beta),
+  delete(Reactants3, (C*M), Beta),
+  delete(Products1, (C*M), Beta),
   
+  write('--> eliminating '), write(M), write(' in loop '), write(Reaction1),
+  write(' <--> '),write(Reaction2), write(' --> '), writeln(Reaction3),
+
   (
    Beta=[]
 
   -> % merge reactions
-   (
-    Products1=Reactants3
-   ->
-    add_merge_reactions(M, Reaction3, Reaction2),
-    (
-     Reactions=[]
-    ->
-     true
-    ;
-     Products2=Reactants4,
-     add_merge_reactions(M, Reaction4, Reaction1),
-     delete_reaction(Reaction4)
-    )
-   ; % 1 and 2 exchanged by pattern matching
-    Products2=Reactants3,
-    add_merge_reactions(M, Reaction3, Reaction1),
-    (
-     Reactions=[]
-    ->
-     true
-    ;
-     Products1=Reactants4,
-     add_merge_reactions(M, Reaction4, Reaction2),
-     delete_reaction(Reaction4)
-    )
-   )
+   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),
+  
+  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
-   (Products1=Reactants3 ; Products2=Reactants3),
-   add_reaction(Kinetics1, Beta, Inhibitors1, Products1, false),
    add_reaction(Kinetics2, Beta, Inhibitors2, Products2, false),
-   add_reaction(Kinetics3, Reactants3, Inhibitors3, Beta, false),
-   (
-    Reactions=[]
-   ->
-    true
-   ;
-    (
-     Products1=Reactants3
-    ->
-     Products2=Reactants4
-    ;
-     Products1=Reactants4
-    ),
-    add_reaction(Kinetics4, Reactants4, Inhibitors4, Beta, false),
-    delete_reaction(Reaction4)
-   )
+   add_reaction(Kinetics3, Beta, Inhibitors3, Products3, false),
+   add_reaction(Kinetics1, Reactants1, Inhibitors1, Beta, false),
+   add_reaction(Kinetics4, Reactants4, Inhibitors4, Beta, false)
   ),
-  write('--> Eliminating '), writeln(M),
-  write(' product of loop '), write(Reaction1), write(' and '),writeln(Reaction2),
-  write(' reactant of '), write(Reaction3), write(' '), writeln(Reaction4),
   
-  delete_reaction(Reaction1),
   delete_reaction(Reaction2),
-  delete_reaction(Reaction3).
+  delete_reaction(Reaction3),
+  delete_reaction(Reaction1),
+  delete_reaction(Reaction4).
+
    
 
+:- devdoc(' eliminate M in two reactions by merging them and not removing them ').
 
-add_merge_reactions(M, Reaction1, Reaction2):-
+add_merge_reaction(M, Reaction1, Reaction2):-
   reaction(Reaction1, Kinetics1, Reactants1, Inhibitors1, Products1),
   delete(Reactants1, (_*M), Reac1),
   delete(Products1, (_*M), Prod1),
@@ -366,7 +409,7 @@ add_merge_reactions(M, Reaction1, Reaction2):-
   delete(Products2, (_*M), Prod2),
   delete(Inhibitors2, M, Inhib2),
 
-  kinetics(Reac1, Inhib1, Kinetics1, Kin1), % M remains in kinetics
+  kinetics(Reac1, Inhib1, Kinetics1, Kin1), % M will remain in kinetics
   kinetics(Reac2, Inhib2, Kinetics2, Kin2),
   append(Reac1, Reac2, Reac),
   append(Prod1, Prod2, Prod),

rate_independence.plt

@@ -27,7 +27,7 @@ test('loop', [true(Reason = loop(_)), setup(max_model), cleanup(clear_model)]) :
 
 test('Max', [setup(max_model), cleanup(clear_model)]) :-
   with_output_to(atom(Result), test_rate_independence),
-  atom_concat('Current model is rate independent', _, Result).
+  atom_concat('The model is rate independent', _, Result).
 
 
 :- end_tests(rate_independence).

rate_independence_test.bc

@@ -4,6 +4,133 @@ prolog('writeln("%%%%%%%% WRONG ANSWER %%%%%%%%%")').
 prolog("writeln(' ')").
 
 
+
+prolog("writeln(' ')").
+prolog('writeln("%%%%%%%% UNDECIDED %%%%%%%%%")').
+prolog("writeln(' ')").
+
+
+
+
+prolog("writeln(' ')").
+clear_model.
+c+d=>e.
+e=>c+f.
+
+test_rate_independence_inputs_sinks.
+prolog('writeln("The model is rate independent: composite loop leakage with and-fork")').
+
+
+prolog("writeln(' ')").
+clear_model.
+a+b=>c+d.
+c+d=>a+b.
+c+d=>e.
+e=>c+f.
+
+test_rate_independence_inputs_sinks.
+prolog('writeln("The model is rate independent: composite loop leakage with and-fork")').
+
+
+
+prolog("writeln(' ')").
+prolog("writeln(' ')").
+prolog('writeln("%%%%%%% CORRECT ANSWER %%%%%%%%")').
+prolog("writeln(' ')").
+
+
+
+
+prolog("writeln(' ')").
+clear_model.
+c+d=>e.
+e=>c.
+e=>f.
+
+test_rate_independence_inputs_sinks.
+prolog('writeln("The model is rate independent: composite loop leakage with and-fork")').
+
+
+
+prolog("writeln(' ')").
+clear_model.
+a+b=>c+d.
+c+d=>e.
+e=>c.
+e=>f.
+
+test_rate_independence_inputs_sinks.
+prolog('writeln("The model is rate independent: partial loop leakage with or-fork")').
+
+prolog("writeln(' ')").
+clear_model.
+a+b=>c+d.
+c+d=>a+b.
+c+d=>e.
+e=>c.
+e=>f.
+
+test_rate_independence_inputs_sinks.
+prolog('writeln("The model is rate independent: partial loop leakage with or-fork")').
+
+
+prolog("writeln(' ')").
+clear_model.
+MA(1) for A=>B.
+MA(2) for B=>C.
+MA(5) for B=>D.
+MA(6) for D=>B.
+
+test_rate_independence_inputs_sinks. % no
+prolog('writeln("The model is rate independent: trivial loop leakage")').
+
+
+prolog("writeln(' ')").
+clear_model.
+MA(k1) for a+b => b.
+MA(k2) for b+c=>d.
+
+test_rate_independence({b,c}, {d}). % no
+prolog('writeln("The model is rate independent: without input a")').
+
+prolog("writeln(' ')").
+clear_model.
+a+b=>c+d.
+c+d=>a+b.
+c+d=>e.
+e=>c+d.
+e=>f+g.
+
+test_rate_independence_inputs_sinks.
+prolog('writeln("The model is rate independent: composite double loop leakage with composite or-fork")').
+
+prolog("writeln(' ')").
+clear_model.
+a+b=>c+d.
+c+d=>a+b.
+c+d=>e.
+e=>c+d+f.
+
+test_rate_independence_inputs_sinks.
+prolog('writeln("Undecided since no: DIV double loop (not by single loop)")').
+
+
+
+prolog("writeln(' ')").
+clear_model.
+input=>a.
+MA(k1) for a => b.
+MA(k2) for b => a.
+MA(k3) for b => c.
+MA(k4) for c => b.
+
+test_rate_independence(input, c). % no
+prolog('writeln("Undecided since no: output in loop equilibrium")').
+
+
+
+
+
 prolog("writeln(' ')").
 clear_model.
 d=>a.
@@ -11,7 +138,7 @@ a=>b.
 b=>2*d.
 a=>c.
 test_rate_independence_inputs_sinks.
-prolog('writeln("ANSWER NO AAAARGH NON CONFLUENCE OF RESULT / REACTION ORDERING AS SHOWN BELOW:")'). 
+prolog('writeln("Undecided since no: DIV")'). 
 
 
 prolog("writeln(' ')").
@@ -21,13 +148,7 @@ a=>b.
 b=>2*d.
 d=>a.
 test_rate_independence_inputs_sinks.
-prolog('writeln("ANSWER NO ")'). 
-
-
-
-prolog("writeln(' ')").
-prolog('writeln("%%%%%%%% DIVERGING SYLVAIN %%%%%%%%%")').
-prolog("writeln(' ')").
+prolog('writeln("Undecided since no: DIV ")'). 
 
 
 prolog("writeln(' ')").
@@ -36,7 +157,7 @@ a=>b.
 b=> a+c.
 list_model.
 test_rate_independence_inputs_sinks.
-prolog('writeln("ANSWER YES DIV ")').
+prolog('writeln("Undecided since no: DIV ")').
 
 
 prolog("writeln(' ')").
@@ -46,7 +167,7 @@ b=> a+c.
 d=> a.
 
 test_rate_independence_inputs_sinks.
-prolog('writeln("ANSWER YES DIV ")').
+prolog('writeln("Undecided since no: DIV ")').
 
 
 
@@ -57,7 +178,7 @@ a=>b.
 b=>2*a.
 
 test_rate_independence_inputs_sinks.
-prolog('writeln("ANSWER YES DIV ")').
+prolog('writeln("Undecided since no: DIV ")').
 
 
 prolog("writeln(' ')").
@@ -68,15 +189,10 @@ b=>2*d.
 d=>a.
 
 test_rate_independence_inputs_sinks.
-prolog('writeln("ANSWER YES DIV ")').
+prolog('writeln("Undecided since no: DIV ")').
 
 
 
-prolog("writeln(' ')").
-prolog('writeln("%%%%%%%% UNDECIDED %%%%%%%%%")').
-prolog("writeln(' ')").
-
-
 prolog("writeln(' ')").
 clear_model.
 a => x+c.
@@ -85,24 +201,9 @@ x+y => z.
 c+z => r.
 
 test_rate_independence_inputs_sinks.
-prolog('writeln("ANSWER YES: max = sum-min")').
+prolog('writeln("The model is rate independent: MAX = sum-min")').
 
 
-prolog("writeln(' ')").
-clear_mo