model_correction.pl 8.16 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
/*
Implement the automatic correction of models loaded from a sbml file
Coder: M. Hemery
*/

:- module(
  correction,
  [
  % Commands
  correct_model/0
  % API
  ]
).

% Only for separate compilation/linting
:- use_module(arithmetic_rules).
:- use_module(doc).
:- use_module(formal_derivation).
:- use_module(reaction_rules).
:- use_module(util).

correct_model :-
  biocham_command,
24
  doc('Improve wellformedness and strictness of a CRN model by removing reactions of null kinetics, splitting reactions that are implicitely reversible in two and adding correct annotation for modifiers.'),
25
  single_model(ModelId),
26
  setof(
27 28 29 30
    Molecule,
    identifier_kind(ModelId, Molecule, object),
    List_Molecule
  ),
31 32 33 34 35 36 37 38 39 40 41
  apply_to_all(List_Molecule, deal_with_null_kinetic),
  apply_to_all(List_Molecule, split_reaction(List_Molecule)),
  apply_to_all(List_Molecule, reverse_reaction),
  apply_to_all(List_Molecule, correct_modifiers(List_Molecule)).


%! apply_to_all(+List_Molecule, +Action)
%
% Extract from the model all the reaction and apply the action to all of them

apply_to_all(List_Molecule, Action) :-
42
  reactions_with_species(List_Molecule, List_Reaction),
43 44
  maplist(Action, List_Reaction).

45

46
%! deal_with_null_kinetic(+Reaction)
47
%
48 49 50
% Detect and modify reactions with identically null kinetics
% If the reaction has rate set to 0, it is surely a parsing error and is delete
% otherwise a new parameter called "p0_i" is introduce and set to 0.
51

52
deal_with_null_kinetic(Expression for Reaction) :-
53
  (
54
    Expression = 0
55 56
  ->
    delete_reaction(Expression for Reaction)
57 58 59 60 61 62 63 64 65
  ;
    zero_like(Value),
    is_present(Value, Expression)
  ->
    parameterize_zeros(Expression, New_Expression),
    delete_reaction(Expression for Reaction),
    add_reaction(New_Expression for Reaction)
  ;
    true
66 67
  ).

68 69 70 71 72 73 74 75 76 77 78 79 80 81
parameterize_zeros(Expression, New_Expression) :-
  (
    zero_like(Value),
    is_present(Value, Expression)
  ->
    new_parameter(New),
    set_parameter(New, Value),
    substitute_once(Value, New, Expression, Tempo_expression),
    parameterize_zeros(Tempo_expression, New_Expression)
  ;
    New_Expression = Expression
  ).


82
%! split_reaction(+List_Mol, +Reaction)
83 84 85
%
% Split a bidirection reaction in two

86
split_reaction(List_Mol, Expression for Reactant => Product) :-
87
  (
88
    ode:substitute_functions(Expression, Expr),
Mathieu Hemery's avatar
Mathieu Hemery committed
89 90
    % distribute(Expr, DistExpr),
    DistExpr = Expr,
91
    (
92
      DistExpr = _A + _B
93
    ->
94 95
      delete_reaction(Expression for Reactant => Product),
      split_reaction_sr(List_Mol, DistExpr for Reactant => Product)
96
    ;
97 98 99 100 101 102 103 104
      DistExpr = _Up - Down,
      % test if it is effectively a bidir. reaction
      member(P, List_Mol),
      models:formal_product(P, Reactant => Product),
      present_in_kinetics(P, Down)
    ->
      delete_reaction(Expression for Reactant => Product),
      split_reaction_sr(List_Mol, DistExpr for Reactant => Product)
105 106
    ;
      true
107
    )
108 109
  ).

110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
split_reaction_sr(List_Mol, Expression + Monomial for Reactant => Product) :-
  !,
  add_reaction(Monomial for Reactant => Product),
  split_reaction_sr(List_Mol, Expression for Reactant => Product).

split_reaction_sr(List_Mol, Expression - Monomial for Reactant => Product) :-
  !,
  (
    Reactant = Reac/Inhib
  ->
    add_reaction(Monomial for Product/Inhib => Reac)
  ;
    add_reaction(Monomial for Product => Reactant)
  ),
  split_reaction_sr(List_Mol, Expression for Reactant => Product).

split_reaction_sr(_List_Mol, Expression for Reactant => Product) :-
  add_reaction(Expression for Reactant => Product).


130 131 132 133 134 135 136 137 138
%! reverse_reaction(+Reaction)
%
% Reverse a reaction if its kinetic is negative

reverse_reaction(Expression for Reactant => Product) :-
  (
    ode:substitute_functions(Expression, - RevExpr)
  ->
    delete_reaction(Expression for Reactant => Product),
Mathieu Hemery's avatar
DNR  
Mathieu Hemery committed
139 140 141 142 143 144 145
    (
      Reactant = Reac/Inhib
    ->
      add_reaction(RevExpr for Product/Inhib => Reac)
    ;
      add_reaction(RevExpr for Product => Reactant)
    )
146 147
  ;
    true
148 149 150 151 152 153 154 155 156 157
  ).


%! correct_modifiers(+List_Molecule, +Reaction)
%
% Ensure that the modifiers are correctly formulated

correct_modifiers(List_Molecule, Expression_raw for Reaction) :-
  ode:substitute_functions(Expression_raw, Expression),
  detect_modifiers(List_Molecule, Expression, Catalysts, Inhibitors),
158 159 160
  % purify reaction
  findall(M, (member(M, List_Molecule), models:formal_catalyst(M, Reaction)), LC),
  remove_all_catalysts(LC, Reaction, Purified_Reaction),
161
  (
162 163 164 165 166
    %  append(Catalysts, Inhibitors, [])
    %->
    %true
    %;
    add_catalysts(Catalysts, Purified_Reaction, ReactionTempo),
167 168 169 170 171 172
    add_inhibitors(Inhibitors, ReactionTempo, NewReaction),
    delete_reaction(Expression_raw for Reaction),
    add_reaction(Expression_raw for NewReaction)
  ).


173 174 175 176 177
%! detect_modifiers(+List_Molecule, +Expression, -Catalyst, -Inhibitor)
%
% Search among the list of molecules for catalysts and inhibitors
% note that as only Expression is look at, reactants are sorted as catalysts

178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
detect_modifiers([], _Expression, [], []).

detect_modifiers([Mol|TailM], Expression, NewC, NewI) :-
  derivate(Expression, Mol, Deriv),
  (
    is_null(Deriv)
  ->
    NewC = TailC,
    NewI = TailI
  ;
    always_negative(Deriv)
  ->
    NewC = TailC,
    NewI = [Mol|TailI]
  ;
    NewC = [Mol|TailC],
    NewI = TailI
  ),
  detect_modifiers(TailM, Expression, TailC, TailI).

198 199 200 201 202

%! add_catalysts(+List_catalyst, +Reaction, -Modified_Reaction)
%
% modify reaction to add catalysts (let untouch if already present)

203
add_catalysts([], Reaction, Reaction).
204

205 206 207 208
add_catalysts([Head|Tail], Reaction, NewReaction) :-
  models:formal_reactant(Head, Reaction),
  !,
  add_catalysts(Tail, Reaction, NewReaction). 
209

210 211 212 213 214 215 216
add_catalysts([Head|Tail], Reactants => Products, NewReactants => NProducts+Head) :-
  add_catalysts(Tail, Reactants => Products, NReactants => NProducts),
  (
    NReactants = Reac / Inhib
  ->
    NewReactants = Reac+Head / Inhib
  ;
217
    NewReactants = NReactants+Head
218 219
  ).

220 221 222 223 224

%! add_inhibitors(+List_inhibitors, +Reaction, -Modified_Reaction)
%
% same as add_catalysts for inhibitors

225
add_inhibitors([], Reaction, Reaction).
226

227 228 229
add_inhibitors([Head|Tail], Reaction, NewReaction) :-
  models:formal_inhibitor(Head, Reaction),
  !,
230 231
  add_inhibitors(Tail, Reaction, NewReaction).

Mathieu Hemery's avatar
Mathieu Hemery committed
232 233 234 235 236 237 238 239 240
add_inhibitors([Head|Tail], Reactants => Products, NewReactants => NProducts) :-
  add_inhibitors(Tail, Reactants => Products, ReactantsTempo => ProductsTempo),
  (
    remove_molecule(ReactantsTempo, Head, NReactants),
    remove_molecule(ProductsTempo, Head, NProducts)
  ;
    NReactants = ReactantsTempo,
    NProducts = ProductsTempo
  ),
241 242 243 244 245
  (
    NReactants = Reac / Inhib
  ->
    NewReactants = Reac / (Head,Inhib)
  ;
Mathieu Hemery's avatar
Mathieu Hemery committed
246
    NewReactants = NReactants / Head
247
  ).
Mathieu Hemery's avatar
Mathieu Hemery committed
248

249 250 251
%! remove_molecule(+Species, +Molecule, -Without)
%
% Remove one occurence of Molecule from the Species input
Mathieu Hemery's avatar
Mathieu Hemery committed
252 253

remove_molecule(Molecule, Molecule, '_') :- !.
254

Mathieu Hemery's avatar
Mathieu Hemery committed
255 256 257 258 259 260 261 262 263 264 265 266 267
remove_molecule(P*Molecule, Molecule, Result) :- 
  (
    P = 1
  ->
    Result = '_'
  ;
    P = 2
  ->
    Result = Molecule
  ;
    Pm is P-1,
    Result = Pm*Molecule
  ),!.
268

Mathieu Hemery's avatar
Mathieu Hemery committed
269
remove_molecule(Species+Molecule, Molecule, Species) :- !.
270

Mathieu Hemery's avatar
Mathieu Hemery committed
271 272 273 274 275 276 277 278 279 280 281 282 283
remove_molecule(Species+P*Molecule, Molecule, Result) :-
  (
    P = 1
  ->
    Result = Species
  ;
    P = 2
  ->
    Result = Species+Molecule
  ;
    Pm is P-1,
    Result = Species+Pm*Molecule
  ),!.
284

Mathieu Hemery's avatar
Mathieu Hemery committed
285 286
remove_molecule(Molecule+Other, Molecule, Result) :- !,
  remove_molecule(Other+Molecule, Molecule, Result).
287

Mathieu Hemery's avatar
Mathieu Hemery committed
288
remove_molecule(P*Molecule+Other, Molecule, Result) :- !,
289 290 291 292 293 294 295 296 297
  remove_molecule(P*Molecule, Molecule, ResultTempo),
  (
    ResultTempo = '_'
  ->
    Result = Other
  ;
    Result = ResultTempo+Other
  ).

Mathieu Hemery's avatar
Mathieu Hemery committed
298 299
remove_molecule(Species+Other, Molecule, Without+Other) :-
  remove_molecule(Species, Molecule, Without).
300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325


%! remove_all_catalysts(+List_Catalyst, +Reaction, -NewReaction)
%
%

remove_all_catalysts([], Reaction, Reaction).

remove_all_catalysts([Head|Tail], R/I=>P, RR/I=>PP) :-
  remove_all_catalysts(Tail, R/I=>P, RT/I=>PT),
  remove_catalyst(RT, PT, Head, RR, PP).

remove_all_catalysts([Head|Tail], R=>P, RR=>PP) :-
  remove_all_catalysts(Tail, R=>P, RT=>PT),
  remove_catalyst(RT, PT, Head, RR, PP).

remove_catalyst(Reac, Prod, Cat, NReac, NProd) :-
  (
    remove_molecule(Reac, Cat, TReac),
    remove_molecule(Prod, Cat, TProd)
  ->
    remove_catalyst(TReac, TProd, Cat, NReac, NProd)
  ;
    NReac = Reac,
    NProd = Prod
  ).