COMSOL 6.4 - Oscillations in Metabolic Reaction Networks

Oscillations in Metabolic Reaction Networks

Oscillating chemical reactions were long thought to simply not exist in homogeneous solution, and even the poster child, the Belousov–Zhabotinsky reaction, met such an initial skepticism, that even though it was discovered in 1951, it took almost 20 years for it to gain widespread fame. Since this seminal discovery, many more oscillating reaction networks have been discovered, and this model studies the oscillatory behavior of a key part of glycolysis, the ubiquitous metabolic pathway common to a large part of the living organisms on earth.

Model Definition

It has been found that yeast cells, if first starved, and then subjected to small amounts of cyanide, and replenished with glucose, starts to consume the glucose at an oscillating rate. The fact that the variation in concentrations of metabolites can be measured macroscopically, means that the phase of the oscillations in individual yeast cells synchronize among billions of cells. This communication between cells occurs via molecular transport across cell membranes. The setup here models this intercellular transport via a fixed flux of glucose into the cell, acetaldehyde out of the cell, and an equilibrium of acetaldehyde across the cell membrane.

The minimal chemical reaction network presented in Ref. 1, which is formed by a subset of the glycolysis, is set up using 11 species in a Reaction Engineering interface. A 0D component, with a Reaction Engineering interface, is used to describe the reactions in the cytoplasm of the yeast cell. The individual reactions, and how they map to the reaction network is presented in Figure 1. A possible extension of this model could be to set up explicit cell membranes and solve the problem using a component in a higher dimension, resolving diffusion spatially.

Figure 1: The relation between the reaction network and the features added to the Reaction Engineering interface.

Abbreviations for the metabolites are introduced for ease of input and brevity. The mappings of short names to their respective compound(s) are listed in Table 1.

Table 1: Table of species names.
SHORT NAME	Chemical SPECIES
S1	Glucose
S2	Fructose-1,6-bisphosphate
S3	Triosephosphate, glyceraldehyde-3-phosphate, dihydroxyacetone phosphate
S4	3-phosphoglycerate
S5	Pyruvate
S6	Acetaldehyde
S6ex	Extracellular acetaldehyde
A2	ADP
A3	ATP
N1	NAD+
N2	NADH

Note how S3 represents a pool of metabolites, this is due to the model using lumped reactions to keep the number of variables low, while still reproducing the key features of the system.

Kinetics

Most reactions in the network are assigned mass action law expressions, the most prominent exception is that of the first reaction, where enzyme inhibition leads to a slightly more complicated expression, where a nonlinear factor enters the rate expression:

(1)

The final rate expressions are shown in Table 2. The influx of glucose into the cell is accounted for using an additional source.

Table 2: Rate expressions.
Rate variable	Rate expression
v1	k1[S1][A1]f([A3])
v2	k2[S2]
v3	k3fwd[A2][N1][S3] - k3rev[A3][N2][S4]
k3fwd	kGAPDHp*kPGKp/v3_denom
k3rev	kGAPDHm*kPGKm/v3_denom
v3_denom	kGAPDHm[N2]+kPGKp(A-[A3])
v4	k4[S4](A-[A3])
v5	k5[S5]
v6	k6[S6][N2]
v7	k7[A3]
v8	k8[S3][N2]
v9	k9[S6ex]

The stoichiometric coefficients associated with each rate expression are given in Figure 1.

Parameters

The parameter values associated with the kinetic model are presented in Table 3.

Table 3: Parameter values.
Parameter name	Value
J0	50.0[mM/min]
k1	550.0[1/mM/min]
Ki	1.0[mM]
k2	9.8/min
kGAPDHp	323.8[1/mM/min]
kPGKp	323.8[1/mM/min]
kPGKm	76411.1[1/mM/min]
k4	80[1/mM/min]
k5	9.7[1/min]
k6	2000[1/mM/min]
k7	28[1/min]
k8	85.7[1/mM/min]
kappa	375[1/min]
phi	0.1
A	4.0[mM]
N	1.0[mM]
n	4