component_contribution.py 11.9 KB
Newer Older
1
"""A wrapper for the GibbeEnergyPredictor in component-contribution."""
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
# The MIT License (MIT)
#
# Copyright (c) 2013 Weizmann Institute of Science
# Copyright (c) 2018 Institute for Molecular Systems Biology,
# ETH Zurich
# Copyright (c) 2018 Novo Nordisk Foundation Center for Biosustainability,
# Technical University of Denmark
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.


29
from typing import Dict, List, Tuple
Elad Noor's avatar
Elad Noor committed
30

31
import numpy as np
32
from component_contribution import GibbsEnergyPredictor
eladnoor's avatar
eladnoor committed
33

eladnoor's avatar
eladnoor committed
34
from . import FARADAY, R, default_I, default_pH, default_pMg, default_T, ureg
35
from .phased_reaction import Reaction
36

eladnoor's avatar
eladnoor committed
37

eladnoor's avatar
eladnoor committed
38
class ComponentContribution(object):
39 40 41
    """A wrapper class for GibbsEnergyPredictor.

    Also holds default conditions for compounds in the different phases.
42
    """
eladnoor's avatar
eladnoor committed
43

44
    predictor = GibbsEnergyPredictor()
eladnoor's avatar
eladnoor committed
45

eladnoor's avatar
eladnoor committed
46
    @ureg.check(None, None, None, '[concentration]', '[temperature]')
47
    def __init__(
eladnoor's avatar
eladnoor committed
48
        self,
49 50 51 52
        p_h: float = default_pH,
        p_mg: float = default_pMg,
        ionic_strength: float = default_I,
        temperature: float = default_T
53
    ) -> object:
54
        """Create a ComponentContribution object.
55

56 57 58 59 60
        :param p_h:
        :param p_mg:
        :param ionic_strength:
        :param temperature:
        """
eladnoor's avatar
eladnoor committed
61
        self.p_h = p_h
eladnoor's avatar
eladnoor committed
62
        self.ionic_strength = ionic_strength
eladnoor's avatar
eladnoor committed
63 64
        self.p_mg = p_mg
        self.temperature = temperature
eladnoor's avatar
eladnoor committed
65

66 67
    @property
    def RT(self) -> float:
68
        """Get the value of RT."""
69
        return R * self.temperature
70

71
    @staticmethod
72
    def standard_dg(reaction: Reaction) -> Tuple[float, float]:
73
        """Calculate the chemical reaction energies of a reaction.
74

75 76 77 78
        :param reaction: the input Reaction object
        :return: a tuple of the dG0 in kJ/mol and standard error. to
        calculate the confidence interval, use the range -1.96 to 1.96 times
        this value
eladnoor's avatar
eladnoor committed
79
        """
80
        residual_reaction, stored_dg = reaction.separate_stored_dg()
81

82 83 84 85 86 87 88 89
        standard_dg, dg_sigma = \
            ComponentContribution.predictor.standard_dg(
                residual_reaction)

        return standard_dg + stored_dg, dg_sigma

    @staticmethod
    def standard_dg_multi(
90
            reactions: List[Reaction]) -> Tuple[np.array, np.array]:
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
        """Calculate the chemical reaction energies of a list of reactions.

        Using the major microspecies of each of the reactants.

        :return: a tuple with the CC estimation of the major microspecies'
        standard formation energy, and the uncertainty matrix.
        """
        rxn_dg_pairs = map(lambda r: r.separate_stored_dg(), reactions)
        residual_reactions, stored_dg = zip(*rxn_dg_pairs)
        stored_dg = np.array(stored_dg)

        standard_dg, dg_sigma = \
            ComponentContribution.predictor.standard_dg_multi(
                residual_reactions)
        return standard_dg + stored_dg, dg_sigma
eladnoor's avatar
eladnoor committed
106

eladnoor's avatar
eladnoor committed
107 108
    def standard_dg_prime(
        self,
109
        reaction: Reaction
eladnoor's avatar
eladnoor committed
110
    ) -> Tuple[float, float]:
111
        """Calculate the transformed reaction energies of a reaction.
112

113 114 115 116
        :param reaction: the input Reaction object
        :return: a tuple of the dG0_prime in kJ/mol and standard error. to
        calculate the confidence interval, use the range -1.96 to 1.96 times
        this value
117
        """
eladnoor's avatar
eladnoor committed
118
        residual_reaction, stored_dg_prime = \
119
            reaction.separate_stored_dg_prime(
eladnoor's avatar
eladnoor committed
120 121 122 123
                p_h=self.p_h,
                ionic_strength=self.ionic_strength,
                temperature=self.temperature
            )
124

125 126 127 128 129 130 131
        standard_dg_prime, dg_sigma = \
            ComponentContribution.predictor.standard_dg_prime(
                residual_reaction,
                p_h=self.p_h,
                ionic_strength=self.ionic_strength,
                temperature=self.temperature
            )
eladnoor's avatar
eladnoor committed
132

133
        return standard_dg_prime + stored_dg_prime, dg_sigma
eladnoor's avatar
eladnoor committed
134

eladnoor's avatar
eladnoor committed
135 136
    def dg_prime(
        self,
137
        reaction: Reaction
138
    ) -> Tuple[float, float]:
139 140 141 142 143 144 145
        """Calculate the dG'0 of a single reaction.

        :param reaction: an object of type Reaction
        :return: a tuple of (dG_r_prime, dG_uncertainty) where dG_r_prime is
        the estimated Gibbs free energy of reaction and dG_uncertainty is the
        standard deviation of estimation. Multiply it by 1.96 to get a 95%
        confidence interval (which is the value shown on eQuilibrator).
146
        """
eladnoor's avatar
eladnoor committed
147
        dg_prime, dg_uncertainty = self.standard_dg_prime(reaction)
148
        dg_prime += self.RT * reaction.dg_correction()
eladnoor's avatar
eladnoor committed
149
        return dg_prime, dg_uncertainty
eladnoor's avatar
eladnoor committed
150

151 152
    def standard_dg_prime_multi(
        self,
153
        reactions: List[Reaction]
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
    ) -> Tuple[np.array, np.array]:
        """Calculate the transformed reaction energies of a list of reactions.

        :return: a tuple of two arrays. the first is a 1D NumPy array
        containing the CC estimates for the reactions' transformed dG0
        the second is a 2D numpy array containing the covariance matrix
        of the standard errors of the estimates. one can use the eigenvectors
        of the matrix to define a confidence high-dimensional space, or use
        U as the covariance of a Gaussian used for sampling
        (where dG0_cc is the mean of that Gaussian).
        """
        rxn_dg_pairs = map(lambda r: r.separate_stored_dg_prime(
            p_h=self.p_h, ionic_strength=self.ionic_strength,
            temperature=self.temperature), reactions)
        residual_reactions, stored_dg_primes = zip(*rxn_dg_pairs)
        stored_dg_primes = np.array(stored_dg_primes)

        standard_dg_prime, dg_sigma = \
            ComponentContribution.predictor.standard_dg_prime_multi(
                residual_reactions,
                p_h=self.p_h,
                ionic_strength=self.ionic_strength,
                temperature=self.temperature
            )
        return standard_dg_prime + stored_dg_primes, dg_sigma

eladnoor's avatar
eladnoor committed
180 181
    def physiological_dg_prime(
        self,
182
        reaction: Reaction
eladnoor's avatar
eladnoor committed
183
    ) -> Tuple[float, float]:
184
        """Calculate the dG'm of a single reaction.
185

186 187
        Assume all aqueous reactants are at 1 mM, gas reactants at 1 mbar and
        the rest at their standard concentration.
188

189
        :param reaction: an object of type Reaction
190 191 192 193
        :return: a tuple (dG_r_prime, dG_uncertainty) where dG_r_prime is
        the estimated Gibbs free energy of reaction and dG_uncertainty is the
        standard deviation of estimation. Multiply it by 1.96 to get a 95%
        confidence interval (which is the value shown on eQuilibrator).
eladnoor's avatar
eladnoor committed
194
        """
eladnoor's avatar
eladnoor committed
195 196
        physiological_dg_prime, dg_uncertainty = self.standard_dg_prime(
            reaction)
eladnoor's avatar
eladnoor committed
197 198
        physiological_dg_prime += \
            self.RT * reaction.physiological_dg_correction()
eladnoor's avatar
eladnoor committed
199
        return physiological_dg_prime, dg_uncertainty
eladnoor's avatar
eladnoor committed
200

201
    def ln_reversibility_index(self, reaction: Reaction) -> float:
202
        """Calculate the reversibility index (ln Gamma) of a single reaction.
203

204
        :return: ln_RI - The reversibility index (in natural log scale).
205
        """
eladnoor's avatar
eladnoor committed
206
        physiological_dg_prime, _ = self.physiological_dg_prime(reaction)
207

208 209 210
        abs_sum_coeff = reaction._GetAbsSumCoeff()
        if abs_sum_coeff == 0:
            return np.inf
eladnoor's avatar
eladnoor committed
211
        ln_RI = (2.0 / abs_sum_coeff) * physiological_dg_prime / self.RT
212
        return ln_RI
eladnoor's avatar
eladnoor committed
213

214
    def standard_e_prime(self, reaction: Reaction) -> Tuple[float, float]:
215 216
        """Calculate the E'0 of a single half-reaction.

217
        :param reaction: a Reaction object
218 219 220 221 222
        :return: a tuple (E0_prime, E0_uncertainty) where E0_prime is
        the estimated standard electrostatic potential of reaction and
        E0_uncertainty is the standard deviation of estimation. Multiply it
        by 1.96 to get a 95% confidence interval (which is the value shown on
        eQuilibrator).
eladnoor's avatar
eladnoor committed
223 224 225 226 227 228 229 230
        """
        n_e = reaction.check_half_reaction_balancing()
        if n_e is None:
            raise ValueError('reaction is not chemically balanced')
        if n_e == 0:
            raise ValueError('this is not a half-reaction, '
                             'electrons are balanced')

eladnoor's avatar
eladnoor committed
231
        dG0_prime, dG0_uncertainty = self.standard_dg_prime(reaction)
eladnoor's avatar
eladnoor committed
232

eladnoor's avatar
eladnoor committed
233
        E0_prime = -dG0_prime / (n_e*FARADAY)  # in Volts
234 235 236 237 238 239
        E0_uncertainty = dG0_uncertainty / abs(n_e*FARADAY)  # in Volts

        return E0_prime, E0_uncertainty

    def physiological_e_prime(
        self,
240
        reaction: Reaction
241 242 243
    ) -> Tuple[float, float]:
        """Calculate the E'0 of a single half-reaction.

244
        :param reaction: a Reaction object
245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266
        :return: a tuple (E0_prime, E0_uncertainty) where E0_prime is
        the estimated standard electrostatic potential of reaction and
        E0_uncertainty is the standard deviation of estimation. Multiply it
        by 1.96 to get a 95% confidence interval (which is the value shown on
        eQuilibrator).
        """
        n_e = reaction.check_half_reaction_balancing()
        if n_e is None:
            raise ValueError('reaction is not chemically balanced')
        if n_e == 0:
            raise ValueError('this is not a half-reaction, '
                             'electrons are balanced')

        dG0_prime, dG0_uncertainty = self.physiological_dg_prime(reaction)

        E0_prime = -dG0_prime / (n_e*FARADAY)  # in Volts
        E0_uncertainty = dG0_uncertainty / abs(n_e*FARADAY)  # in Volts

        return E0_prime, E0_uncertainty

    def e_prime(
        self,
267
        reaction: Reaction
268 269 270
    ) -> Tuple[float, float]:
        """Calculate the E'0 of a single half-reaction.

271
        :param reaction: a Reaction object
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
        :return: a tuple (E0_prime, E0_uncertainty) where E0_prime is
        the estimated standard electrostatic potential of reaction and
        E0_uncertainty is the standard deviation of estimation. Multiply it
        by 1.96 to get a 95% confidence interval (which is the value shown on
        eQuilibrator).
        """
        n_e = reaction.check_half_reaction_balancing()
        if n_e is None:
            raise ValueError('reaction is not chemically balanced')
        if n_e == 0:
            raise ValueError('this is not a half-reaction, '
                             'electrons are balanced')

        dG0_prime, dG0_uncertainty = self.dg_prime(reaction)

        E0_prime = -dG0_prime / (n_e*FARADAY)  # in Volts
        E0_uncertainty = dG0_uncertainty / abs(n_e*FARADAY)  # in Volts
eladnoor's avatar
eladnoor committed
289

eladnoor's avatar
eladnoor committed
290
        return E0_prime, E0_uncertainty
eladnoor's avatar
eladnoor committed
291

292 293
    def dg_analysis(
        self,
294
        reaction: Reaction
295 296 297
    ) -> List[Dict[str, object]]:
        """Get the analysis of the component contribution estimation process.

298
        :param reaction: a Reaction object.
299 300
        :return: the analysis results as a list of dictionaries
        """
301
        return ComponentContribution.predictor.get_dg_analysis(reaction)
eladnoor's avatar
eladnoor committed
302

303
    def is_using_group_contribution(self, reaction: Reaction) -> bool:
304 305
        """Check whether group contribution is needed to get this reactions' dG.

306
        :param reaction: a Reaction object.
307 308
        :return: true iff group contribution is needed
        """
eladnoor's avatar
eladnoor committed
309 310 311
        return ComponentContribution.predictor.is_using_group_contribution(
            reaction
        )