diff --git a/src/enzax/examples/linear.py b/src/enzax/examples/linear.py index b2dff2c..dc911db 100644 --- a/src/enzax/examples/linear.py +++ b/src/enzax/examples/linear.py @@ -29,28 +29,79 @@ log_drain=jnp.array([]), ) structure = KineticModelStructure( - S=jnp.array([[-1, 0, 0], [1, -1, 0], [0, 1, -1], [0, 0, 1]]), - water_stoichiometry=jnp.array([0, 0, 0]), + S=jnp.array( + [[-1, 0, 0], [1, -1, 0], [0, 1, -1], [0, 0, 1]], dtype=jnp.float64 + ), balanced_species=jnp.array([1, 2]), - rate_to_enzyme_ix=[[0], [1], [2]], - rate_to_km_ixs=jnp.array([[0, 1], [2, 3], [4, 5]]), - species_to_metabolite_ix=jnp.array([0, 0, 1, 1]), - rate_to_subunits=jnp.array([1, 1, 1]), - rate_to_tc_ix=[[0], [1], []], - rate_to_dc_ixs_activation=[[0], [], []], - rate_to_dc_ixs_inhibition=[[], [1], []], - rate_to_drain_ix=[[], [], []], - drain_sign=[], - dc_to_species_ix=jnp.array([2, 1]), - ki_to_species_ix=jnp.array([1]), - rate_to_ki_ixs=[[], [0], []], + unbalanced_species=jnp.array([0, 3]), ) + unparameterised_model = UnparameterisedKineticModel( structure, [ - AllostericReversibleMichaelisMenten, - AllostericReversibleMichaelisMenten, - ReversibleMichaelisMenten, + AllostericReversibleMichaelisMenten( + kcat_ix=0, + enzyme_ix=0, + km_ix=jnp.array([0, 1], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([-1, 1], dtype=jnp.int16), + reactant_to_dgf=jnp.array([0, 0], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0], dtype=jnp.int16), + ix_substrate=jnp.array([0], dtype=jnp.int16), + ix_product=jnp.array([1], dtype=jnp.int16), + ix_reactants=jnp.array([0, 1], dtype=jnp.int16), + product_reactant_positions=jnp.array([1], dtype=jnp.int16), + product_km_positions=jnp.array([1], dtype=jnp.int16), + water_stoichiometry=jnp.array(0.0), + tc_ix=0, + ix_dc_inhibition=jnp.array([], dtype=jnp.int16), + ix_dc_activation=jnp.array([0], dtype=jnp.int16), + species_activation=jnp.array([2], dtype=jnp.int16), + species_inhibition=jnp.array([], dtype=jnp.int16), + subunits=1, + ), + AllostericReversibleMichaelisMenten( + kcat_ix=1, + enzyme_ix=1, + km_ix=jnp.array([2, 3], dtype=jnp.int16), + ki_ix=jnp.array([0]), + reactant_stoichiometry=jnp.array([-1, 1], dtype=jnp.int16), + reactant_to_dgf=jnp.array([0, 1], dtype=jnp.int16), + ix_ki_species=jnp.array([1]), + substrate_km_positions=jnp.array([0], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0], dtype=jnp.int16), + ix_substrate=jnp.array([1], dtype=jnp.int16), + ix_product=jnp.array([2], dtype=jnp.int16), + ix_reactants=jnp.array([1, 2], dtype=jnp.int16), + product_reactant_positions=jnp.array([1], dtype=jnp.int16), + product_km_positions=jnp.array([1], dtype=jnp.int16), + water_stoichiometry=jnp.array(0.0), + tc_ix=1, + ix_dc_inhibition=jnp.array([1], dtype=jnp.int16), + ix_dc_activation=jnp.array([], dtype=jnp.int16), + species_activation=jnp.array([], dtype=jnp.int16), + species_inhibition=jnp.array([1], dtype=jnp.int16), + subunits=1, + ), + ReversibleMichaelisMenten( + kcat_ix=2, + enzyme_ix=2, + km_ix=jnp.array([4, 5], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + ix_substrate=jnp.array([2], dtype=jnp.int16), + ix_product=jnp.array([3], dtype=jnp.int16), + ix_reactants=jnp.array([2, 3], dtype=jnp.int16), + reactant_to_dgf=jnp.array([1, 1], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([-1, 1], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0], dtype=jnp.int16), + product_reactant_positions=jnp.array([1], dtype=jnp.int16), + product_km_positions=jnp.array([1], dtype=jnp.int16), + water_stoichiometry=jnp.array(0.0), + ), ], ) model = KineticModel(parameters, unparameterised_model) diff --git a/src/enzax/examples/methionine.py b/src/enzax/examples/methionine.py index 09635a7..1334b54 100644 --- a/src/enzax/examples/methionine.py +++ b/src/enzax/examples/methionine.py @@ -145,8 +145,8 @@ log_dissociation_constant=jnp.log( jnp.array( [ - 0.000316641, # amet MAT3 0.00059999, # met-L MAT3 + 0.000316641, # amet MAT3 1.98e-05, # amet GNMT1 0.000228576, # mlthf GNMT1 9.30e-05, # amet CBS1 @@ -178,22 +178,8 @@ [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], # nadp [0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0], # nadph [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], # cyst-L - ] - ), - water_stoichiometry=jnp.array( - [ - 0.0, # met-L source - -1.0, # MAT1 METAT - -1.0, # MAT3 METAT - 0.0, # METH Gen - 0.0, # GNMT1 - -1.0, # AHC - 0.0, # MS - 0.0, # BHMT - 1.0, # CBS - 0.0, # MTHFR - 0.0, # PROT - ] + ], + dtype=jnp.float64, ), balanced_species=jnp.array( [ @@ -204,147 +190,175 @@ 11, # 5mthf ] ), - rate_to_enzyme_ix=[ - [], # met-L source - [0], # MAT1 METAT - [1], # MAT3 METAT - [2], # METH Gen - [3], # GNMT1 - [4], # AHC - [5], # MS - [6], # BHMT - [7], # CBS - [8], # MTHFR - [9], # PROT - ], - rate_to_km_ixs=[ - [], # met-L source - [0, 1], # MAT1 METAT - [2, 3], # MAT3 METAT - [4], # METH Gen - [5, 6], # GNMT1 - [7, 8, 9], # AHC - [10, 11], # MS - [12, 13], # BHMT - [14, 15], # CBS - [16, 17], # MTHFR - [18], # PROT - ], - species_to_metabolite_ix=jnp.arange(19), - rate_to_subunits=jnp.array( + unbalanced_species=jnp.array( [ - 0, # met-L source - 1, # MAT1 METAT - 2, # MAT3 METAT - 1, # METH Gen - 4, # GNMT1 - 1, # AHC - 1, # MS - 1, # BHMT - 2, # CBS - 2, # MTHFR - 1, # PROT + 1, # atp + 2, # pi + 3, # ppi + 6, # gly + 7, # sarcs + 9, # adn + 10, # thf + 12, # mlthf + 13, # glyb + 14, # dmgly + 15, # ser-L + 16, # nadp + 17, # nadph + 18, # cyst-L ] ), - rate_to_tc_ix=[ - [], # met-L source - [], # MAT1 METAT - [0], # MAT3 METAT - [], # METH Gen - [1], # GNMT1 - [], # AHC - [], # MS - [], # BHMT - [2], # CBS - [3], # MTHFR - [], # PROT - ], - rate_to_dc_ixs_activation=[ - [], # met-L source - [], # MAT1 METAT - [0, 1], # MAT3 METAT - [], # METH Gen - [2], # GNMT1 - [], # AHC - [], # MS - [], # BHMT - [4], # CBS - [6], # MTHFR - [], # PROT - ], - rate_to_dc_ixs_inhibition=[ - [], # met-L source - [], # MAT1 METAT - [], # MAT3 METAT - [], # METH Gen - [3], # GNMT1 - [], # AHC - [], # MS - [], # BHMT - [], # CBS - [5], # MTHFR - [], # PROT - ], - rate_to_drain_ix=[ - [0], # met-L source - [], # MAT1 METAT - [], # MAT3 METAT - [], # METH Gen - [], # GNMT1 - [], # AHC - [], # MS - [], # BHMT - [], # CBS - [], # MTHFR - [], # PROT - ], - drain_sign=jnp.array([1.0]), - dc_to_species_ix=jnp.array( - [ - 4, # amet -> MAT3 - 0, # met-L -> MAT3 - 4, # amet -> GNMT1 - 12, # mlthf -| GNMT1 - 4, # amet -> CBS1 - 4, # amet -| MTHFR1 - 5, # ahcys -> MTHFR1 - ] - ), - ki_to_species_ix=jnp.array( - [ - 4, # amet -| MAT1 - 5, # ahcys -| METH-GEN - 5, # ahcys -| GNMT1 - ] - ), - rate_to_ki_ixs=[ - [], # met-L source - [0], # MAT1 METAT - [], # MAT3 METAT - [1], # METH Gen - [2], # GNMT1 - [], # AHC - [], # MS - [], # BHMT - [], # CBS - [], # MTHFR - [], # PROT - ], ) unparameterised_model = UnparameterisedKineticModel( structure, [ - Drain, # met-L source - IrreversibleMichaelisMenten, # MAT1 - AllostericIrreversibleMichaelisMenten, # MAT3 - IrreversibleMichaelisMenten, # METH - AllostericIrreversibleMichaelisMenten, # GNMT1 - ReversibleMichaelisMenten, # AHC - IrreversibleMichaelisMenten, # MS - IrreversibleMichaelisMenten, # BHMT - AllostericIrreversibleMichaelisMenten, # CBS - AllostericIrreversibleMichaelisMenten, # MTHFR - IrreversibleMichaelisMenten, # PROT + Drain(sign=jnp.array(1.0), drain_ix=0), # met-L source + IrreversibleMichaelisMenten( # MAT1 + kcat_ix=0, + enzyme_ix=0, + km_ix=jnp.array([0, 1], dtype=jnp.int16), + ki_ix=jnp.array([0], dtype=jnp.int16), + reactant_stoichiometry=jnp.array( + [-1.0, -1.0, 1.0, 1.0, 1.0], dtype=jnp.int16 + ), + ix_substrate=jnp.array([0, 1], dtype=jnp.int16), + ix_ki_species=jnp.array([4], dtype=jnp.int16), + substrate_km_positions=jnp.array([0, 1], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0, 1], dtype=jnp.int16), + ), + AllostericIrreversibleMichaelisMenten( # MAT3 + kcat_ix=1, + enzyme_ix=1, + km_ix=jnp.array([2, 3], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array( + [-1.0, -1.0, 1.0, 1.0, 1.0], dtype=jnp.int16 + ), + ix_substrate=jnp.array([0, 1], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0, 1], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0, 1], dtype=jnp.int16), + subunits=2, + tc_ix=0, + ix_dc_inhibition=jnp.array([], dtype=jnp.int16), + ix_dc_activation=jnp.array([0, 1], dtype=jnp.int16), + species_inhibition=jnp.array([], dtype=jnp.int16), + species_activation=jnp.array([0, 4], dtype=jnp.int16), + ), + IrreversibleMichaelisMenten( # METH + kcat_ix=2, + enzyme_ix=2, + km_ix=jnp.array([4], dtype=jnp.int16), + ki_ix=jnp.array([1], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([-1, 1], dtype=jnp.int16), + ix_substrate=jnp.array([4], dtype=jnp.int16), + ix_ki_species=jnp.array([5], dtype=jnp.int16), + substrate_km_positions=jnp.array([0], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0], dtype=jnp.int16), + ), + AllostericIrreversibleMichaelisMenten( # GNMT1 + kcat_ix=3, + enzyme_ix=3, + km_ix=jnp.array([5, 6], dtype=jnp.int16), + ki_ix=jnp.array([2], dtype=jnp.int16), + reactant_stoichiometry=jnp.array( + [-1.0, 1.0, -1.0, 1.0], dtype=jnp.int16 + ), + ix_substrate=jnp.array([4, 6], dtype=jnp.int16), + ix_ki_species=jnp.array([5], dtype=jnp.int16), + substrate_km_positions=jnp.array([0, 1], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0, 2], dtype=jnp.int16), + subunits=4, + tc_ix=1, + ix_dc_activation=jnp.array([2], dtype=jnp.int16), + ix_dc_inhibition=jnp.array([3], dtype=jnp.int16), + species_inhibition=jnp.array([12], dtype=jnp.int16), + species_activation=jnp.array([4], dtype=jnp.int16), + ), + ReversibleMichaelisMenten( # AHC + kcat_ix=4, + enzyme_ix=4, + km_ix=jnp.array([7, 8, 9], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([-1.0, 1.0, 1.0], dtype=jnp.int16), + ix_substrate=jnp.array([5], dtype=jnp.int16), + ix_product=jnp.array([8, 9], dtype=jnp.int16), + ix_reactants=jnp.array([5, 8, 9], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0], dtype=jnp.int16), + product_km_positions=jnp.array([1, 2], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0], dtype=jnp.int16), + product_reactant_positions=jnp.array([1, 2], dtype=jnp.int16), + water_stoichiometry=jnp.array(-1.0), + reactant_to_dgf=jnp.array([5, 8, 9], dtype=jnp.int16), + ), + IrreversibleMichaelisMenten( # MS + kcat_ix=5, + enzyme_ix=5, + km_ix=jnp.array([10, 11], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([1, -1, 1, -1], dtype=jnp.int16), + ix_substrate=jnp.array([8, 11], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0, 1], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([1, 3], dtype=jnp.int16), + ), + IrreversibleMichaelisMenten( # BHMT + kcat_ix=6, + enzyme_ix=6, + km_ix=jnp.array([12, 13], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([1, -1, -1, 1], dtype=jnp.int16), + ix_substrate=jnp.array([8, 13], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0, 1], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([1, 2], dtype=jnp.int16), + ), + AllostericIrreversibleMichaelisMenten( # CBS1 + kcat_ix=7, + enzyme_ix=7, + km_ix=jnp.array([14, 15], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([-1, -1, 1], dtype=jnp.int16), + ix_substrate=jnp.array([8, 15], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0, 1], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0, 1], dtype=jnp.int16), + subunits=2, + tc_ix=2, + ix_dc_activation=jnp.array([4], dtype=jnp.int16), + ix_dc_inhibition=jnp.array([], dtype=jnp.int16), + species_inhibition=jnp.array([4], dtype=jnp.int16), + species_activation=jnp.array([], dtype=jnp.int16), + ), + AllostericIrreversibleMichaelisMenten( # MTHFR + kcat_ix=8, + enzyme_ix=8, + km_ix=jnp.array([16, 17], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([1, -1, 1, -1], dtype=jnp.int16), + ix_substrate=jnp.array([12, 17], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0, 1], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([1, 3], dtype=jnp.int16), + subunits=2, + tc_ix=3, + ix_dc_activation=jnp.array([6], dtype=jnp.int16), + ix_dc_inhibition=jnp.array([5], dtype=jnp.int16), + species_inhibition=jnp.array([4], dtype=jnp.int16), + species_activation=jnp.array([5], dtype=jnp.int16), + ), + IrreversibleMichaelisMenten( # PROT + kcat_ix=9, + enzyme_ix=9, + km_ix=jnp.array([18], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([-1.0], dtype=jnp.int16), + ix_substrate=jnp.array([0], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0], dtype=jnp.int16), + ), ], ) model = KineticModel(parameters, unparameterised_model) diff --git a/src/enzax/kinetic_model.py b/src/enzax/kinetic_model.py index 8c62689..fa306b1 100644 --- a/src/enzax/kinetic_model.py +++ b/src/enzax/kinetic_model.py @@ -1,9 +1,10 @@ -from abc import abstractmethod, ABC import equinox as eqx import jax.numpy as jnp -from jaxtyping import Array, Float, Int, Scalar, ScalarLike, jaxtyped +from jaxtyping import Array, Float, Int, PyTree, Scalar, ScalarLike, jaxtyped from typeguard import typechecked +from enzax.rate_equations import RateEquation + def ragged_jax_index(lol: list[list[int]]) -> list[Int[Array, " _"]]: """Convert a list of integer lists into a list of jax int arrays. @@ -39,169 +40,28 @@ class KineticModelStructure(eqx.Module): """Structural information about a kinetic model.""" S: Float[Array, " s r"] = eqx.field(static=True) - water_stoichiometry: Float[Array, " r"] = eqx.field(static=True) balanced_species: Int[Array, " n_balanced"] = eqx.field(static=True) - rate_to_enzyme_ix: list[Int[Array, " _"]] = eqx.field(static=True) - rate_to_km_ixs: list[Int[Array, " _"]] = eqx.field(static=True) - rate_to_ki_ixs: list[Int[Array, " _"]] = eqx.field(static=True) - rate_to_tc_ix: list[Int[Array, " _"]] = eqx.field(static=True) - rate_to_drain_ix: list[Int[Array, " _"]] = eqx.field(static=True) - drain_sign: Float[Array, " n_drain"] = eqx.field(static=True) - rate_to_dc_ixs_inhibition: list[Int[Array, " _"]] = eqx.field(static=True) - rate_to_dc_ixs_activation: list[Int[Array, " _"]] = eqx.field(static=True) - rate_to_subunits: Int[Array, " r"] = eqx.field(static=True) - species_to_metabolite_ix: Int[Array, " s"] = eqx.field(static=True) - ki_to_species_ix: Int[Array, " n_ki"] = eqx.field(static=True) - dc_to_species_ix: Int[Array, " n_dc"] = eqx.field(static=True) unbalanced_species: Int[Array, " n_unbalanced"] = eqx.field(static=True) - rate_to_reactants: list[Int[Array, " _"]] = eqx.field(static=True) - rate_to_substrates: list[Int[Array, " _"]] = eqx.field(static=True) - rate_to_products: list[Int[Array, " _"]] = eqx.field(static=True) - rate_to_substrate_reactant_positions: list[Int[Array, " _"]] = eqx.field( - static=True - ) - rate_to_product_reactant_positions: list[Int[Array, " _"]] = eqx.field( - static=True - ) - rate_to_stoichs: list[Float[Array, " _"]] = eqx.field(static=True) - - def __init__( - self, - S, - water_stoichiometry, - balanced_species, - rate_to_enzyme_ix, - rate_to_km_ixs, - rate_to_ki_ixs, - rate_to_tc_ix, - rate_to_drain_ix, - drain_sign, - rate_to_dc_ixs_inhibition, - rate_to_dc_ixs_activation, - rate_to_subunits, - species_to_metabolite_ix, - ki_to_species_ix, - dc_to_species_ix, - ): - self.S = jnp.array(S, dtype=jnp.float64) - self.water_stoichiometry = jnp.array( - water_stoichiometry, dtype=jnp.float64 - ) - self.balanced_species = jnp.array(balanced_species, dtype=jnp.int16) - self.rate_to_enzyme_ix = ragged_jax_index(rate_to_enzyme_ix) - self.rate_to_km_ixs = ragged_jax_index(rate_to_km_ixs) - self.rate_to_ki_ixs = ragged_jax_index(rate_to_ki_ixs) - self.rate_to_tc_ix = ragged_jax_index(rate_to_tc_ix) - self.rate_to_drain_ix = ragged_jax_index(rate_to_drain_ix) - self.drain_sign = jnp.array(drain_sign) - self.rate_to_dc_ixs_inhibition = ragged_jax_index( - rate_to_dc_ixs_inhibition - ) - self.rate_to_dc_ixs_activation = ragged_jax_index( - rate_to_dc_ixs_activation - ) - self.rate_to_subunits = rate_to_subunits - self.species_to_metabolite_ix = species_to_metabolite_ix - self.ki_to_species_ix = ki_to_species_ix - self.dc_to_species_ix = dc_to_species_ix - self.unbalanced_species = jnp.array( - [ - i - for i in range(self.S.shape[0]) - if i not in self.balanced_species - ] - ) - self.rate_to_stoichs = [ - jnp.array( - [coeff for coeff in coeffs if coeff != 0], dtype=jnp.float64 - ) - for coeffs in self.S.T - ] - self.rate_to_reactants = [ - jnp.array( - [i for i, coeff in enumerate(coeffs) if coeff != 0], - dtype=jnp.int16, - ) - for coeffs in self.S.T - ] - self.rate_to_substrates = [ - jnp.array( - [i for i, coeff in enumerate(coeffs) if coeff < 0], - dtype=jnp.int16, - ) - for coeffs in self.S.T - ] - self.rate_to_products = [ - jnp.array( - [i for i, coeff in enumerate(coeffs) if coeff > 0], - dtype=jnp.int16, - ) - for coeffs in self.S.T - ] - self.rate_to_substrate_reactant_positions = [ - jnp.array( - [i for i, coeff in enumerate(coeffs) if coeff < 0], - dtype=jnp.int16, - ) - for coeffs in self.rate_to_stoichs - ] - self.rate_to_product_reactant_positions = [ - jnp.array( - [i for i, coeff in enumerate(coeffs) if coeff > 0], - dtype=jnp.int16, - ) - for coeffs in self.rate_to_stoichs - ] - - -class RateEquation(ABC, eqx.Module): - """Class representing an abstract rate equation.""" - - @abstractmethod - def __init__( - self, - parameters: KineticModelParameters, - structure: KineticModelStructure, - ix: int, - ): - """Signature for the __init__ method of a rate equation. - - A rate equation is initialised from a set of parameters, a structure and a positive integer index. - """ - ... - - @abstractmethod - def __call__(self, conc: Float[Array, " n"]) -> Scalar: - """Signature for the __call__ method of a rate equation. - - A rate equation takes in a one dimensional array of positive float-valued concentrations and returns a scalar flux. - """ - ... class UnparameterisedKineticModel(eqx.Module): """A kinetic model without parameter values.""" structure: KineticModelStructure - rate_equation_classes: list[type[RateEquation]] + rate_equations: list[RateEquation] class KineticModel(eqx.Module): """A parameterised kinetic model.""" - parameters: KineticModelParameters + parameters: PyTree structure: KineticModelStructure rate_equations: list[RateEquation] def __init__(self, parameters, unparameterised_model): self.parameters = parameters self.structure = unparameterised_model.structure - self.rate_equations = [ - cls(self.parameters, self.structure, ix) - for ix, cls in enumerate( - unparameterised_model.rate_equation_classes - ) - ] + self.rate_equations = unparameterised_model.rate_equations def flux( self, conc_balanced: Float[Array, " n_balanced"] @@ -218,7 +78,7 @@ def flux( conc = conc.at[self.structure.unbalanced_species].set( jnp.exp(self.parameters.log_conc_unbalanced) ) - t = [f(conc) for ix, f in enumerate(self.rate_equations)] + t = [f(conc, self.parameters) for f in self.rate_equations] out = jnp.array(t) return out diff --git a/src/enzax/rate_equations.py b/src/enzax/rate_equations.py index c95e768..a6d913c 100644 --- a/src/enzax/rate_equations.py +++ b/src/enzax/rate_equations.py @@ -1,273 +1,342 @@ """Module containing rate equations for enzyme-catalysed reactions.""" -from abc import abstractmethod +from abc import ABC, abstractmethod +from equinox import Module from jax import numpy as jnp -from jaxtyping import Array, Float, Int, Scalar +from jaxtyping import Array, Float, Int, PyTree, Scalar -from enzax.kinetic_model import ( - KineticModelParameters, - KineticModelStructure, - RateEquation, -) + +ConcArray = Float[Array, " n"] + + +class RateEquation(Module, ABC): + @abstractmethod + def __call__(self, conc: ConcArray, parameters: PyTree) -> Scalar: ... + + +def drain_function(sign: Scalar, log_v: Scalar) -> Scalar: + return sign * jnp.exp(log_v) class Drain(RateEquation): - log_v: Scalar sign: Scalar + drain_ix: int - def __init__( - self, - parameters: KineticModelParameters, - structure: KineticModelStructure, - ix: int, - ): - self.log_v = parameters.log_drain[ - jnp.array(structure.rate_to_drain_ix[ix][0]) - ] - self.sign = structure.drain_sign[ - jnp.array(structure.rate_to_drain_ix[ix])[0] - ] - - def __call__(self, conc: Float[Array, " n"]) -> Scalar: - """Get flux of a drain reaction. + def __call__(self, conc: ConcArray, parameters: PyTree) -> Scalar: + log_v = parameters.log_drain[self.drain_ix] + return drain_function(self.sign, log_v) - :param conc: A 1D array of non-negative numbers representing concentrations of the species that the reaction produces and consumes. - """ - return self.sign * jnp.exp(self.log_v) +def numerator_mm( + conc: ConcArray, + km: Float[Array, " n"], + ix_substrate: Int[Array, " n_substrate"], + substrate_km_positions: 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 + """ + return jnp.prod((conc[ix_substrate] / km[substrate_km_positions])) class MichaelisMenten(RateEquation): - """Abstract base class for Michaelis Menten rate equations. + """Base class for Michaelis Menten rate equations. - Subclasses need to implement the method free_enzyme_ratio (and also the __call__ method). + Subclasses need to implement the __call__ method. """ - dgf: Float[Array, " n"] - log_km: Float[Array, " n"] - log_enzyme: Scalar - log_kcat: Scalar - log_ki: Float[Array, " n_ki"] - temperature: Scalar - stoich: Float[Array, " n"] + 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"] - water_stoichiometry: Scalar substrate_km_positions: Int[Array, " n_substrate"] substrate_reactant_positions: Int[Array, " n_substrate"] - def __init__( - self, - parameters: KineticModelParameters, - structure: KineticModelStructure, - ix: int, - ): - ix_dgf = structure.species_to_metabolite_ix[ - structure.rate_to_reactants[ix] - ] - self.dgf = parameters.dgf[ix_dgf] - self.log_km = parameters.log_km[structure.rate_to_km_ixs[ix]] - self.log_enzyme = parameters.log_enzyme[ - structure.rate_to_enzyme_ix[ix][0] - ] - self.log_kcat = parameters.log_kcat[structure.rate_to_enzyme_ix[ix][0]] - self.log_ki = parameters.log_ki[structure.rate_to_ki_ixs[ix]] - self.temperature = parameters.temperature - self.stoich = structure.rate_to_stoichs[ix] - self.ix_substrate = structure.rate_to_substrates[ix] - self.substrate_km_positions = jnp.arange(len(self.ix_substrate)) - self.substrate_reactant_positions = ( - structure.rate_to_substrate_reactant_positions[ix] + 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]) + + +def free_enzyme_ratio_imm( + conc: ConcArray, + km: Float[Array, " n"], + ki: Float[Array, " n_ki"], + ix_substrate: Int[Array, " n_substrate"], + substrate_km_positions: Int[Array, " n_substrate"], + substrate_reactant_positions: Int[Array, " n_substrate"], + ix_ki_species: Int[Array, " n_ki"], + reactant_stoichiometry: Float[Array, " n"], +) -> Scalar: + """Free enzyme ratio for irreversible Michaelis Menten reactions.""" + return 1.0 / ( + jnp.prod( + ((conc[ix_substrate] / km[substrate_km_positions]) + 1) + ** jnp.abs(reactant_stoichiometry[substrate_reactant_positions]) ) - self.ix_ki_species = structure.ki_to_species_ix[ - structure.rate_to_ki_ixs[ix] - ] - self.water_stoichiometry = structure.water_stoichiometry[ix] - - @property - def km(self): - return jnp.exp(self.log_km) - - @property - def kcat(self): - return jnp.exp(self.log_kcat) - - @property - def ki(self): - return jnp.exp(self.log_ki) - - @property - def enzyme(self): - return jnp.exp(self.log_enzyme) - - def numerator(self, conc: Float[Array, " n"]) -> 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 - """ - return jnp.prod( - (conc[self.ix_substrate] / self.km[self.substrate_km_positions]) - ) - - @abstractmethod - def free_enzyme_ratio(self, conc: Float[Array, " n"]) -> Scalar: ... + + jnp.sum(conc[ix_ki_species] / ki) + ) class IrreversibleMichaelisMenten(MichaelisMenten): - def free_enzyme_ratio(self, conc: Float[Array, " n"]) -> Scalar: - return 1.0 / ( - jnp.prod( - ( - ( - conc[self.ix_substrate] - / self.km[self.substrate_km_positions] - ) - + 1 - ) - ** jnp.abs(self.stoich[self.substrate_reactant_positions]) - ) - + jnp.sum(conc[self.ix_ki_species] / self.ki) + def __call__(self, conc: Float[Array, " n"], parameters: PyTree) -> Scalar: + """Get flux of a reaction with irreversible Michaelis Menten kinetics.""" + kcat = self.get_kcat(parameters) + enzyme = self.get_enzyme(parameters) + km = self.get_km(parameters) + ki = self.get_ki(parameters) + numerator = numerator_mm( + conc=conc, + km=km, + ix_substrate=self.ix_substrate, + substrate_km_positions=self.substrate_km_positions, ) + free_enzyme_ratio = free_enzyme_ratio_imm( + conc=conc, + km=km, + ki=ki, + ix_substrate=self.ix_substrate, + substrate_km_positions=self.substrate_km_positions, + substrate_reactant_positions=self.substrate_reactant_positions, + ix_ki_species=self.ix_ki_species, + reactant_stoichiometry=self.reactant_stoichiometry, + ) + return kcat * enzyme * numerator * free_enzyme_ratio - def __call__(self, conc: Float[Array, " n"]) -> Scalar: - """Get flux of a reaction with irreversible Michaelis Menten kinetics. - :param conc: A 1D array of non-negative numbers representing concentrations of the species that the reaction produces and consumes. +def get_reversibility( + conc: Float[Array, " n"], + water_stoichiometry: Scalar, + dgf: Float[Array, " n_reactant"], + temperature: Scalar, + reactant_stoichiometry: Float[Array, " n_reactant"], + ix_reactants: Int[Array, " n_reactant"], +) -> Scalar: + """Get the reversibility of a reaction. - """ - saturation: Scalar = self.numerator(conc) * self.free_enzyme_ratio(conc) - out = self.kcat * self.enzyme * saturation - return out + Hard coded water dgf is taken from . + + """ + 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)) + return out + + +def free_enzyme_ratio_rmm( + conc: ConcArray, + km: Float[Array, " n_reactant"], + ki: Float[Array, " n_ki"], + reactant_stoichiometry: Float[Array, " n_reactant"], + ix_substrate: Int[Array, " n_substrate"], + ix_product: Int[Array, " n_product"], + substrate_km_positions: Int[Array, " n_substrate"], + product_km_positions: Int[Array, " n_product"], + substrate_reactant_positions: Int[Array, " n_substrate"], + product_reactant_positions: Int[Array, " n_product"], + ix_ki_species: Int[Array, " n_ki"], +) -> Scalar: + return 1.0 / ( + -1.0 + + jnp.prod( + ((conc[ix_substrate] / km[substrate_km_positions]) + 1.0) + ** jnp.abs(reactant_stoichiometry[substrate_reactant_positions]) + ) + + jnp.prod( + ((conc[ix_product] / km[product_km_positions]) + 1.0) + ** jnp.abs(reactant_stoichiometry[product_reactant_positions]) + ) + + jnp.sum(conc[ix_ki_species] / ki) + ) class ReversibleMichaelisMenten(MichaelisMenten): 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 __init__( - self, - parameters: KineticModelParameters, - structure: KineticModelStructure, - ix: int, - ): - super().__init__(parameters, structure, ix) - self.ix_product = structure.rate_to_products[ix] - self.ix_reactants = structure.rate_to_reactants[ix] - self.product_reactant_positions = ( - structure.rate_to_product_reactant_positions[ix] - ) - - def reversibility( - self, conc: Float[Array, " n"], water_stoichiometry: Scalar - ) -> Scalar: - """Get the reversibility of a reaction. - - Hard coded water dgf is taken from . - - """ - RT = self.temperature * 0.008314 - dgf_water = -150.9 - dgr = self.stoich @ self.dgf + water_stoichiometry * dgf_water - quotient = self.stoich @ jnp.log(conc[self.ix_reactants]) - out = 1.0 - jnp.exp(((dgr + RT * quotient) / RT)) - return out - - def free_enzyme_ratio(self, conc: Float[Array, " n"]) -> Scalar: - return 1.0 / ( - -1.0 - + jnp.prod( - ( - ( - conc[self.ix_substrate] - / self.km[self.substrate_reactant_positions] - ) - + 1.0 - ) - ** jnp.abs(self.stoich[self.substrate_reactant_positions]) - ) - + jnp.prod( - ( - ( - conc[self.ix_product] - / self.km[self.product_reactant_positions] - ) - + 1.0 - ) - ** jnp.abs(self.stoich[self.product_reactant_positions]) - ) - + jnp.sum(conc[self.ix_ki_species] / self.ki) - ) + def get_dgf(self, parameters: PyTree): + return parameters.dgf[self.reactant_to_dgf] - def __call__(self, conc: Float[Array, " n"]) -> Scalar: + def __call__(self, conc: Float[Array, " n"], parameters: PyTree) -> 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. """ - saturation: Scalar = self.numerator(conc) * self.free_enzyme_ratio(conc) - out = ( - self.reversibility(conc, self.water_stoichiometry) - * self.kcat - * self.enzyme - * saturation + kcat = self.get_kcat(parameters) + enzyme = self.get_enzyme(parameters) + km = self.get_km(parameters) + ki = self.get_ki(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, + ) + numerator = numerator_mm( + conc=conc, + km=km, + ix_substrate=self.ix_substrate, + substrate_km_positions=self.substrate_km_positions, ) - return out + free_enzyme_ratio = free_enzyme_ratio_rmm( + conc=conc, + km=km, + ki=ki, + reactant_stoichiometry=self.reactant_stoichiometry, + ix_substrate=self.ix_substrate, + ix_product=self.ix_product, + substrate_km_positions=self.substrate_km_positions, + product_km_positions=self.product_km_positions, + substrate_reactant_positions=self.substrate_reactant_positions, + product_reactant_positions=self.product_reactant_positions, + ix_ki_species=self.ix_ki_species, + ) + return reversibility * kcat * enzyme * numerator * free_enzyme_ratio + + +def get_allosteric_effect( + conc: Float[Array, " n_reactant"], + free_enzyme_ratio: Scalar, + tc: Scalar, + dc_inhibition: Float[Array, " n_inhibition"], + dc_activation: Float[Array, " n_activation"], + species_inhibition: Int[Array, " n_inhibition"], + species_activation: Int[Array, " n_activation"], + subunits: int, +) -> Scalar: + qnum = 1 + jnp.sum(conc[species_inhibition] / dc_inhibition) + qdenom = 1 + jnp.sum(conc[species_activation] / dc_activation) + out = 1.0 / (1 + tc * (free_enzyme_ratio * qnum / qdenom) ** subunits) + return out class AllostericRateLaw(MichaelisMenten): - """Mixin class providing the method allosteric_effect.""" + """Mixin class for allosteric rate laws.""" 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"] - tc: Scalar - dc_activation: Float[Array, " n_activation"] - dc_inhibition: Float[Array, " n_inhibition"] - def __init__(self, parameters, structure, ix): - ix_dc_activation = jnp.array( - structure.rate_to_dc_ixs_activation[ix], dtype=jnp.int16 - ) - ix_dc_inhibition = jnp.array( - structure.rate_to_dc_ixs_inhibition[ix], dtype=jnp.int16 - ) - super().__init__(parameters, structure, ix) - self.subunits = structure.rate_to_subunits[ix] # type: ignore - self.species_activation = structure.dc_to_species_ix[ix_dc_activation] - self.species_inhibition = structure.dc_to_species_ix[ix_dc_inhibition] - self.tc = jnp.exp( - parameters.log_transfer_constant[structure.rate_to_tc_ix[ix][0]] - ) - self.dc_activation = jnp.exp( - parameters.log_dissociation_constant[ix_dc_activation] - ) - self.dc_inhibition = jnp.exp( - parameters.log_dissociation_constant[ix_dc_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 allosteric_effect(self, conc: Float[Array, " n"]) -> Scalar: - qnum = 1 + jnp.sum(conc[self.species_inhibition] / self.dc_inhibition) - qdenom = 1 + jnp.sum(conc[self.species_activation] / self.dc_activation) - out = 1.0 / ( - 1 - + self.tc - * (self.free_enzyme_ratio(conc) * qnum / qdenom) ** self.subunits + def get_dc_inhibition(self, parameters: PyTree) -> Scalar: + return jnp.exp( + parameters.log_dissociation_constant[self.ix_dc_inhibition] ) - return out class AllostericIrreversibleMichaelisMenten( AllostericRateLaw, IrreversibleMichaelisMenten ): - def __call__(self, conc: Float[Array, " n"]) -> Scalar: - out = super().__call__(conc) * self.allosteric_effect(conc) - return out + def __call__(self, conc: Float[Array, " n"], parameters: PyTree) -> Scalar: + km = self.get_km(parameters) + ki = self.get_ki(parameters) + tc = self.get_tc(parameters) + dc_activation = self.get_dc_activation(parameters) + dc_inhibition = self.get_dc_inhibition(parameters) + free_enzyme_ratio = free_enzyme_ratio_imm( + conc=conc, + km=km, + ki=ki, + ix_substrate=self.ix_substrate, + substrate_km_positions=self.substrate_km_positions, + substrate_reactant_positions=self.substrate_reactant_positions, + ix_ki_species=self.ix_ki_species, + reactant_stoichiometry=self.reactant_stoichiometry, + ) + allosteric_effect = get_allosteric_effect( + conc=conc, + free_enzyme_ratio=free_enzyme_ratio, + tc=tc, + dc_inhibition=dc_inhibition, + dc_activation=dc_activation, + species_inhibition=self.species_inhibition, + species_activation=self.species_activation, + subunits=self.subunits, + ) + non_allosteric_rate = super().__call__(conc, parameters) + return non_allosteric_rate * allosteric_effect class AllostericReversibleMichaelisMenten( AllostericRateLaw, ReversibleMichaelisMenten ): - def __call__(self, conc: Float[Array, " n"]) -> Scalar: - return super().__call__(conc) * self.allosteric_effect(conc) + def __call__(self, conc: ConcArray, parameters: PyTree) -> Scalar: + km = self.get_km(parameters) + ki = self.get_ki(parameters) + tc = self.get_tc(parameters) + dc_activation = self.get_dc_activation(parameters) + dc_inhibition = self.get_dc_inhibition(parameters) + free_enzyme_ratio = free_enzyme_ratio_rmm( + conc=conc, + km=km, + ki=ki, + reactant_stoichiometry=self.reactant_stoichiometry, + ix_substrate=self.ix_substrate, + ix_product=self.ix_product, + substrate_km_positions=self.substrate_km_positions, + product_km_positions=self.product_km_positions, + substrate_reactant_positions=self.substrate_reactant_positions, + product_reactant_positions=self.product_reactant_positions, + ix_ki_species=self.ix_ki_species, + ) + allosteric_effect = get_allosteric_effect( + conc=conc, + free_enzyme_ratio=free_enzyme_ratio, + tc=tc, + dc_inhibition=dc_inhibition, + dc_activation=dc_activation, + species_inhibition=self.species_inhibition, + species_activation=self.species_activation, + subunits=self.subunits, + ) + non_allosteric_rate = super().__call__(conc, parameters) + return non_allosteric_rate * allosteric_effect + + +def main(): + _ = IrreversibleMichaelisMenten( + kcat_ix=0, + enzyme_ix=0, + km_ix=jnp.array([0, 1]), + ki_ix=jnp.array([]), + reactant_stoichiometry=jnp.array([-1, 1]), + ix_substrate=jnp.array([0]), + ix_ki_species=jnp.array([]), + substrate_km_positions=jnp.array([0]), + substrate_reactant_positions=jnp.array([0]), + ) + + +if __name__ == "__main__": + main() diff --git a/tests/test_rate_equations.py b/tests/test_rate_equations.py index 73606b6..9816a69 100644 --- a/tests/test_rate_equations.py +++ b/tests/test_rate_equations.py @@ -1,6 +1,5 @@ """Unit tests for rate equations.""" -import pytest from jax import numpy as jnp from enzax.rate_equations import ( AllostericIrreversibleMichaelisMenten, @@ -8,22 +7,110 @@ IrreversibleMichaelisMenten, ReversibleMichaelisMenten, ) -from enzax.examples import linear +from enzax.kinetic_model import KineticModelParameters - -@pytest.mark.parametrize( - ["rate_equation", "expected_rate"], - [ - (IrreversibleMichaelisMenten, 0.01105354), - (ReversibleMichaelisMenten, -0.01188049), - (AllostericIrreversibleMichaelisMenten, 0.00474257), - (AllostericReversibleMichaelisMenten, -0.00535258), - ], +EXAMPLE_CONC = jnp.array([0.5, 0.2, 0.1, 0.3]) +EXAMPLE_PARAMETERS = KineticModelParameters( + log_kcat=jnp.array([-0.1]), + log_enzyme=jnp.log(jnp.array([0.3])), + dgf=jnp.array([-3, -1.0]), + log_km=jnp.array([0.1, -0.2]), + log_ki=jnp.array([1.0]), + log_conc_unbalanced=jnp.log(jnp.array([0.5, 0.3])), + temperature=jnp.array(310.0), + log_transfer_constant=jnp.array([-0.2, 0.3]), + log_dissociation_constant=jnp.array([-0.1, 0.2]), + log_drain=jnp.array([]), ) -def test_rate_equations_linear_one(rate_equation, expected_rate): - """Check output of different rate equations for rate 1 of linear model.""" - conc = jnp.array([0.5, 0.2, 0.1, 0.3]) - ix = 1 - f = rate_equation(linear.parameters, linear.structure, ix) - rate = f(conc[linear.structure.rate_to_reactants[ix]]) + + +def test_irreversible_michaelis_menten(): + expected_rate = 0.08455524 + f = IrreversibleMichaelisMenten( + kcat_ix=0, + enzyme_ix=0, + km_ix=jnp.array([0, 1], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([-1, 1], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0], dtype=jnp.int16), + ix_substrate=jnp.array([0], dtype=jnp.int16), + ) + rate = f(EXAMPLE_CONC, EXAMPLE_PARAMETERS) + assert jnp.isclose(rate, expected_rate) + + +def test_reversible_michaelis_menten(): + expected_rate = 0.04342889 + f = ReversibleMichaelisMenten( + kcat_ix=0, + enzyme_ix=0, + km_ix=jnp.array([0, 1], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([-1, 1], dtype=jnp.int16), + reactant_to_dgf=jnp.array([0, 0], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0], dtype=jnp.int16), + ix_substrate=jnp.array([0], dtype=jnp.int16), + ix_product=jnp.array([1], dtype=jnp.int16), + ix_reactants=jnp.array([0, 1], dtype=jnp.int16), + product_reactant_positions=jnp.array([1], dtype=jnp.int16), + product_km_positions=jnp.array([1], dtype=jnp.int16), + water_stoichiometry=jnp.array(0.0), + ) + rate = f(EXAMPLE_CONC, EXAMPLE_PARAMETERS) + assert jnp.isclose(rate, expected_rate) + + +def test_allosteric_irreversible_michaelis_menten(): + expected_rate = 0.05608589 + f = AllostericIrreversibleMichaelisMenten( + kcat_ix=0, + enzyme_ix=0, + km_ix=jnp.array([0, 1], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([-1, 1], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0], dtype=jnp.int16), + ix_substrate=jnp.array([0], dtype=jnp.int16), + tc_ix=0, + ix_dc_inhibition=jnp.array([], dtype=jnp.int16), + ix_dc_activation=jnp.array([0], dtype=jnp.int16), + species_activation=jnp.array([2], dtype=jnp.int16), + species_inhibition=jnp.array([], dtype=jnp.int16), + subunits=1, + ) + rate = f(EXAMPLE_CONC, EXAMPLE_PARAMETERS) + assert jnp.isclose(rate, expected_rate) + + +def test_allosteric_reversible_michaelis_menten(): + expected_rate = 0.03027414 + f = AllostericReversibleMichaelisMenten( + kcat_ix=0, + enzyme_ix=0, + km_ix=jnp.array([0, 1], dtype=jnp.int16), + ki_ix=jnp.array([], dtype=jnp.int16), + reactant_stoichiometry=jnp.array([-1, 1], dtype=jnp.int16), + reactant_to_dgf=jnp.array([0, 0], dtype=jnp.int16), + ix_ki_species=jnp.array([], dtype=jnp.int16), + substrate_km_positions=jnp.array([0], dtype=jnp.int16), + substrate_reactant_positions=jnp.array([0], dtype=jnp.int16), + ix_substrate=jnp.array([0], dtype=jnp.int16), + ix_product=jnp.array([1], dtype=jnp.int16), + ix_reactants=jnp.array([0, 1], dtype=jnp.int16), + product_reactant_positions=jnp.array([1], dtype=jnp.int16), + product_km_positions=jnp.array([1], dtype=jnp.int16), + water_stoichiometry=jnp.array(0.0), + tc_ix=0, + ix_dc_inhibition=jnp.array([], dtype=jnp.int16), + ix_dc_activation=jnp.array([0], dtype=jnp.int16), + species_activation=jnp.array([2], dtype=jnp.int16), + species_inhibition=jnp.array([], dtype=jnp.int16), + subunits=1, + ) + rate = f(EXAMPLE_CONC, EXAMPLE_PARAMETERS) assert jnp.isclose(rate, expected_rate)