From c8ddebdd7aff6b39596c30ac18d61b9c57a10013 Mon Sep 17 00:00:00 2001 From: teddygroves Date: Wed, 27 Nov 2024 12:21:30 +0100 Subject: [PATCH] WIP better ui --- src/enzax/kinetic_model.py | 69 +++- src/enzax/rate_equation.py | 20 +- src/enzax/rate_equations/__init__.py | 15 +- src/enzax/rate_equations/generalised_mwc.py | 187 +++++----- src/enzax/rate_equations/michaelis_menten.py | 358 +++++++++++-------- 5 files changed, 389 insertions(+), 260 deletions(-) diff --git a/src/enzax/kinetic_model.py b/src/enzax/kinetic_model.py index bb8e7fc..9f60d28 100644 --- a/src/enzax/kinetic_model.py +++ b/src/enzax/kinetic_model.py @@ -15,15 +15,59 @@ class KineticModelStructure(eqx.Module): """Structural information about a kinetic model.""" S: Float[Array, " s r"] - balanced_species: Int[Array, " n_balanced"] - unbalanced_species: Int[Array, " n_unbalanced"] + species: list[str] + reactions: list[str] + balanced_species: list[str] + species_to_dgf_ix: Int[Array, " s"] + balanced_species_ix: Int[Array, " b"] + unbalanced_species_ix: Int[Array, " u"] + + def __init__( + self, + S, + species, + reactions, + balanced_species, + species_to_dgf_ix, + ): + self.S = S + self.species = species + self.reactions = reactions + self.balanced_species = balanced_species + self.species_to_dgf_ix = species_to_dgf_ix + self.balanced_species_ix = jnp.array( + [i for i, s in enumerate(species) if s in balanced_species], + dtype=jnp.int16, + ) + self.unbalanced_species_ix = jnp.array( + [i for i, s in enumerate(species) if s not in balanced_species], + dtype=jnp.int16, + ) + + +class RateEquationKineticModelStructure(KineticModelStructure): + rate_equations: list[RateEquation] + + def __init__( + self, + S, + species, + reactions, + balanced_species, + species_to_dgf_ix, + rate_equations, + ): + super().__init__( + S, species, reactions, balanced_species, species_to_dgf_ix + ) + self.rate_equations = rate_equations class KineticModel(eqx.Module, ABC): """Abstract base class for kinetic models.""" parameters: PyTree - structure: KineticModelStructure + structure: KineticModelStructure = eqx.field(static=True) @abstractmethod def flux( @@ -44,7 +88,7 @@ def dcdt( """ # Noqa: E501 sv = self.structure.S @ self.flux(conc) - return sv[self.structure.balanced_species] + return sv[self.structure.balanced_species_ix] class RateEquationModel(KineticModel): @@ -64,10 +108,17 @@ def flux( """ # Noqa: E501 conc = jnp.zeros(self.structure.S.shape[0]) - conc = conc.at[self.structure.balanced_species].set(conc_balanced) - conc = conc.at[self.structure.unbalanced_species].set( + conc = conc.at[self.structure.balanced_species_ix].set(conc_balanced) + conc = conc.at[self.structure.unbalanced_species_ix].set( jnp.exp(self.parameters.log_conc_unbalanced) ) - t = [f(conc, self.parameters) for f in self.rate_equations] - out = jnp.array(t) - return out + flux_list = [] + for i, rate_equation in enumerate(self.structure.rate_equations): + ipt = rate_equation.get_input( + self.parameters, + i, + self.structure.S, + self.structure.species_to_dgf_ix, + ) + flux_list.append(rate_equation(conc, ipt)) + return jnp.array(flux_list) diff --git a/src/enzax/rate_equation.py b/src/enzax/rate_equation.py index aa89443..e3dbe24 100644 --- a/src/enzax/rate_equation.py +++ b/src/enzax/rate_equation.py @@ -1,10 +1,11 @@ """Module containing rate equations for enzyme-catalysed reactions.""" from abc import ABC, abstractmethod -from equinox import Module +import numpy as np +from equinox import Module from jaxtyping import Array, Float, PyTree, Scalar - +from numpy.typing import NDArray ConcArray = Float[Array, " n"] @@ -12,8 +13,19 @@ class RateEquation(Module, ABC): """Abstract definition of a rate equation. - A rate equation is an equinox [Module](https://docs.kidger.site/equinox/api/module/module/) with a `__call__` method that takes in a 1 dimensional array of concentrations and an arbitrary PyTree of parameters, returning a scalar value representing a single flux. + A rate equation is an equinox [Module](https://docs.kidger.site/equinox/api/module/module/) with a `__call__` method that takes in a 1 dimensional array of concentrations and an arbitrary PyTree of other inputs, returning a scalar value representing a single flux. """ # Noqa: E501 @abstractmethod - def __call__(self, conc: ConcArray, parameters: PyTree) -> Scalar: ... + def get_input( + self, + parameters: PyTree, + rxn_ix: int, + S: NDArray[np.float64], + species_to_dgf_ix: NDArray[np.int16], + ) -> PyTree: ... + + @abstractmethod + def __call__( + self, conc: ConcArray, rate_equation_input: PyTree + ) -> Scalar: ... diff --git a/src/enzax/rate_equations/__init__.py b/src/enzax/rate_equations/__init__.py index a0da47d..95f03d1 100644 --- a/src/enzax/rate_equations/__init__.py +++ b/src/enzax/rate_equations/__init__.py @@ -1,19 +1,18 @@ from enzax.rate_equations.michaelis_menten import ( ReversibleMichaelisMenten, IrreversibleMichaelisMenten, - MichaelisMenten, -) -from enzax.rate_equations.generalised_mwc import ( - AllostericReversibleMichaelisMenten, - AllostericIrreversibleMichaelisMenten, ) + +# from enzax.rate_equations.generalised_mwc import ( +# AllostericReversibleMichaelisMenten, +# AllostericIrreversibleMichaelisMenten, +# ) from enzax.rate_equations.drain import Drain AVAILABLE_RATE_EQUATIONS = [ ReversibleMichaelisMenten, IrreversibleMichaelisMenten, - MichaelisMenten, - AllostericReversibleMichaelisMenten, - AllostericIrreversibleMichaelisMenten, + # AllostericReversibleMichaelisMenten, + # AllostericIrreversibleMichaelisMenten, Drain, ] diff --git a/src/enzax/rate_equations/generalised_mwc.py b/src/enzax/rate_equations/generalised_mwc.py index 90f867e..224f588 100644 --- a/src/enzax/rate_equations/generalised_mwc.py +++ b/src/enzax/rate_equations/generalised_mwc.py @@ -1,94 +1,93 @@ -from jax import numpy as jnp -from jaxtyping import Array, Float, Int, PyTree, Scalar - -from enzax.rate_equation import ConcArray -from enzax.rate_equations.michaelis_menten import ( - IrreversibleMichaelisMenten, - MichaelisMenten, - ReversibleMichaelisMenten, -) - - -def generalised_mwc_effect( - conc_inhibitor: Float[Array, " n_inhibition"], - dc_inhibitor: Float[Array, " n_inhibition"], - conc_activator: Float[Array, " n_activation"], - dc_activator: Float[Array, " n_activation"], - free_enzyme_ratio: Scalar, - tc: Scalar, - subunits: int, -) -> Scalar: - """Get the allosteric effect on a rate. - - The equation is generalised Monod Wyman Changeux model as presented in Popova and Sel'kov 1975: https://doi.org/10.1016/0014-5793(75)80034-2. - - """ # noqa: E501 - qnum = 1 + jnp.sum(conc_inhibitor / dc_inhibitor) - qdenom = 1 + jnp.sum(conc_activator / dc_activator) - out = 1.0 / (1 + tc * (free_enzyme_ratio * qnum / qdenom) ** subunits) - return out - - -class GeneralisedMWC(MichaelisMenten): - """Mixin class for allosteric rate laws, assuming generalised MWC kinetics. - - See Popova and Sel'kov 1975 for the rate law: https://doi.org/10.1016/0014-5793(75)80034-2. - - Note that it is assumed there is a free_enzyme_ratio method available - that is why this is a subclass of MichaelisMenten rather than RateEquation. - """ # noqa: E501 - - subunits: int - tc_ix: int - ix_dc_activation: Int[Array, " n_activation"] - ix_dc_inhibition: Int[Array, " n_inhibition"] - species_activation: Int[Array, " n_activation"] - species_inhibition: Int[Array, " n_inhibition"] - - def get_tc(self, parameters: PyTree) -> Scalar: - return jnp.exp(parameters.log_transfer_constant[self.tc_ix]) - - def get_dc_activation(self, parameters: PyTree) -> Scalar: - return jnp.exp( - parameters.log_dissociation_constant[self.ix_dc_activation], - ) - - def get_dc_inhibition(self, parameters: PyTree) -> Scalar: - return jnp.exp( - parameters.log_dissociation_constant[self.ix_dc_inhibition], - ) - - def allosteric_effect(self, conc: ConcArray, parameters: PyTree) -> Scalar: - return generalised_mwc_effect( - conc_inhibitor=conc[self.species_inhibition], - conc_activator=conc[self.species_activation], - free_enzyme_ratio=self.free_enzyme_ratio(conc, parameters), - tc=self.get_tc(parameters), - dc_inhibitor=self.get_dc_inhibition(parameters), - dc_activator=self.get_dc_activation(parameters), - subunits=self.subunits, - ) - - -class AllostericIrreversibleMichaelisMenten( - GeneralisedMWC, IrreversibleMichaelisMenten -): - """A reaction with irreversible Michaelis Menten kinetics and allostery.""" - - def __call__(self, conc: Float[Array, " n"], parameters: PyTree) -> Scalar: - """The flux of an irreversible allosteric Michaelis Menten reaction.""" - allosteric_effect = self.allosteric_effect(conc, parameters) - non_allosteric_rate = super().__call__(conc, parameters) - return non_allosteric_rate * allosteric_effect - - -class AllostericReversibleMichaelisMenten( - GeneralisedMWC, - ReversibleMichaelisMenten, -): - """A reaction with reversible Michaelis Menten kinetics and allostery.""" - - def __call__(self, conc: ConcArray, parameters: PyTree) -> Scalar: - """The flux of an allosteric reversible Michaelis Menten reaction.""" - allosteric_effect = self.allosteric_effect(conc, parameters) - non_allosteric_rate = super().__call__(conc, parameters) - return non_allosteric_rate * allosteric_effect +# from jax import numpy as jnp +# from jaxtyping import Array, Float, Int, PyTree, Scalar + +# from enzax.rate_equation import ConcArray +# from enzax.rate_equations.michaelis_menten import ( +# IrreversibleMichaelisMenten, +# ReversibleMichaelisMenten, +# ) + + +# def generalised_mwc_effect( +# conc_inhibitor: Float[Array, " n_inhibition"], +# dc_inhibitor: Float[Array, " n_inhibition"], +# conc_activator: Float[Array, " n_activation"], +# dc_activator: Float[Array, " n_activation"], +# free_enzyme_ratio: Scalar, +# tc: Scalar, +# subunits: int, +# ) -> Scalar: +# """Get the allosteric effect on a rate. + +# The equation is generalised Monod Wyman Changeux model as presented in Popova and Sel'kov 1975: https://doi.org/10.1016/0014-5793(75)80034-2. + +# """ # noqa: E501 +# qnum = 1 + jnp.sum(conc_inhibitor / dc_inhibitor) +# qdenom = 1 + jnp.sum(conc_activator / dc_activator) +# out = 1.0 / (1 + tc * (free_enzyme_ratio * qnum / qdenom) ** subunits) +# return out + + +# class GeneralisedMWC: +# """Mixin class for allosteric rate laws, assuming generalised MWC kinetics. + +# See Popova and Sel'kov 1975 for the rate law: https://doi.org/10.1016/0014-5793(75)80034-2. + +# Note that it is assumed there is a free_enzyme_ratio method available - that is why this is a subclass of MichaelisMenten rather than RateEquation. +# """ # noqa: E501 + +# subunits: int +# tc_ix: int +# ix_dc_activation: Int[Array, " n_activation"] +# ix_dc_inhibition: Int[Array, " n_inhibition"] +# species_activation: Int[Array, " n_activation"] +# species_inhibition: Int[Array, " n_inhibition"] + +# def get_tc(self, parameters: PyTree) -> Scalar: +# return jnp.exp(parameters.log_transfer_constant[self.tc_ix]) + +# def get_dc_activation(self, parameters: PyTree) -> Scalar: +# return jnp.exp( +# parameters.log_dissociation_constant[self.ix_dc_activation], +# ) + +# def get_dc_inhibition(self, parameters: PyTree) -> Scalar: +# return jnp.exp( +# parameters.log_dissociation_constant[self.ix_dc_inhibition], +# ) + +# def allosteric_effect(self, conc: ConcArray, parameters: PyTree) -> Scalar: +# return generalised_mwc_effect( +# conc_inhibitor=conc[self.species_inhibition], +# conc_activator=conc[self.species_activation], +# free_enzyme_ratio=self.free_enzyme_ratio(conc, parameters), +# tc=self.get_tc(parameters), +# dc_inhibitor=self.get_dc_inhibition(parameters), +# dc_activator=self.get_dc_activation(parameters), +# subunits=self.subunits, +# ) + + +# class AllostericIrreversibleMichaelisMenten( +# GeneralisedMWC, IrreversibleMichaelisMenten +# ): +# """A reaction with irreversible Michaelis Menten kinetics and allostery.""" + +# def __call__(self, conc: Float[Array, " n"], parameters: PyTree) -> Scalar: +# """The flux of an irreversible allosteric Michaelis Menten reaction.""" +# allosteric_effect = self.allosteric_effect(conc, parameters) +# non_allosteric_rate = super().__call__(conc, parameters) +# return non_allosteric_rate * allosteric_effect + + +# class AllostericReversibleMichaelisMenten( +# GeneralisedMWC, +# ReversibleMichaelisMenten, +# ): +# """A reaction with reversible Michaelis Menten kinetics and allostery.""" + +# def __call__(self, conc: ConcArray, parameters: PyTree) -> Scalar: +# """The flux of an allosteric reversible Michaelis Menten reaction.""" +# allosteric_effect = self.allosteric_effect(conc, parameters) +# non_allosteric_rate = super().__call__(conc, parameters) +# return non_allosteric_rate * allosteric_effect diff --git a/src/enzax/rate_equations/michaelis_menten.py b/src/enzax/rate_equations/michaelis_menten.py index 55cf1d0..91011f8 100644 --- a/src/enzax/rate_equations/michaelis_menten.py +++ b/src/enzax/rate_equations/michaelis_menten.py @@ -1,110 +1,108 @@ -from abc import abstractmethod +import equinox as eqx +import numpy as np from jax import numpy as jnp from jaxtyping import Array, Float, Int, PyTree, Scalar +from numpy.typing import NDArray + +from enzax.rate_equation import RateEquation + + +class IrreversibleMichaelisMentenInput(eqx.Module): + kcat: Scalar + enzyme: Scalar + ix_ki_species: NDArray[np.int16] + ki: Float[Array, " n_ki"] + ix_substrate: NDArray[np.int16] + substrate_kms: Float[Array, " n_substrate"] + substrate_stoichiometry: NDArray[np.float64] + + +def get_irreversible_michaelis_menten_input( + parameters: PyTree, + rxn_ix: int, + S: NDArray[np.float64], + species_to_dgf_ix: NDArray[np.int16], + ci_ix: NDArray[np.int16], +) -> IrreversibleMichaelisMentenInput: + Sj = S[:, rxn_ix] + ix_substrate = np.argwhere(Sj < 0.0) + return IrreversibleMichaelisMentenInput( + kcat=jnp.exp(parameters.log_kcat[rxn_ix]), + enzyme=jnp.exp(parameters.log_enzyme[rxn_ix]), + ix_substrate=ix_substrate, + substrate_kms=jnp.exp(parameters.log_km[rxn_ix]), + substrate_stoichiometry=Sj[ix_substrate], + ix_ki_species=ci_ix, + ki=parameters.jnp.exp(parameters.log_ki[rxn_ix]), + ) + -from enzax.rate_equation import RateEquation, ConcArray +class ReversibleMichaelisMentenInput(eqx.Module): + kcat: Scalar + enzyme: Scalar + ki: Float[Array, " n_ki"] + substrate_kms: Float[Array, " n_substrate"] + product_kms: Float[Array, " n_product"] + dgf: Float[Array, " n_reactant"] + temperature: Scalar + ix_ki_species: NDArray[np.int16] + ix_reactant: NDArray[np.int16] + ix_substrate: NDArray[np.int16] + ix_product: NDArray[np.int16] + reactant_stoichiometry: NDArray[np.float64] + substrate_stoichiometry: NDArray[np.float64] + product_stoichiometry: NDArray[np.float64] + water_stoichiometry: float + + +def get_reversible_michaelis_menten_input( + parameters: PyTree, + rxn_ix: int, + S: NDArray[np.float64], + species_to_dgf_ix: NDArray[np.int16], + ci_ix: NDArray[np.int16], + water_stoichiometry: float, +) -> ReversibleMichaelisMentenInput: + Sj = S[:, rxn_ix] + ix_reactant = np.argwhere(Sj != 0.0) + ix_substrate = np.argwhere(Sj < 0.0) + ix_product = np.argwhere(Sj > 0.0) + return ReversibleMichaelisMentenInput( + kcat=jnp.exp(parameters.log_kcat[rxn_ix]), + enzyme=jnp.exp(parameters.log_enzyme[rxn_ix]), + substrate_kms=jnp.exp(parameters.log_km[rxn_ix][0]), + product_kms=jnp.exp(parameters.log_km[rxn_ix][1]), + ki=jnp.exp(parameters.log_ki[rxn_ix]), + dgf=jnp.exp(parameters.dgf[species_to_dgf_ix][ix_reactant]), + temperature=parameters.temperature, + ix_ki_species=ci_ix, + ix_reactant=ix_reactant, + ix_substrate=ix_substrate, + ix_product=ix_product, + reactant_stoichiometry=Sj[ix_reactant], + substrate_stoichiometry=Sj[ix_substrate], + product_stoichiometry=Sj[ix_product], + water_stoichiometry=water_stoichiometry, + ) def numerator_mm( - conc: ConcArray, - km: Float[Array, " n"], - ix_substrate: Int[Array, " n_substrate"], - substrate_km_positions: Int[Array, " n_substrate"], + substrate_conc: Float[Array, " n_substrate"], + substrate_kms: Int[Array, " n_substrate"], ) -> Scalar: """Get the product of each substrate's concentration over its km. This quantity is the numerator in a Michaelis Menten reaction's rate equation """ # Noqa: E501 - return jnp.prod((conc[ix_substrate] / km[substrate_km_positions])) - - -class MichaelisMenten(RateEquation): - """Base class for Michaelis Menten rate equations. - - Subclasses need to implement the __call__ and free_enzyme_ratio methods. - - """ - - kcat_ix: int - enzyme_ix: int - km_ix: Int[Array, " n"] - ki_ix: Int[Array, " n_ki"] - reactant_stoichiometry: Float[Array, " n"] - ix_substrate: Int[Array, " n_substrate"] - ix_ki_species: Int[Array, " n_ki"] - substrate_km_positions: Int[Array, " n_substrate"] - substrate_reactant_positions: Int[Array, " n_substrate"] - - def get_kcat(self, parameters: PyTree) -> Scalar: - return jnp.exp(parameters.log_kcat[self.kcat_ix]) - - def get_km(self, parameters: PyTree) -> Scalar: - return jnp.exp(parameters.log_km[self.km_ix]) - - def get_ki(self, parameters: PyTree) -> Scalar: - return jnp.exp(parameters.log_ki[self.ki_ix]) - - def get_enzyme(self, parameters: PyTree) -> Scalar: - return jnp.exp(parameters.log_enzyme[self.enzyme_ix]) - - @abstractmethod - def free_enzyme_ratio( - self, - conc: ConcArray, - parameters: PyTree, - ) -> Scalar: ... - - -def free_enzyme_ratio_imm( - conc_sub: Float[Array, " n_substrate"], - km_sub: Float[Array, " n_substrate"], - stoich_sub: Float[Array, " n_substrate"], - ki: Float[Array, " n_ki"], - conc_inhibitor: Float[Array, " n_ki"], -) -> Scalar: - """Free enzyme ratio for irreversible Michaelis Menten reactions.""" - return 1.0 / ( - jnp.prod(((conc_sub / km_sub) + 1) ** jnp.abs(stoich_sub)) - + jnp.sum(conc_inhibitor / ki) - ) - - -class IrreversibleMichaelisMenten(MichaelisMenten): - """A reaction with irreversible Michaelis Menten kinetics.""" - - def free_enzyme_ratio(self, conc, parameters): - return free_enzyme_ratio_imm( - conc_sub=conc[self.ix_substrate], - km_sub=self.get_km(parameters)[self.substrate_km_positions], - ki=self.get_ki(parameters), - conc_inhibitor=conc[self.ix_ki_species], - stoich_sub=self.reactant_stoichiometry[ - self.substrate_reactant_positions - ], - ) - - def __call__(self, conc: Float[Array, " n"], parameters: PyTree) -> Scalar: - """Get flux of a reaction with irreversible Michaelis Menten kinetics.""" # noqa: E501 - kcat = self.get_kcat(parameters) - enzyme = self.get_enzyme(parameters) - km = self.get_km(parameters) - numerator = numerator_mm( - conc=conc, - km=km, - ix_substrate=self.ix_substrate, - substrate_km_positions=self.substrate_km_positions, - ) - free_enzyme_ratio = self.free_enzyme_ratio(conc, parameters) - return kcat * enzyme * numerator * free_enzyme_ratio + return jnp.prod((substrate_conc / substrate_kms)) def get_reversibility( - conc: Float[Array, " n"], - water_stoichiometry: Scalar, + reactant_conc: Float[Array, " n_reactant"], dgf: Float[Array, " n_reactant"], temperature: Scalar, - reactant_stoichiometry: Float[Array, " n_reactant"], - ix_reactants: Int[Array, " n_reactant"], + reactant_stoichiometry: NDArray[np.float64], + water_stoichiometry: float, ) -> Scalar: """Get the reversibility of a reaction. @@ -113,81 +111,151 @@ def get_reversibility( """ RT = temperature * 0.008314 dgf_water = -150.9 - dgr = reactant_stoichiometry @ dgf + water_stoichiometry * dgf_water - quotient = reactant_stoichiometry @ jnp.log(conc[ix_reactants]) - out = 1.0 - jnp.exp(((dgr + RT * quotient) / RT)) + dgr = ( + reactant_stoichiometry.T @ dgf + water_stoichiometry * dgf_water + ).flatten() + quotient = (reactant_stoichiometry.T @ jnp.log(reactant_conc)).flatten() + out = 1.0 - jnp.exp((dgr + RT * quotient) / RT)[0] return out +def free_enzyme_ratio_imm( + substrate_conc: Float[Array, " n_substrate"], + substrate_km: Float[Array, " n_substrate"], + ki: Float[Array, " n_ki"], + inhibitor_conc: Float[Array, " n_ki"], + substrate_stoichiometry: NDArray[np.float64], +) -> Scalar: + """Free enzyme ratio for irreversible Michaelis Menten reactions.""" + return 1.0 / ( + jnp.prod( + ((substrate_conc / substrate_km) + 1) + ** jnp.abs(substrate_stoichiometry) + ) + + jnp.sum(inhibitor_conc / ki) + ) + + def free_enzyme_ratio_rmm( - conc_sub: Float[Array, " n_substrate"], - km_sub: Float[Array, " n_substrate"], - stoich_sub: Float[Array, " n_substrate"], - conc_prod: Float[Array, " n_product"], - km_prod: Float[Array, " n_prod"], - stoich_prod: Float[Array, " n_prod"], - conc_inhibitor: Float[Array, " n_ki"], + substrate_conc: Float[Array, " n_substrate"], + product_conc: Float[Array, " n_product"], + substrate_kms: Float[Array, " n_substrate"], + product_kms: Float[Array, " n_product"], + inhibitor_conc: Float[Array, " n_ki"], ki: Float[Array, " n_ki"], + substrate_stoichiometry: NDArray[np.float64], + product_stoichiometry: NDArray[np.float64], ) -> Scalar: """The free enzyme ratio for a reversible Michaelis Menten reaction.""" return 1.0 / ( -1.0 - + jnp.prod(((conc_sub / km_sub) + 1.0) ** jnp.abs(stoich_sub)) - + jnp.prod(((conc_prod / km_prod) + 1.0) ** jnp.abs(stoich_prod)) - + jnp.sum(conc_inhibitor / ki) + + jnp.prod( + ((substrate_conc / substrate_kms) + 1.0) + ** jnp.abs(substrate_stoichiometry) + ) + + jnp.prod( + ((product_conc / product_kms) + 1.0) + ** jnp.abs(product_stoichiometry) + ) + + jnp.sum(inhibitor_conc / ki) ) -class ReversibleMichaelisMenten(MichaelisMenten): +class IrreversibleMichaelisMenten(RateEquation): + """A reaction with irreversible Michaelis Menten kinetics.""" + + competitive_inhibitor_ix: NDArray[np.int16] + + def get_input( + self, + parameters: PyTree, + rxn_ix: int, + S: NDArray[np.float64], + species_to_dgf_ix: NDArray[np.int16], + ): + return get_irreversible_michaelis_menten_input( + parameters=parameters, + rxn_ix=rxn_ix, + S=S, + species_to_dgf_ix=species_to_dgf_ix, + ci_ix=self.competitive_inhibitor_ix, + ) + + def __call__( + self, + conc: Float[Array, " n"], + imm_input: IrreversibleMichaelisMentenInput, + ) -> Scalar: + """Get flux of a reaction with irreversible Michaelis Menten kinetics.""" # noqa: E501 + numerator = numerator_mm( + substrate_conc=conc[imm_input.ix_substrate], + substrate_kms=imm_input.substrate_kms, + ) + fer = free_enzyme_ratio_imm( + substrate_conc=conc[imm_input.ix_substrate], + substrate_km=imm_input.substrate_kms, + ki=imm_input.ki, + inhibitor_conc=conc[imm_input.ix_ki_species], + substrate_stoichiometry=imm_input.substrate_stoichiometry, + ) + return imm_input.kcat * imm_input.enzyme * numerator * fer + + +class ReversibleMichaelisMenten(RateEquation): """A reaction with reversible Michaelis Menten kinetics.""" - ix_product: Int[Array, " n_product"] - ix_reactants: Int[Array, " n_reactant"] - product_reactant_positions: Int[Array, " n_product"] - product_km_positions: Int[Array, " n_product"] - water_stoichiometry: Scalar - reactant_to_dgf: Int[Array, " n_reactant"] - - def _get_dgf(self, parameters: PyTree): - return parameters.dgf[self.reactant_to_dgf] - - def free_enzyme_ratio(self, conc: ConcArray, parameters: PyTree) -> Scalar: - return free_enzyme_ratio_rmm( - conc_sub=conc[self.ix_substrate], - km_sub=self.get_km(parameters)[self.substrate_reactant_positions], - stoich_sub=self.reactant_stoichiometry[ - self.substrate_reactant_positions - ], - conc_prod=conc[self.ix_product], - km_prod=self.get_km(parameters)[self.product_reactant_positions], - stoich_prod=self.reactant_stoichiometry[ - self.product_reactant_positions - ], - conc_inhibitor=conc[self.ix_ki_species], - ki=self.get_ki(parameters), + competitive_inhibitor_ix: NDArray[np.int16] = eqx.field( + default_factory=lambda: np.array([], dtype=np.int16), static=True + ) + water_stoichiometry: float = eqx.field( + default_factory=lambda: 0.0, static=True + ) + + def get_input( + self, + parameters: PyTree, + rxn_ix: int, + S: NDArray[np.float64], + species_to_dgf_ix: NDArray[np.int16], + ): + return get_reversible_michaelis_menten_input( + parameters=parameters, + rxn_ix=rxn_ix, + S=S, + species_to_dgf_ix=species_to_dgf_ix, + ci_ix=self.competitive_inhibitor_ix, + water_stoichiometry=self.water_stoichiometry, ) - def __call__(self, conc: Float[Array, " n"], parameters: PyTree) -> Scalar: + def __call__( + self, + conc: Float[Array, " n"], + rmm_input: ReversibleMichaelisMentenInput, + ) -> Scalar: """Get flux of a reaction with reversible Michaelis Menten kinetics. :param conc: A 1D array of non-negative numbers representing concentrations of the species that the reaction produces and consumes. """ # noqa: E501 - kcat = self.get_kcat(parameters) - enzyme = self.get_enzyme(parameters) - reversibility = get_reversibility( - conc=conc, - water_stoichiometry=self.water_stoichiometry, - dgf=self._get_dgf(parameters), - temperature=parameters.temperature, - reactant_stoichiometry=self.reactant_stoichiometry, - ix_reactants=self.ix_reactants, + rev = get_reversibility( + reactant_conc=conc[rmm_input.ix_reactant], + reactant_stoichiometry=rmm_input.reactant_stoichiometry, + dgf=rmm_input.dgf, + temperature=rmm_input.temperature, + water_stoichiometry=rmm_input.water_stoichiometry, ) numerator = numerator_mm( - conc=conc, - km=self.get_km(parameters), - ix_substrate=self.ix_substrate, - substrate_km_positions=self.substrate_km_positions, + substrate_conc=conc[rmm_input.ix_substrate], + substrate_kms=rmm_input.substrate_kms, + ) + fer = free_enzyme_ratio_rmm( + substrate_conc=conc[rmm_input.ix_substrate], + product_conc=conc[rmm_input.ix_product], + inhibitor_conc=conc[rmm_input.ix_ki_species], + substrate_kms=rmm_input.substrate_kms, + product_kms=rmm_input.product_kms, + substrate_stoichiometry=rmm_input.substrate_stoichiometry, + product_stoichiometry=rmm_input.product_stoichiometry, + ki=rmm_input.ki, ) - free_enzyme_ratio = self.free_enzyme_ratio(conc, parameters) - return reversibility * kcat * enzyme * numerator * free_enzyme_ratio + return rev * rmm_input.kcat * rmm_input.enzyme * numerator * fer