reduce.pl

:- module(
  reduce,
  [
    reduce_model/1,
    reduce_model/0,
    delete_molecule/1,
    merge_molecules/2
  ]
).

:- use_module(ctl).

% Only for separate compilation/linting
:- use_module(nusmv).
:- use_module(util).

:- doc('This section describes commands to reduce a reaction model by deleting species, deleting reactions, merging species and merging reactions. We refer to \\cite{GSF10bi} for the use of these operations for detecting model reduction relationships in BioModels, and to \\cite{GFS14dam} for the theory of these graphical operations which are equivalent to subgraph epimorphisms (SEPI). We also provide a model reduction command which automatically finds a model reduction preserving the CTL specification of the behavior.').

:- devcom('FF: These commands should be generalized, 
first, to deletions and merge of molecules and reactions (easy SEPI linear for deletions quadratic for mergings) done!

second, to influence systems and hybrid reaction-influence systems (great SEPI unification),  See CMSB 2017 paper ;-)

and third, to hybrid attractors, CTL, FOLTL, trace specifications of the behavior.').

:- devdoc('\\section{Commands}').


reduce_model(Query) :-
    biocham_command,
    type(Query, ctl),
    doc(
	          'Deletes reactions as long as the specification of the behavior given by a CTL formula passed as
    argument remains satisfied in the current initial state.'
	      ),
    (
	      ctl_truth(Query,true)
     ->
	       reduce_model(Query, Removed),
	       format('removed ~p\n',[Removed])
     ;
     print('Error: the specification is not satisfied in the current model.')
    ).

reduce_model :-
    biocham_command,
    doc('Same as above using the current CTL
    specification of the behavior.'),
    findall(Formula, item([kind: ctl_spec, item: Formula]), Specs),
    join_with_op(Specs, '/\\', Spec),
    reduce_model(Spec).

:- doc('\\begin{example}').
:- biocham_silent(clear_model).
:- biocham(a=>b).
:- biocham(b=>c).
:- biocham(c=>d).
:- biocham(present(b)).
:- biocham(make_absent_not_present).
:- biocham(generate_ctl).
:- biocham(reduce_model).
:- doc('\\end{example}').  

:- doc("\\begin{example}
Reduction of the MAPK model with respect to the output reachability property only:
all dephosphorylation reactions can be removed, resulting in a different bifurcation diagram with full memory effect.
\\trace{
biocham: load('library:examples/mapk/mapk.bc').
biocham: reduce_model(reachable('MAPK~{p1,p2}')).
removed ['RAF-RAFK'=>'RAF'+'RAFK','RAF~{p1}'+'RAFPH'=>'RAF~{p1}-RAFPH','RAF~{p1}-RAFPH'=>'RAF~{p1}'+'RAFPH','MEK-RAF~{p1}'=>'MEK'+'RAF~{p1}','MEK~{p1}-RAF~{p1}'=>'MEK~{p1}'+'RAF~{p1}','MEKPH'+'MEK~{p1}'=>'MEK~{p1}-MEKPH','MEK~{p1}-MEKPH'=>'MEKPH'+'MEK~{p1}','MEKPH'+'MEK~{p1,p2}'=>'MEK~{p1,p2}-MEKPH','MEK~{p1,p2}-MEKPH'=>'MEKPH'+'MEK~{p1,p2}','MAPK-MEK~{p1,p2}'=>'MAPK'+'MEK~{p1,p2}','MAPK~{p1}-MEK~{p1,p2}'=>'MAPK~{p1}'+'MEK~{p1,p2}','MAPKPH'+'MAPK~{p1}'=>'MAPK~{p1}-MAPKPH','MAPK~{p1}-MAPKPH'=>'MAPKPH'+'MAPK~{p1}','MAPKPH'+'MAPK~{p1,p2}'=>'MAPK~{p1,p2}-MAPKPH','MAPK~{p1,p2}-MAPKPH'=>'MAPKPH'+'MAPK~{p1,p2}','RAF~{p1}-RAFPH'=>'RAF'+'RAFPH','MEK~{p1}-MEKPH'=>'MEK'+'MEKPH','MEK~{p1,p2}-MEKPH'=>'MEK~{p1}'+'MEKPH','MAPK~{p1}-MAPKPH'=>'MAPK'+'MAPKPH','MAPK~{p1,p2}-MAPKPH'=>'MAPK~{p1}'+'MAPKPH']
biocham: dose_response('RAFK',1e-5,1e-3, time:200, show:'MAPK~{p1,p2}').
}
\\end{example}
").

:- devcom(
	          'FF: todo with options sepi_reduce_model Deletes and/or merges molecules and/or reactions (aka subgraph epimorphism, SEPI reduction) as long as the specification of the behavior given by a CTL formula passed as
    argument remains satisfied in the current initial state.
').

%	       delete_molecules(Query, MoleculesRemoved),
%	       format('molecules removed ~p\n',[MoleculesRemoved]),
%	       merge_molecules(Query, MoleculesMerged),
%	       format('removed ~p\n',[MoleculesMerged]),
%	       delete_reactions(Query, ReactionsRemoved),
%	       format('removed ~p\n',[ReactionsRemoved]),
%	       merge_reactions(Query, ReactionsMerged),
%	       format('removed ~p\n',[ReactionsMerged])


delete_molecule(M):-
  biocham_command,
  doc('Deletes a molecule and its aliases from the current model, by replacing them by a parameter set to the initial value of the concentration.'),
  get_initial_concentration(M, C),
  \+ (
    item([kind: initial_state, id: Id, key:M]),
    \+ (
      delete_item(Id)
    )
  ),
  set_parameter(M, C),
  \+ (
    item([kind: reaction, id: Id, item: Reaction]),
    reaction(Reaction, [
      reactants: Reactants,
      inhibitors: Inhibitors,
      products: Products
    ]),
    (member(_ * M, Reactants) ; member(_ * M, Inhibitors) ; member(_ * M, Products)),
    \+ (
      delete_item(Id),
      rewrite_molecule_in_reaction(M, Reaction, SimplifiedReaction),
      add_item([kind: reaction, item: SimplifiedReaction])
    )
  ).


merge_molecules(M,N):-
  biocham_command,
  type(M, object),
  type(N, object),
  doc('
    merges two molecules by creating an alias making the second molecule canonical, and with initial concentration set to the sum of the initial concentrations of the merged molecules.
  '),
  devcom('FF should be the semantics of alias which should be synonymous'),
  get_initial_concentration(M, CM),
  get_initial_concentration(N, CN),
  alias(M=N),
  canonical(N),
  C is CM+CN,
  present(N,C).


:- devdoc('\\section{Public API}').


reduce_model(Query,Removed):-
    findall(Item, item([kind: reaction, item: Item]), Reactions),
    reduce_model(Reactions,Query, Removed).


:- devdoc('\\section{Internal predicates}').

reduce_model([R|Reactions],Query, Result) :-
  delete_item([kind: reaction, item: R]),
  item([kind: option, item: option(nusmv_initial_states: Init)]),
  item([kind: option, item: option(boolean_semantics: Bool)]),
  check_ctl_impl(Query, Init, no, Bool, true),
  !,
  Result=[R|Removed],
  reduce_model(Reactions,Query,Removed).

reduce_model([R|Reactions],Query, Removed) :-
  add_item([kind: reaction, item: R]),
  reduce_model(Reactions,Query,Removed).

reduce_model([],_,[]).