-
Notifications
You must be signed in to change notification settings - Fork 3
/
split_reversible_reactions.py
210 lines (178 loc) · 10.4 KB
/
split_reversible_reactions.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
""" Transform models.
:Author: Jonathan Karr <karr@mssm.edu>
:Author: Arthur Goldberg <Arthur.Goldberg@mssm.edu>
:Date: 2018-06-19
:Copyright: 2018, Karr Lab
:License: MIT
"""
from .core import Transform
from obj_tables.math.expression import ObjTablesTokenCodes
from wc_lang import Model, Reaction, DfbaObjReaction, RateLawDirection, FluxBounds, DfbaObjectiveExpression
from wc_onto import onto
from wc_utils.util.ontology import are_terms_equivalent
import copy
import math
import re
class SplitReversibleReactionsTransform(Transform):
""" Split reversible reactions in submodels into separate forward and backward reactions """
class Meta(object):
id = 'SplitReversibleReactions'
label = 'Split reversible reactions into separate forward and backward reactions'
def __init__(self, excluded_frameworks=None):
"""
Args:
excluded_frameworks (:obj:`list` of :obj:`pronto.term.Term`, optional): submodels using
these modeling integration frameworks are not transformed
"""
self.excluded_frameworks = excluded_frameworks
def run(self, model):
""" Split reversible reactions in submodels into separate forward and backward reactions
Args:
model (:obj:`Model`): model definition
Returns:
:obj:`Model`: same model definition, but with reversible reactions split into separate
forward and backward reactions, their flux bounds adjusted accordingly, and the
dFBA objective expression adjusted accordingly
"""
for submodel in model.submodels:
# skip submodels which use an excluded framework
if self.excluded_frameworks is not None:
if any([are_terms_equivalent(submodel.framework, excluded_framework)
for excluded_framework in self.excluded_frameworks]):
continue
for rxn in list(submodel.reactions):
if rxn.reversible:
# remove reversible reaction
model.reactions.remove(rxn)
submodel.reactions.remove(rxn)
# create separate forward and reverse reactions
rxn_for = submodel.reactions.create(
model=model,
id='{}_forward'.format(rxn.id),
name='{} (forward)'.format(rxn.name),
reversible=False,
evidence=rxn.evidence,
conclusions=rxn.conclusions,
identifiers=rxn.identifiers,
comments=rxn.comments,
references=rxn.references,
)
rxn_bck = submodel.reactions.create(
model=model,
id='{}_backward'.format(rxn.id),
name='{} (backward)'.format(rxn.name),
reversible=False,
evidence=rxn.evidence,
conclusions=rxn.conclusions,
identifiers=rxn.identifiers,
comments=rxn.comments,
references=rxn.references,
)
rxn.evidence = []
rxn.conclusions = []
rxn.identifiers = []
rxn.references = []
# copy participants and negate for backward reaction
for part in rxn.participants:
rxn_for.participants.append(part)
part_back = part.species.species_coefficients.get_one(coefficient=-1 * part.coefficient)
if part_back:
rxn_bck.participants.append(part_back)
else:
rxn_bck.participants.create(species=part.species, coefficient=-1 * part.coefficient)
rxn.participants = []
# copy rate laws
law_for = rxn.rate_laws.get_one(direction=RateLawDirection.forward)
law_bck = rxn.rate_laws.get_one(direction=RateLawDirection.backward)
if law_for:
law_for.reaction = rxn_for
law_for.direction = RateLawDirection.forward
law_for.id = law_for.gen_id()
if law_bck:
law_bck.reaction = rxn_bck
law_bck.direction = RateLawDirection.forward
law_bck.id = law_bck.gen_id()
# copy flux bounds
if are_terms_equivalent(submodel.framework, onto['WC:dynamic_flux_balance_analysis']):
if rxn.flux_bounds:
if not math.isnan(rxn.flux_bounds.min) and not math.isnan(rxn.flux_bounds.max):
# assume flux_bounds.min <= flux_bounds.max
assert rxn.flux_bounds.min <= rxn.flux_bounds.max, \
f"min flux bound greater than max in {rxn.id}"
# Mapping of flux bounds to backward and forward reactions
# Principles:
# lower bounds must be set, and cannot be negative because that would imply reversible
# upper bounds may be NaN if not previously set
# the forward reaction uses existing bounds that are positive
# the backward reaction uses existing bounds that are negative,
# swapping min and max and negating signs
# NaNs
# backward rxn forward rxn
# ------------ -----------
# min max min max min max
# ----- ----- ---- ---- ---- ----
# NaN NaN 0 NaN 0 NaN
# Values
# backward rxn forward rxn
# ------------ -----------
# min/max min max min max
# --------------- ---- ---- ---- ----
# min <= max <= 0 -max -min 0 0
# min <= 0 <= max 0 -min 0 max
# 0 <= min <= max 0 0 min max
backward_min = 0.
backward_max = float('NaN')
forward_min = 0.
forward_max = float('NaN')
if not math.isnan(rxn.flux_bounds.min):
if rxn.flux_bounds.min < 0:
backward_max = -rxn.flux_bounds.min
else:
backward_max = 0.
forward_min = rxn.flux_bounds.min
if not math.isnan(rxn.flux_bounds.max):
if 0 < rxn.flux_bounds.max:
forward_max = rxn.flux_bounds.max
else:
forward_max = 0.
backward_min = -rxn.flux_bounds.max
rxn_bck.flux_bounds = FluxBounds(min=backward_min,
max=backward_max,
units=rxn.flux_bounds.units)
rxn_for.flux_bounds = FluxBounds(min=forward_min,
max=forward_max,
units=rxn.flux_bounds.units)
# transform dFBA objective expression
# each dFBA objective expression is transformed for each reaction it uses
if rxn.dfba_obj_expression:
dfba_obj_expr = rxn.dfba_obj_expression
parsed_expr = dfba_obj_expr._parsed_expression
# create a new dFBA objective expression
# 1. use parsed_expr._obj_tables_tokens to recreate the expression and
# the objects it uses, while substituting the split reactions for the reversible reaction
new_obj_expr_elements = []
all_reactions = {Reaction: {},
DfbaObjReaction: {}}
for ot_token in parsed_expr._obj_tables_tokens:
if (ot_token.code == ObjTablesTokenCodes.obj_id and
issubclass(ot_token.model_type, Reaction) and
ot_token.model_id == rxn.id):
new_obj_expr_elements.append(f'({rxn_for.id} - {rxn_bck.id})')
all_reactions[Reaction][rxn_for.id] = rxn_for
all_reactions[Reaction][rxn_bck.id] = rxn_bck
continue
if (ot_token.code == ObjTablesTokenCodes.obj_id and
issubclass(ot_token.model_type, (Reaction, DfbaObjReaction))):
new_obj_expr_elements.append(ot_token.token_string)
all_reactions[ot_token.model_type][ot_token.model_id] = ot_token.model
continue
new_obj_expr_elements.append(ot_token.token_string)
new_obj_expr = ' '.join(new_obj_expr_elements)
# 2. create a new DfbaObjectiveExpression
dfba_obj_expr, error = DfbaObjectiveExpression.deserialize(new_obj_expr, all_reactions)
assert error is None, str(error)
rxn.dfba_obj_expression = None
rxn_for.dfba_obj_expression = dfba_obj_expr
rxn_bck.dfba_obj_expression = dfba_obj_expr
submodel.dfba_obj.expression = dfba_obj_expr
return model