Skip to content

Yeast glycolysis oscillations in a flow reactor

The glycolysis of a yeast extract in a flow reactor, whose NADH oscillates with a period of minutes.

  • Name: glycolysis_oscillations
  • Field: biology/systems_biology
  • Formalism: ordinary differential equations
  • References: Nielsen et al. 1998; Raman 2021

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('glycolysis_oscillations')
entry = describe('glycolysis_oscillations')
trajectory = create_simulation('glycolysis_oscillations').prepare().run().trajectory_

Factory

@register_system(
    field='biology/systems_biology',
    tags=('glycolysis', 'metabolic oscillations', 'nadh', 'biomodels', 'ode'),
    references=('Nielsen et al. 1998', 'Raman 2021'),
    defaults={'until': 300.0, 'dt': 0.5},
)
def glycolysis_oscillations(*, flow: float = 0.011, V4: float = 10.0) -> System:
    """Return the glycolysis of a yeast extract in a flow reactor, whose NADH oscillates with a period of minutes.

    The model of Nielsen et al. (1998), BioModels BIOMD0000000042, has fifteen metabolites, the inhibition of
    phosphofructokinase by `ATP` and its activation by `AMP`, and the exchange with the feed of the reactor at the rate
    `flow`. With the published parameters the reactor settles on sustained oscillations of `NADH` and of all the
    other metabolites. Raman (2021, exercise 5.7) varies `V4`, `K4GAP`, `K4NAD`, `k8f`, and `k8b` by half to study
    their effect on `NAD` and `NADH`.

    Args:
        flow:
            The flow rate of the reactor, per minute.

        V4:
            The largest rate of glyceraldehyde phosphate dehydrogenase, in mM per minute.

    Returns:
        system:
            The system, with the concentrations `extract.<metabolite>`, such as `extract.NADH`, in mM.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import glycolysis_oscillations
        >>> glycolysis_oscillations(flow=0.02).elements[0].get_value('flow')
        0.02
    """
    variables = {
        'flow': DeclaredVariableSpec(flow, ParameterSpec(domain=RealSet.closed_open(0, oo))),
        'V4': DeclaredVariableSpec(V4, ParameterSpec(domain=RealSet.closed_open(0, oo))),
    }
    return System(elements=[GlycolysisElement('extract', variables=variables)])

Components

Elements:

  • GlycolysisElement, from DeclaredElement
class GlycolysisElement(DeclaredElement):
    """Yeast glycolysis in a flow reactor, the fifteen metabolites of the model of Nielsen and coworkers.

    The model of Nielsen, Sørensen, Hynne, and Busse (1998), BioModels BIOMD0000000042, follows glucose `GLC` through
    hexokinase, phosphofructokinase with its inhibition by the ratio of `ATP` to `AMP`, aldolase, glyceraldehyde
    phosphate dehydrogenase, phosphoglycerate kinase, pyruvate kinase, pyruvate decarboxylase, alcohol dehydrogenase,
    adenylate kinase, and a glycerol branch to `P`. The reactor is fed with `ATP`, `ADP`, `NADH`, `NAD`, and glucose
    at fixed concentrations and washed out at the rate `flow`. The concentrations, in mM, and the rate constants, in mM
    and per minute, are declared variables named as in the source, and `variables` replaces any of them.

    Args:
        label:
            The label of the cell extract.

        visual:
            How the element is drawn, or `None` for `default_visual_`.

        variables:
            The declared states and parameters that replace those of the source, by name, or `None` for none.

    Examples:
        >>> from skcomplex.systems.biology.systems_biology import GlycolysisElement
        >>> GlycolysisElement('extract').get_value('V4')
        10.0
    """

    def _declare_variables(self) -> Mapping[str, DeclaredVariableSpec]:
        """Return the concentrations and the rate constants of the source, updated by `variables`."""
        nonnegative = RealSet.closed_open(0, oo)
        declared = {
            name: DeclaredVariableSpec(value, StateSpec(domain=nonnegative)) for name, value in CONCENTRATIONS.items()
        }
        declared |= {
            name: DeclaredVariableSpec(value, ParameterSpec(domain=nonnegative)) for name, value in CONSTANTS.items()
        }
        return {**declared, **(self.variables or {})}

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the rates of the metabolites, the reactions of glycolysis and the exchange with the feed.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the component.
        """
        v1 = s.V1 * s.ATP * s.GLC / ((s.K1GLC + s.GLC) * (s.K1ATP + s.ATP))
        inhibition = s.K2 * (1 + s.k2 * (s.ATP / s.AMP) ** 2)
        v2 = s.V2 * s.ATP * s.F6P**2 / ((inhibition + s.F6P**2) * (s.K2ATP + s.ATP))
        v3 = s.k3f * s.FBP - s.k3b * s.GAP**2
        v4 = s.V4 * s.NAD * s.GAP / ((s.K4GAP + s.GAP) * (s.K4NAD + s.NAD))
        v5 = s.k5f * s.DPG * s.ADP - s.k5b * s.PEP * s.ATP
        v6 = s.V6 * s.ADP * s.PEP / ((s.K6PEP + s.PEP) * (s.K6ADP + s.ADP))
        v7 = s.V7 * s.PYR / (s.K7PYR + s.PYR)
        v8 = s.k8f * s.NADH * s.ACA - s.k8b * s.NAD * s.EtOH
        v9 = s.k9f * s.AMP * s.ATP - s.k9b * s.ADP**2
        v10 = s.k10 * s.F6P
        reactions = {
            'ATP': -v1 - v2 + v5 + v6 - v9,
            'ADP': v1 + v2 - v5 - v6 + 2 * v9,
            'AMP': -v9,
            'GLC': -v1,
            'F6P': v1 - v2 - v10,
            'FBP': v2 - v3,
            'GAP': 2 * v3 - v4,
            'NAD': -v4 + v8,
            'NADH': v4 - v8,
            'DPG': v4 - v5,
            'PEP': v5 - v6,
            'PYR': v6 - v7,
            'ACA': v7 - v8,
            'EtOH': v8,
            'P': v10,
        }
        laws: list[BaseLaw] = []
        for name, rate in reactions.items():
            state = getattr(s, name)
            exchange = s.flow * (FEEDS.get(name, 0.0) - state)
            laws.append(RateLaw(state, rate + exchange))
        return laws

Equations

extract.ATP: diff(ATP_extract(t), t) == k9b_extract*ADP_extract(t)^2 - k9f_extract*AMP_extract(t)*ATP_extract(t) + k5f_extract*ADP_extract(t)*DPG_extract(t) - k5b_extract*ATP_extract(t)*PEP_extract(t) - flow_extract*(ATP_extract(t) - 3.5) - V2_extract*ATP_extract(t)*F6P_extract(t)^2/((K2_extract*(k2_extract*ATP_extract(t)^2/AMP_extract(t)^2 + 1) + F6P_extract(t)^2)*(K2ATP_extract + ATP_extract(t))) - V1_extract*ATP_extract(t)*GLC_extract(t)/((K1ATP_extract + ATP_extract(t))*(K1GLC_extract + GLC_extract(t))) + V6_extract*ADP_extract(t)*PEP_extract(t)/((K6ADP_extract + ADP_extract(t))*(K6PEP_extract + PEP_extract(t)))
extract.ADP: diff(ADP_extract(t), t) == -2*k9b_extract*ADP_extract(t)^2 + 2*k9f_extract*AMP_extract(t)*ATP_extract(t) - k5f_extract*ADP_extract(t)*DPG_extract(t) + k5b_extract*ATP_extract(t)*PEP_extract(t) - flow_extract*(ADP_extract(t) - 1.1) + V2_extract*ATP_extract(t)*F6P_extract(t)^2/((K2_extract*(k2_extract*ATP_extract(t)^2/AMP_extract(t)^2 + 1) + F6P_extract(t)^2)*(K2ATP_extract + ATP_extract(t))) + V1_extract*ATP_extract(t)*GLC_extract(t)/((K1ATP_extract + ATP_extract(t))*(K1GLC_extract + GLC_extract(t))) - V6_extract*ADP_extract(t)*PEP_extract(t)/((K6ADP_extract + ADP_extract(t))*(K6PEP_extract + PEP_extract(t)))
extract.AMP: diff(AMP_extract(t), t) == k9b_extract*ADP_extract(t)^2 - k9f_extract*AMP_extract(t)*ATP_extract(t) - flow_extract*AMP_extract(t)
extract.GLC: diff(GLC_extract(t), t) == -flow_extract*(GLC_extract(t) - 50.0) - V1_extract*ATP_extract(t)*GLC_extract(t)/((K1ATP_extract + ATP_extract(t))*(K1GLC_extract + GLC_extract(t)))
extract.F6P: diff(F6P_extract(t), t) == -flow_extract*F6P_extract(t) - k10_extract*F6P_extract(t) - V2_extract*ATP_extract(t)*F6P_extract(t)^2/((K2_extract*(k2_extract*ATP_extract(t)^2/AMP_extract(t)^2 + 1) + F6P_extract(t)^2)*(K2ATP_extract + ATP_extract(t))) + V1_extract*ATP_extract(t)*GLC_extract(t)/((K1ATP_extract + ATP_extract(t))*(K1GLC_extract + GLC_extract(t)))
extract.FBP: diff(FBP_extract(t), t) == k3b_extract*GAP_extract(t)^2 + V2_extract*ATP_extract(t)*F6P_extract(t)^2/((K2_extract*(k2_extract*ATP_extract(t)^2/AMP_extract(t)^2 + 1) + F6P_extract(t)^2)*(K2ATP_extract + ATP_extract(t))) - flow_extract*FBP_extract(t) - k3f_extract*FBP_extract(t)
extract.GAP: diff(GAP_extract(t), t) == -2*k3b_extract*GAP_extract(t)^2 + 2*k3f_extract*FBP_extract(t) - flow_extract*GAP_extract(t) - V4_extract*GAP_extract(t)*NAD_extract(t)/((K4GAP_extract + GAP_extract(t))*(K4NAD_extract + NAD_extract(t)))
extract.NAD: diff(NAD_extract(t), t) == k8f_extract*ACA_extract(t)*NADH_extract(t) - k8b_extract*EtOH_extract(t)*NAD_extract(t) - flow_extract*(NAD_extract(t) - 4.0) - V4_extract*GAP_extract(t)*NAD_extract(t)/((K4GAP_extract + GAP_extract(t))*(K4NAD_extract + NAD_extract(t)))
extract.NADH: diff(NADH_extract(t), t) == -k8f_extract*ACA_extract(t)*NADH_extract(t) + k8b_extract*EtOH_extract(t)*NAD_extract(t) - flow_extract*(NADH_extract(t) - 0.24) + V4_extract*GAP_extract(t)*NAD_extract(t)/((K4GAP_extract + GAP_extract(t))*(K4NAD_extract + NAD_extract(t)))
extract.DPG: diff(DPG_extract(t), t) == -k5f_extract*ADP_extract(t)*DPG_extract(t) + k5b_extract*ATP_extract(t)*PEP_extract(t) - flow_extract*DPG_extract(t) + V4_extract*GAP_extract(t)*NAD_extract(t)/((K4GAP_extract + GAP_extract(t))*(K4NAD_extract + NAD_extract(t)))
extract.PEP: diff(PEP_extract(t), t) == k5f_extract*ADP_extract(t)*DPG_extract(t) - k5b_extract*ATP_extract(t)*PEP_extract(t) - flow_extract*PEP_extract(t) - V6_extract*ADP_extract(t)*PEP_extract(t)/((K6ADP_extract + ADP_extract(t))*(K6PEP_extract + PEP_extract(t)))
extract.PYR: diff(PYR_extract(t), t) == -flow_extract*PYR_extract(t) + V6_extract*ADP_extract(t)*PEP_extract(t)/((K6ADP_extract + ADP_extract(t))*(K6PEP_extract + PEP_extract(t))) - V7_extract*PYR_extract(t)/(K7PYR_extract + PYR_extract(t))
extract.ACA: diff(ACA_extract(t), t) == -k8f_extract*ACA_extract(t)*NADH_extract(t) + k8b_extract*EtOH_extract(t)*NAD_extract(t) - flow_extract*ACA_extract(t) + V7_extract*PYR_extract(t)/(K7PYR_extract + PYR_extract(t))
extract.EtOH: diff(EtOH_extract(t), t) == k8f_extract*ACA_extract(t)*NADH_extract(t) - k8b_extract*EtOH_extract(t)*NAD_extract(t) - flow_extract*EtOH_extract(t)
extract.P: diff(P_extract(t), t) == k10_extract*F6P_extract(t) - flow_extract*P_extract(t)