of 21
1
Title:
Data
-
driven model of glycolysis identifies the
role
of
allostery
in
maintaining ATP
homeostasis
Authors:
Mangyu Choe
1
,2
, Tal Einav
4
,6
, Rob Phillips
4
,5
, Denis
V.
Titov
1
,2,3
*
.
Affiliations
:
1
Department of Nutritional
Sciences and Toxicology, University of California; Berkeley, CA
5
94720, USA
.
2
Department of Molecular and Cell Biology, University of California; Berkeley, CA
94720,
USA
.
3
Center for Computational Biology, University of California; Berkeley, CA
94720, USA
.
4
Division of Biology and Biological Engineering
, California Institute of Technology;
10
Pasadena, CA 91125, USA
.
5
Department of Physics, California Institute of Technology; Pasadena, CA 91125, USA
.
6
Basic Sciences Division and Computational Bio
logy Program, Fred Hutchinson Cancer
Research Center, Seattle, WA 98109, USA
.
*Corresponding author. Email: titov@berkeley.edu
15
Abstract
:
The
specific
role
s
of allostery in regulating metabolism
are
not well understood
. Here,
we develop a data
-
driven mathematical model of mammalian glycolysis that uses enzyme rate
equations and coupled ordinary differential equations.
The
key component
s
of our model are
the
rate equation
s
for allosterically regulated
enzymes
based on
th
e
Monod
-
Wyman
-
Changeux
20
model
that we
derive
using
a
rigorous analysis of
thousands of
in vitro
kinetic measurements.
The resulting model
recapitulates the properties of glycolysis observed in live cells and shows
that the specific function of allosteric regulation is to maintain high
and stable
concentration
s
of
ATP, while glycolysis without allosteric regulation is fully capable of
produci
ng ATP and
ensuring that
ATP hydrolysis
generates energy
.
Our data
-
based modeling approach
provides a
25
roadmap for
a better
understanding of
the role of allostery in
metabolism regulation
.
One
-
Sentence Summary:
The glycolysis
model based on
allosteric
enzy
me
rate equations
recapitulates properties of glycolysis observed in live cells
.
30
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
2
Main Text:
Glycolysis is conserved across all domains of life. This key pathway harnesses the breakdown of
glucose to produce energy in the form of ATP and precursors for t
he biosynthesis of amino
acids, nucleotides, carbohydrates, and lipids.
Glycolysis contributes to ATP homeostasis by both
directly producing ATP
by a process referred to as fermentation or aerobic glycolysis
(
Fig.
1
A)
5
and by producing pyruvate
that
can be
used
as a substrate
for additional ATP
generation
by
mitochondrial respiration
.
A
erobic glycolysis
is the most active metabolic pathway in
proliferat
ing mammalian cells
(
1
)
(
an
observation known as the Warburg Effect
(
2
)
)
, where it
can satisfy
all
the ATP demand in the absence of respiration
(
3
,
4
)
.
The net reaction of aerobic
glycolysis converts extracellular glucose to two molecules of extracellular lactic acid, w
hile the
10
only
net intracellular reaction is the production of ATP from ADP and inorganic phosphate
(
Fig.
1
A)
.
The reliance of proliferating cells on aerobic glycolysis for ATP production
and the lack of
intracellular products except for ATP
makes glycolysis a convenient self
-
contained system for
s
tudying ATP homeostasis.
Glycolytic ATP production is regulated by
mass action and a constellation of allosteric effectors
15
that tune enzyme activity. Most of our knowledge about the allosteric regulation of glycolysis,
and metabolism more broadly, is based
on careful characterization of purified enzymes.
However, the
properties
of
metabolic pathways
are
the result of non
-
linear interactions between
dozens of enzymes and metabolites, making it difficult to fully understand the global effects of
regulation by
studying enzymes in isolation. The production of glycolytic ATP provides a clear
20
example of where our understanding falls short. While many allosteric regulators
of glycolytic
enzymes
have been identified, we do not know what specific
properties
of glycol
ytic ATP
production these molecules regulate.
Properties of glycolytic ATP production include matching
ATP supply and demand, maintaining ATP, ADP and phosphate concentrations such that ATP
hydrolysis generates energy, maintaining
a
high concentration of A
TP relative to ADP and AMP
25
(
5
,
6
)
,
and
enabling rapid adjustment of ATP production in response to changes in demand given
the
half
-
life of intracellular ATP
as low as
~1
-
10
sec
ond
s
(
7
)
. In addition, glycolytic ATP
production has to coordinate with
respiratory ATP production and with biosynthetic pathways.
Currently
,
we do not know which of the latter properties of glycolysis require which allosteric
regulators
,
and, in fact,
we do not know which of the
se
properties require allosteric regulation
at
30
all
.
Due to the complexity of glycolytic ATP production, a better understanding of its regulation
require
s
the use of mathematical modeling
. This has been an active ar
ea of research for several
decades, with many groups developing mathematical models of glycolysis
(
8
21
)
. These efforts
range from simple idealized models involving
a subset of
enzymes
(
8
10
,
12
,
16
)
to full
-
scale
35
models
containing all the
enzymes
(
11
,
13
15
,
17
21
)
.
Previous studies
had
two
key
limitations
that we
wanted to
address
in this report. First,
a
systematic
analysis
of the requirement of
allosteric regulators for various properties of glycolysis has
not
been performed
. Second,
enzyme
rate equations
used in the published models
, especially for the allosterically regulated enzymes,
were not
systematically
derived using
comprehensive data from multiple
in vitro
kinetic studies
.
40
Here we describe the development, analysis, and testing of a comprehensive model for the
mammalian glycolytic pathway.
Our model uses enzyme rate equations to describe glycolysis
activity using a
system of ordinary differential equations
(ODEs)
. The
key components
of our
model are the newly derived rate equations for allosterically regulated enzymes hexokinase
(HK), phosphofructokinase
(PFK), glyceraldehyde
-
3
-
phosphate dehydrogenase (GAPDH)
,
and
45
py
ruvate kinase (PK). We use the Monod
-
Wyman
-
Changeux model
(
22
26
)
in combination with
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
3
rigorous statistical approaches to identify
the relevant rate equations and to estimate
the
kinetic
constants using a manually curated dataset of thousands of
data points
from dozens of
publications.
Our model
recapitulates
all
the properties of glycoly
tic ATP homeostasis
described
above
and
accurately predicts
absolute
concentrations of glycolytic intermediates and isotope
tracing patterns in live cells.
We use
d
the
model to investigate the
role
of allosteric regulation
in
5
maintaining ATP homeostasis
. Our analysis suggests
that
allosteric reg
ulation is only required
for
the
maintenance of high ATP levels in relation to ADP and AMP, while it is dispensable
for
other properties such as matching ATP supply and demand,
ensuring
that ATP hydrolysis
generates energy, and enabling rapid adjustment of
ATP production in response to changes in
demand.
Our
analysis
identifies
the
specific mechanism that allows allosteric regulation to
10
maintain high ATP levels by controlling the rate of HK and PFK and shows
that
allosteric
regulation is absolutely required
for this property of glycolysis.
Results
M
odel
development
overview
W
e
have
developed a comprehensive
mathematical
model
of aerobic glycolysis to understand the
15
role of allostery in regulating ATP homeostasis
.
In this section, we provide an overview of th
e
model development while
most of
the
technical
details
are
described
in a comprehensive
Supplementary Text
, and Materials and Methods
.
Our model includes key allosteric regulators
based on a thorough literature search.
Most glycolytic enzymes are encoded by several
homologous genes that produce enzyme isoforms with tissue
-
specific expression and distinct
20
kinetic and allosteric properties.
We chose to focus on
glycolytic
enzyme
isoforms that are most
abundant in proliferating cell
lines
to facilitate experimental testing
(i.e., HK1, PFKP, and PKM2
isoforms of HK, PFK, and PK)
.
We sought to develop a core model of glycolysis that supports
AT
P
homeostasis
in the absence of input from other pathways, and hence we only included
allosteric regulators that are products or substrates of glycolytic enzymes.
Our model
converts a
25
qualitative
schematic of allosteric regulation (
Fig.
1
B
)
into a precise mathematical language.
As
input,
the
model
uses
i
) the
extracellular
concentration
s
of glucose and lactate,
ii
) cellular
concentrations of cofactor pools and glycolytic enzymes, and
iii
) the rate of cell
ular
ATP
consumption (
Fig.
1
C
).
We highlight these inputs
separately from the model as these inputs are
not controlled by glycolysis but depend on
the
cell type
or
experimental conditions and thus
30
cannot be predicted by
a
glycolysis model.
The model takes these inputs and uses
enzyme
rate
equations
assembled into
a system of
coupled ordinary differential equations (ODEs)
to
calculate
the
concentration of
every
glycolytic intermediate and
the rate of every reaction.
Our model can
calculate
both
steady
-
state behavior
and
dynamical
responses
to perturbations. In effe
ct, our
model
uses
in vitro
enzyme kinetics to
predict any measurable property of glycolysis and many
35
properties that cannot currently be measured
.
In total, o
ur model contains
>
1
5
0
parameters
,
including
kinetic and thermodynamic constants
,
and
estimates of enzyme and cofactor pool concentrations
, yet all of these parameters are tightly
constrained by experimental data. Kinetic constants were estimated
from
in vitro
enzyme
kinetics data
as described briefly below and in full detail in the Supplem
entary
Text
.
40
Thermodynamic constants are taken from the eQuilibrator database
(
27
)
. Enzyme and cofactor
pool concentrations are estimate
d from proteomics and metabolomics data, respectively
, as
described in Materials and Methods
. We numerically simulate the
model
using
DifferentialEquations.jl library
(
2
8
)
written in the Julia Programming Language
(
29
)
. Our code is
heavily optimized so that it takes ~1
-
10
milliseconds to calculate the results of the model under
45
given conditions
using
a single core of a
modern computer
processor
.
Optimized
code allows us
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
4
to run
the mod
el
under millions of different conditions to systematically explore the
regulation of
glycolytic ATP homeostasis
.
Fig.
1
.
Overview of
the
mammalian glycolysis model.
(A)
Coarse
-
grained description of
aerobic
glycolysis highlighting its main function to transform glucose into ATP.
(
B
)
Qualitative
5
schematic of glycolysis showing the chain of enzymes (
allosterically regulated enzymes in teal
and the other enzymes in
black) that convert substrates into products (gray).
Allosteric activators
(green)
and inhibitors
(red)
that are included in the model are highlighted
.
(
C
)
Schematic of
the
glycolysis model
including inputs and outputs.
(D)
Kinetic rate equations are shown for GPI and
PFK
and plots are actual GPI and PFK rates calculated by the respective equations
.
Note the
10
dramatic activation of the PFK rate in the presence of inorganic phosphate (Pi) and A
DP.
The defining
feature of our model is the use of newly derived kinetic rate equations to describe
the activity of
four
allosterically regulated enzymes
(i.e., HK
1
, PFK
P
, GAPDH, PK
M2
)
using
the Monod
-
Wyman
-
Changeux (MWC) model. MWC is a powerful
model
for describing the
activity of allosterically regulated enzymes
(
22
26
)
,
which assum
es that
allosteric
enzymes
exist
15
in
two or more conformations with different kinetic properties. Both the binding of substrates
and allosteric regulators can modify
the
kinetic properties
of an MWC enzyme
by stabilizing one
conformation over another. We us
e statistical approaches (i.e., regularization and cross
-
validation) to identify the simplest
MWC
kinetic rate equation that can adequately describe the
available data, allowing us to avoid overfitting and parameter identifiability issues common for
20
fittin
g complex equations to finite data
(
30
)
.
MWC equation describing
the rate of
PFK
P
in the
presence of substrates and regulators is shown in
Fig.
1
D
.
We described t
he
rest of the
glycolytic
enzymes
that are not believed to be allosterically regulated
using
standard kinetic rate equations
derived
from
quasi
-
steady
-
state or rapid equilibrium approximations.
For the eight enzymes and
transporters with one substrate and on
e product
, we used reversible
Michaelis
-
Menten
equations
25
(see equation describing the rate of glucose
-
6
-
phosphate isomerase (GPI) in
Fig.
1
D
)
and
esti
mated their kinetic constants by averaging values
from at least three publications per constant
to verify their consistency and accuracy. For the three
remaining
enzymes with more than one
substrate or product (i.e.,
aldolase (
ALDO
)
,
phosphoglycerate kinas
e (
PGK
)
,
and lactate
dehydrogenase (
LDH
)
), we
fitted
their more complex rate equations
to
manually curated
in vitro
30
kinetics
datasets
containing 350
-
700 data points per enzyme. We chose to fit data from
ALDO,
PGK,
and
LDH
instead of averaging over published kinetic constants because
many
publications
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
5
describing
ALDO, PGK,
and
LDH
activity
used different rate equations, and hence the published
kinetic constants are not directly comparable. A comprehensive description of the
derivations
and fitting of all enzyme rate equations is reported in the
Supplementary Text
, and we hope this
compendium of information will serve as a great resource for future investigations of these
enzymes.
5
Finally, we want to highlight that o
ur model u
ses several assumptions that have to be considered
when interpreting its predictions. First, the model assumes that the activity of enzymes in living
cells is accurately described using
in vitro
activity of purified enzymes. Second, the model
assumes that
enzymes and metabolites are well
-
mixed in the cytosol of the cell. Third, the model
only describes the effect of regulators that are included in the model. Fourth, the model assumes
10
that other pathways (e.g., pentose phosphate pathway, mitochondrial pyruva
te consumption) do
not affect the concentration of glycolytic intermediates or glycolytic reaction rates. Deviations
between the model and experimental observations could be due to any number of these
assumptions not being valid under the conditions of a p
articular experiment. On the other hand,
if the model accurately predicts experimental data, it suggests that these assumptions
are valid
15
under
given
experimental conditions.
Glycolysis model recapitulates key properties of ATP homeostasis
After constructi
ng the model
, we
first
investigated whether
it
can
recapitulate
the three key
properties
of ATP homeostasis
observed in living cells. First, a minimal working model must
match
ATP
supply and demand
. Second, it should maintain
a
high ATP concentration
such
that
20
most of the adenine nucleotide pool is in the form of ATP a
t
a
wide
range of
ATP turnover
s
as
observed
in vivo
(
5
,
6
)
.
Third, the model should maintain the mass
-
action ratio of ATP
hydrolysis reaction (i.e., the ratio of product concentrations [ADP]•[Phosphate] to substrate
concentrations [ATP]) far from equilibrium, which is what allows ATP to drive energy
-
demanding proce
sses in the cell.
25
We
performed simulations to determine
whether the glycolysis model can
recapitulate the key
properties of ATP homeostasis
. We report all ATPase rates in terms of percent of HK1 activity
,
which is the slowest enzyme and
glycolysis cannot p
roceed faster than the maximal rate of its
slowest enzyme (maximal rates are determined by the product of intracellular enzyme
concentration and
V
max
). First, we found that the glycolysis model produces as much ATP as is
30
consumed and rapidly adjusts to ste
pwise increase or decrease of ATP consumption, as is
necessary under physiological conditions (
Fig
.
2
A). Second, we showed that
>
90% of adenine
nucleo
tide pool was maintained as
ATP
and ATP
concentrations
vari
ed
by
<
10% when we
introduced physiologically relevant 2
-
fold stepwise changes in ATP consumption
rate
(
Fig
.
2
B).
Thus, as has been observed in muscle and heart
in vivo
(
5
,
6
)
, our model maintains nearly
35
constant ATP levels in response to large changes in ATP turnover. The latter is the key pr
operty
of feedback regulation
of
glycolysis that is required to ensure that cellular ATP
-
consuming
enzymes with
K
M
~1
-
100 μM are not affected by changes in ATP turnover. Finally, we showed
that the glycolysis model maintains
ATP, ADP and inorganic phosphat
e concentrations such that
the ATP hydrolysis reaction
is
10
9
-
10
11
-
fold away from equilibrium, equivalent to 20
-
25 k
B
T
40
(
Fig
.
2
C), which is
similar to
what has been measured in cells
(
7
)
and what allows ATP
hydrolysis to drive e
nergy
-
demanding processes. We reran the model 10,000 times using
randomly drawn model parameters and initial conditions to show that the
ability of the model to
maintain
most of the adenine nucleotide pool in the form of
ATP
levels
is largely insensitive t
o
initial conditions and to
changes
in
model parameters
within the range of experimental errors for
45
each parameter
(
Fig
.
S
1
). The coefficient of variat
ion (CV) of the
area under the curve of ATP
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
6
concentrations at a range of ATPase rates
(
1
5
%)
is lower than the average CV of a model
parameter (24%),
which suggests that maintenance of
high
ATP levels is the robust property of
the model
(
Fig
.
S
1
B).
Several kinetics models of glycolysis have been reported over the past several decades, yet none
of these reports investigated whether
those
model
s
can recapitulate
properties of
ATP
5
homeostasis except for matching ATP supply and demand
. To perform a head
-
to
-
head
comparison, we implemented a recent mechanistic model of yeast glycolysis that we will refer to
as the van Heerden model
(
19
)
. Upon repeating the simulations in
Fig
.
2
A
-
C
(
Fig
.
S
2
), the
v
an
Heerden model produced ATP in response to ATP consumption and maintained
the
high energy
of ATP hydrolysis
(
Fig
.
S
2
A, C). However, the van Heerden model could only support ATP
10
levels that are orders of magnitude lower than our model and exhibited volatility in ATP
concentration spanning almost two orders of magnitude in response
to 2
-
fold stepwise changes in
ATP demand compared to the <10% change exhibited by our model under the same conditions
(
Fig
.
S
2
B). Thus, the van Heerde
n model does not capture the feedback regulation that is needed
to maintain
most of the adenine nucleotide pool in the form of ATP
and
ensure
stable ATP levels
15
in response to changes in ATPase rate as
observed
in vivo
(
5
,
6
)
.
In summary, our glycolysis model
is the first reported model that
recapitulates all the key
properties of ATP homeostasis such as
matching
ATP
supply
and demand
,
maintaining
stable
ATP levels
such that most of the adenine nucleotide pool is in the form of ATP,
and
maintaining
high energy of ATP hydrolysis over a large range of ATP consumption rates.
20
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
7
Fig
.
2
. Glycolysis
model
recapitulates qualitative and quantitative properties of glycolysis
observed in live cells.
(
A
-
C)
Model simulations showing changes in
(A)
ATP production rate,
(B)
ATP concentration
(total adenine pool size
is labeled
with dashed grey line)
, and
(C)
energy
released during ATP hydrolysis in response to a 2
-
fold stepwise increase or decrease of ATPase
rate.
ATPase e
nergy is calculated as a natural logarithm of the
disequilibrium
ratio
(i.e.,
mass
-
5
action ratio
divided by
the
equilibrium constant
)
for the ATPase reaction
.
(
D)
Comparison of
measured glycolytic intermediate concentrations and model predictions.
Black empty circles
indicate the
experimentally determined
metabolite concentrations.
Colored l
ines with ribbon
indicate
the
median and 95% C
I of model predictions
.
T
he
shift
from left to right
for
each
metabolite level
prediction corresponds to
increasing
rates of ATP consumption
from 2% to 20%
10
of glycolysis Vmax
(see
inset for ADP
) to highlight the effect of ATP consumption on
metabolite concentrations
.
(
E)
Comparison of
[
U
-
13
C
]
Glucose and
[
U
-
13
C
]
Lactate tracing data
with model prediction
.
Circles connected with dotted lines are data and lines with ribbons are
model predictions.
ATPase rate = 15% of glycolysis Vmax was chosen fo
r model predictions.
Model predictions in (D) and (E) are
from
simulations
where bootstrapped model parameter
15
combinations could match ATP supply and demand
,
which were > 7
9
% and >9
3
% of
simulations for (D) and (E), respectively
.
Julia code to reproduce this figure is available at
https://github.com/DenisTitovLab/Glycolysis.jl
.
Glycolysis model recapitulates
live cell
data
Given that ou
r model could reproduce the wealth of measurements of each glycolytic enzyme
in
20
vitro
and
recapitulate the key properties of ATP homeostasis
, we set out to determine whether
the
model could predict the results of experiments in live cells.
To directly compare model predictions to data from bulk cell measurements, we had to adjust the
model output in two ways. First, we corrected for the fact that metabolites inside the cell can
either be bound to proteins or exist in free form in solution. C
oncentration terms in enzyme rate
25
equations refer specifically to free metabolite concentrations, so we converted free metabolite
concentrations to free + bound concentrations (using the known concentrations and binding
constants of all enzymes) before com
paring prediction to bulk cell lysis measurements. An
important caveat of this calculation is that we only consider glycolytic enzymes, which
neglects
the binding of glycolytic intermediates to other enzymes. Second, we corrected for the fact that
30
glycolyt
ic enzymes are localized to the cytosol (i.e., the part of the cytoplasm that is devoid of
organelles) and that a
significant
fraction of cytosol is occupied by macromolecules and not
water. The latter effectively means that all the enzymes and metabolites
of glycolysis are
concentrated in a volume that is ~50% of cellular volume assuming 70% of the cell is cytosol
and 70% of
the
cytosol is water.
35
We first compared the model predictions of bound and free metabolites with absolute
intracellular metabolite co
ncentrations determined by liquid chromatography
-
mass spectrometry
(LC
-
MS). Mammalian cells can synthesize ATP using
the
mitochondrial electron transport chain
(ETC) or
aerobic
glycoly
sis
. Since our current model only describes ATP
homeostasis
by aerobic
g
lycolysis
, the most appropriate comparison would be with cells that don’t have a functional
40
ETC like red blood cells (which do not have mitochondria and get their ATP exclusively from
aerobic
glycolysis) or cells treated with ETC inhibitors. In our experim
ents, we opted for the
latter, performing our experiments in the presence
of
ETC inhibitors. Yet, in addition, we also
used published data from red blood cells as well as cancer cell lines because the model
predictions suggested that the effect of glycolys
is rate on the concentrations of most
45
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
8
intermediates is smaller than the confidence intervals of model predictions
.
Indeed, we found
that
95% confidence intervals of model predictions of
most of the 18 glycolytic intermediate
concentrations at a range of AT
P
consumption
rates
overlapped
within experimental values (
Fig
.
2
D). The latter result is notable given that our model contains no direct information a
bout
intracellular metabolite levels and is free to predict concentrations from
0
to +∞
and even below
5
zero if
the
model is implemented incorrectly
.
Next, we compared [U
-
13
C]Glucose and [U
-
13
C]Lactate labeling kinetics
predicted by the model
and
measure
d
in experiments
. For
these
measurements
, we exchanged normal media for media
containing [U
-
13
C]Glucose (the input for the glycolysis pathway) or [U
-
13
C]Lactate (the final
product) and then lysed cells at different time intervals to estimate the rate and fra
ction at which
10
13
C from [U
-
13
C]Glucose or [U
-
13
C]Lactate is incorporated into glycolytic intermediates. We
observed that the glycolysis model accurately recapitulated the kinetics of
13
C labeling fraction
of intermediates after switching to [U
-
13
C]Glucose
or [U
-
13
C]Lactate
-
containing media (
Fig
.
2
E).
It is particularly noteworthy that the model can recapitulate labeling from [U
-
13
C]Lactate where
only in
tracellular lactate and pyruvate are labeled, as this indicates that our model can accurately
15
capture the extent of
reversibility
of glycolysis reactions.
Allosteric regulation of HK
1
and PFK
P
is required for the maintenance of ATP levels
Having establishe
d that the model recapitulates the key properties of ATP homeostasis and
accurately predicts measurement in live cells, w
e
next
used the model to
investigate
the
function
of alloster
ic regulation
in
maintaining
ATP
homeostasis
. We hypothesized that specific allosteric
20
regulators could serve in a number of critical roles including
matching
ATP
supply and demand
,
maintaining high and stable ATP levels
such that the majority of the adenine nucleotide pool is
in the form of ATP
, ra
pidly responding to acute perturbations of the ATPase rate, or maintaining
high energy of ATP hydrolysis.
We first removed all allosteric regulators to see if the resulting simple pathway remained
25
functional. Specifically, we made the HK
1
, PFK
P
, GAPDH, and
PK
M2
enzymes behave as
Michaelis
-
Menten
-
like enzymes with kinetic parameters corresponding to their faster (active)
MWC conformation. We also removed competitive inhibition of the HK catalytic site by G6P,
which proceeds through standard non
-
allosteric co
mpetitive inhibition. Surprisingly, the model
without any regulation was fully capable of
matching
ATP
supply and demand
and was further
30
able to maintain
the
high energy of ATP hydrolysis
(
Fig.
3
A, C). However, without allosteric
regulation, we observed a complete breakdown of ATP concentration maintenance, where ATP
levels were almost 1000
-
fold lower, and a small 2
-
fold increase and decrease of ATPase r
ate led
to an almost 10
-
fold change in ATP concentration compared to <10% change in ATP
concentration for the complete model (
Fig.
3
B).
Allosteric reg
ulation allowed glycolysis to
35
maintain >90% of the adenine pool in the form of ATP while in the absence of allosteric
regulation <1% of the adenine pool was in the form of ATP.
Systematic investigation of model
behavior at different ATPase rates showed tha
t the model without regulation is not capable of
maintaining
most of the adenine nucleotide pool in the form of
ATP levels at any ATP
ase
value
,
which is contrary to the constant ATP concentrations observed both
in vivo
(
5
,
6
)
and in the
40
complete model (
Fig.
3
D, vertical drop
-
off indicates where the ATP production by the model
cannot match ATPase rate anymore). We also note that the unregulated model can support higher
ATP production rates hinting at the potential
tradeoffs of maintaining hi
gh ATP level using
allosteric regulation
(
Fig.
3
E, solid red line is shifted to the right of the solid blue line).
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
9
We next removed the regulation of s
pecific enzymes one by one, computationally dissecting the
pathway to better understand the role of each enzyme (
Fig.
3
D
-
F,
Fig
.
S
3
).
Specifically, we
removed a particular allosteric regulator from the
relevant
kinetic rate equations by setting its
binding constant for both active and inactive MWC conformation to ∞.
We found that the
allosteric regulation of HK
1
and PFK
P
is responsible for maintaining
high and stable
ATP
5
levels, whereas removing the allosteric re
gulation of GAPDH and PK
M2
had no discernable
effect (
Fig.
3
D). Digging further, we removed each of the allosteric activators and inhibitors of
HK
1
an
d PFK
P
and found that these regulators work together to ensure the robust maintenance
of ATP levels, so that
removing any single regulator only led to a partial loss of ATP
maintenance capacity (
Fig
.
S
3
). In gene
ral,
removing inhibitors (i.e., G6P for HK
1
and ATP for
10
PFK
P
) led to worse ATP maintenance at low ATPase rates (
Fig.
3
E
,
Fig
.
S
3
E
-
H
)
, while
removing activators (i.e., phosphate for both HK1 and PFKP and ADP for PFKP) led to poorer
ATP maintenance at high ATPase rates (
Fig.
3
F
,
Fig
.
S
3
K
-
N
).
Fig.
3
.
Allosteric
regulation is required for
the
maintenance of high and stable ATP levels.
15
(
A
-
C)
Model simulations showing the effect of removal of allosteric regulation of HK1, PFKP,
GAPDH
,
and PKM2 on
(A)
ATP production rate,
(B)
ATP concentration, and
(C)
energy
release
d during ATP hydrolysis [Left axis] in response to a 2
-
fold stepwise increase or decrease
of ATPase rate [Right axis]. Energy is calculated as a natural logarithm of the
disequilibrium
ratio (i.e.,
mass
-
action ratio
divided by
the
equilibrium constant
)
for
the ATPase reaction
.
(
D
-
F)
20
Steady
-
state
ATP concentrations at a range of ATPase rate of the glycolysis model with and
without
(D)
allosteric regulation of enzymes,
(E)
allosteric activators
,
and
(F)
allosteric
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
10
inhibitor
s
.
Dashed grey line indicates the total adenine pool size (i.e., ATP + ADP + AMP).
Julia
code to reproduce this figure is available at
https://github.com/DenisTitovLab/Glycolysis.jl
.
We performed a global sensitivity analysis of our model to systematically explore the role of all
model parameters in
matching
ATP
supply and demand
, maintaining
high
ATP
levels
, and
maintaining high
energy
of
ATP hydrolysis. The goal
of
global sensitivity analysis
is to calculate
5
the contribution of each model parameter and parameter interactions to the variance of a specific
model output
(
31
)
. We used the area under the curve (AUC) for the ratio of ATP production to
ATP consum
ption, the ratio of ATP concentration to adenine pool size, and the energy of ATP
hydrolysis at a log range of ATPase values as proxies for the model’s ability to
match
ATP
supply and demand
, maintain high ATP level
s
and
maintain high
energy
of
ATP hydroly
sis,
10
respectively (
Fig
.
S
4
A, D, G). We first estimated the variance of the AUC proxies by randomly
varying each parameter of the model independently from a uniform distribution spanning a 9
-
fold range from 3 times lowe
r to 3 times higher than
the
default model parameter
values
.
The
coefficient of variation
for all three AUC proxies was in the narrow range of 0.13
-
0.4 in response
to random changes in parameter values spanning 9
-
fold, indicating that these are robust
15
prop
erties of the model
(
Fig
.
S
4
B, E, H)
.
The coefficient of variation for maintaining
high
ATP
levels was about 3
-
fold larger than for
maintaining high e
nergy of ATP hydrolysis
or for
matching ATP supply and demand
, suggesting that the former is more sensitive to changes in
specific
model parameters
.
Two metrics for each model parameter are typically reported for
variance
-
based
global sensitivity analysis
.
The first
-
order effect sensitivity index
S
1
reports the
20
fraction of variance of the model output that will be removed if the corresponding parameter is
fixed. The total
-
order sensitivity index
S
T
reports the variance that will be left if values of all but
the corresponding parameter are fixed. Larger values of
S
1
and
S
T
indicate that a given parameter
is important.
G
lobal sensitivity analysis
showed that the
K
M
and
V
max
for HK1
had
the highest
S
1
and
S
T
for all three model outputs, as expected since HK1 is the rate
-
limiting enzyme in the
25
pathway. The second largest
S
1
and
S
T
for maintaining
high
ATP levels
were given
by the
kinetic
parameters that control
allosteric regulation of HK1 and PFKP
while the s
econd largest
S
1
and
S
T
for
matching ATP
supply and demand
and
for maintaining high
energy
of ATP hydrolysis
were
simply
the
K
M
and
V
max
of PFKP
or
the second
slowest
enzyme
in the model
after HK1
based on
the product of enzyme concentration and
V
max
(
Fig
.
S
4
C, F, I).
30
In summary,
analysis of our
model
propose
s
that the function of allosteric regulation of HK
1
and
PFK
P
is to maintain high and stable A
TP levels,
while
a pathway composed of non
-
allosteric,
unregulated enzymes is sufficient to rapidly produce ATP in response to ATP consumption and
maintain high energy of ATP hydrolysis.
Mechanism of high
and stable
ATP level maintenance by HK
1
and PFK
P
al
losteric
35
regulation.
We next wanted to gain insight into the mechanism by which allosteric regulation of HK
1
and
PFK
P
allows glycolysis to maintain high and stable ATP levels.
Here, we first provide an
intuitive explanation of the mechanism that we identif
ied and then show simulations supporting
this mechanism. The concentration of ATP inside the cell is limited by adenine nucleotide pool
40
size. To maintain
a
high
ATP level, glycolysis
has
to
convert
most of the ADP
into ATP
or
,
equivalently
,
maintain
most of the adenine nucleotide pool in the form of ATP
.
The latter
automatically
leads
to
the ability of glycolysis to maintain
near
-
constant
ATP
levels
in a narrow
range of <10% (or <1%)
at
a
range of ATP turnover rates where glycolysis can main
tain >90%
(or >99%) of adenine nucleotide pool as ATP. The
relative level
of ATP and ADP concentrations
45
depends on the kinetic properties of ATP
-
consuming enzymes HK1 and PFKP and ATP
-
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
11
producing enzymes PGK and PKM2.
ATP consumption rate by the cell is of c
ourse another key
determinant of ATP and ADP levels but this reaction is not directly controlled by glycolysis.
Reactions in a metabolic pathway can proceed at
the same
arbitrary rate of 1 while being far
from equilibrium (i.e.,
forward reaction
rate of 1
and
reverse
rate of 0) or close to equilibrium
(i.e.,
forward
rate of 101 and
reverse
rate of 100).
Under the conditions when
glycolysis
proceed
s
5
in the
net
forward direction, the concentrations of products relative to substrates for each enzyme
will be hi
ghest if
the
enzyme
reaction
is close to equilibrium and lowest if
the
enzyme
reaction
is
far
from
equilibrium. Therefore, to maintain
a
high level of ATP in relation to ADP glycolysis
has to maintain
ATP
-
consuming enzymes (i.e., ATP is a substrate and ADP
is the product) HK1
and PFKP as far from equilibrium as possible and ATP
-
producing enzymes (i.e., ADP is a
10
substrate and ATP is the product) PGK and PKM2 as close to equilibrium as possible
(
Fig
.
4
A)
.
The above is the general consequence of the structure of
the
glycolysis pathway that has both
ATP
-
consuming and ATP
-
producing reaction and is independent on the specific kinetic
properties of enzymes.
Simulations show that
the
glycolysis
model behaves according to the
principle described above where
HK
1
and PFK
P
rea
ction
are kept much further from equilibrium
15
(forward to reverse rate ratio of 10
6
-
10
8
)
than all other glycolytic enzymes at a range of ATPase
rates that support high levels of ATP (
Fig
.
4
B)
.
Fig
.
4
. Allosteric regulation maintains high ATP levels by keeping ATP
-
consuming HK and
PFK reactions far from equilibrium at a range of ATP turnovers.
(A)
Schematic of the
20
glycolysis pathway showing that
in order to maintain high ATP level
HK and PFK reactions
should be
far from equilibrium and PGK and PK
close to equilibrium
.
(B)
Ratio of forward to
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
12
reverse rate for all enzymes in the full model at a range
of ATPase that support high ATP level
showing that HK and PFK are maintained far from equilibrium
. Data for each enzyme is
displayed at a 2
-
20% range of increasing ATPase rates from left to right as shown in more detail
in the inset
.
(
C)
Schematic of a simplified glycolysis
-
like pathway containing two enzymes
.
(D)
Simulation of the two enzyme pathway at 9
-
fold range of the 6 parameter controlling this
5
pathway showing that only ~0.1% of parameter combination support conditions when ATP >
A
DP
.
(E)
Ratio of forward to reverse rates for Enzyme 1, Enzyme 2 and ATPase of two enzyme
model for simulations in Panel D that support ATP > ADP
.
(F)
Ratio
s
of
!
!"#
$%&
(
,
!
!"#
$%&
)
and
!
!"#
*+,"-.
of two enzyme model for simulations
in Panel D that support ATP > ADP
.
(G)
Simulation of two enzyme model
(
"
./
$%&
(
=
"
./
$%&
)
=
100
and
"
./
*+,"-.
=
1000
)
shows that
10
maximal ATP level is achieved when rate of Enzyme 1 is exactly equal to rate of ATPase and
rate of enzyme 2 is
higher than ATPase rate.
Julia code to reproduce this figure is available at
https://github.com/DenisTitovLab/Glycolysis.jl
.
To gain a better understanding of how glycolysis maintains
high
ATP
levels
, we used a
simplified model of glycolysis, referred to as two
-
enzyme model, containing ATP
-
consuming
15
Enzyme 1, ATP
-
producing Enzyme 2 and ATPase
(
Fig
.
4
C):
Enzyme 1:
&
+
&()
,
+
&-)
Enzyme 2:
,
+
)
/01
234
+
2
&-)
6
+
2
&()
ATPase:
&()
+
7
)
8
)
/01
234
+
&-)
Note that metabolite B is phosphorylated, which conserves phosphate in each step. We describe
20
the rates of Enzymes 1 and 2 and ATPase using
a simple mass
-
action driven equation with two
parameters being
the maximal enzyme activit
y
!
!"#
$%&
(superscripts always denote the enzyme
under consideration) and the equilibrium concentration of reaction
"
./
$%&
:
9234
$%&
=
!
!"#
$%&
:
1
Γ
$%&
"
./
$%&
=
(
1
)
where
Γ
$%&
is the mass action ratio of corresponding reaction:
Γ
$%&
(
=
[
,
]
[
&-)
]
[
&
]
[
&()
]
Γ
$%&
)
=
[
6
]
[
&()
]
)
[
,
]
[
&-)
]
)
[
)
/01
234
]
Γ
*+,"-.
=
[
&-)
]
[
)
/01
234
]
[
&()
]
The two
-
enzyme model reduces the number of parameters from > 100
to 6. Only two of the
25
parameters of the two
-
enzyme model can be controlled by allostery
(
!
!"#
$%&
(
and
!
!"#
$%&
)
),
further
simplifying the analysis of the role of allostery in regulating this model
.
W
e
first
searched for values
of the six pa
rameters of the two
-
enzyme model
that would support
ATP
>
ADP.
We performed 100,000 simulations using random values of all 6 parameters drawn
from uniform distribution spanning a 9
-
fold range from 3 times lower to 3 times higher than
then
30
default model par
ameters
!
!"#
$%&
(
=
!
!"#
$%&
)
=
!
!"#
*+,"-.
=
0
.
001
,
"
./
$%&
(
=
"
./
$%&
)
=
100
, and
"
./
*+,"-.
=
1000
(
Fig
.
4
D). Only ~ 0.1% of parameter combinations
could
maintain ATP >
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint
13
ADP. By contrast, ~25% of simulations could
match
ATP
supply and demand
or
maintained
energy
of
ATP
hydrolysis above
10 k
B
T, sugg
esting that maintenance
of
ATP > ADP
is a non
-
trivial task that requires
specific combinations of the 6 parameters
(
Fig
.
S
5
A, D).
We investigated the
parameter combinations that allowed the two
-
enzyme model to maintain
ATP > ADP
.
The simulations that supported ATP > ADP
also
maintained Enzyme 1 and ATPase
5
reaction far from equilibrium and Enzyme 2 reaction close to equilibrium (
Fig
.
4
E)
as is
described for the full glycolysis model above
. We next investigated the
ratios of
parameter
s
that
allowed the two
-
enzyme model to maintain ATP > ADP, and we found
a simple relationship
!
!"#
$%&
)
>
!
!"#
$%&
(
=
!
!"#
*+,"-.
(
Fig
.
4
F).
We confirmed this relationship in
simulations where we
fixed
!
!"#
$%&
(
=
!
!"#
*+,"-.
or
!
!"#
$%&
)
=
!
!"#
*+,"-.
and varied the rate of the other enzyme (
Fig
.
4
G).
10
On the other hand, any condition where
!
!"#
$%&
(
>
!
!"#
*+,"-.
and
!
!"#
$%&
)
>
!
!"#
*+,"-.
was sufficient
to
match
ATP
supply and demand
or
to
maintain high energy of ATP hydrolysis
(
Fig
.
S
5
B, C, E,
F). Overall, our analysis suggests that to maintain high and stable ATP levels, allosteric
regulation has to keep HK
1 and
PFK
P
reactions
far from equilibrium by
controlling
the
apparent
V
max
rates of these
enzymes
such that they are
exactly equal to each other
and to
the apparent
15
V
max
of
ATPase rate
while being
lower than the rates of all other glycolytic enzymes at a range of
ATP
ase rates.
Excess activity of enzymes relative to HK and PFK is required for the maintenance of ATP
levels
We next explored what attributes of the glycolysis pathway other than allosteric regulation of
20
HK
1
and PFK
P
might be required for maintaining
most of
the adenine pool in the form of ATP
.
Analysis of the two
-
enzyme model showed that the activity of ATP
-
producing Enzyme 2 should
be as high as possible compared to the activity of ATP
-
consuming Enzyme 1 and ATPase (
Fig
.
4
G). We investigated whether a similar mechanism works in the full glycolysis model. First, we
looked at the maximal activity of glycolytic enzymes. In agreement with the two
-
enzyme model
,
25
most glycolytic enzymes have a large excess of activity compared to HK1 and PFKP
(
Fig
.
5
A)
.
In fact, most enzymes are present at 10
-
fold higher concentrations than is required for the
maximal activity of the glycolysis pathway
limited
by HK1. Such high levels of expression
repres
ent a significant investment of resources by the cell, given that expression of enzymes like
GAPDH and PKM2 can approach 1% of the proteome.
30
Using our model, we tested whether this high level of enzyme expression is required to maintain
high ATP levels.
I
n
creasing
or decreasing
the concentrations of all enzymes from TPI to LDH by
20
-
fold
more than default model values
led to no changes in ATP maintenance
(
Fig
.
5
B)
.
However, upon decreasing the expression of these enzymes
more than 20
-
fold
, the ability of the
model to maintain
more than 50%
of the adenine pool in the form of
ATP levels eventually
35
collapsed (
Fig
.
5
B)
without significantly affecting the ability of the model to match ATP supply
and demand
(
Fig
.
5
C
)
. In this way, our results suggest a mechanistic explanation for the
seemingly wasteful large excess expression of enzyme of glycolysis relative to HK and PFK.
.
CC-BY 4.0 International license
available under a
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made
The copyright holder for this preprint
this version posted December 30, 2022.
;
https://doi.org/10.1101/2022.12.28.522046
doi:
bioRxiv preprint