4.2 Model construction
4.2.4 A kinetic model for central carbon metabolism
Our aim was to construct a kinetic model which can account for both anaerobic and aerobic metabolism, and can be used together with the datasets presented in Chapter 2. As discussed in the previous sections, the modified stoichiometric network represented in Panel A in Fig. 4.4 was constructed with the primary aim to be used within the kinetic model which is presented in Panel B.
The conceptually simplest approach to construct a kinetic model based on a stoichiometric network would be to assign an appropriate kinetic expression to each reaction in the network. However, this way is seldom chosen in practice, since each kinetic expression adds to the complexity of the model, usually mirrored in the added number of parameters.
Hence, in the kinetic model, we lumped together some reactions which are distinct in the stoichiometric network. On the other hand, we wished to represent the glycolytic pathway, the core of central carbon metabolism, in more detail than in the stoichiometric network.
The above imply that the stoichiometry of the kinetic model to be presented is not the same as the stoichiometric network we used to determine the flux distribution based on data. However the two models are constructed such that:
(i) A flux distribution in the stoichiometric network can be uniquely mapped to a flux distribution in the kinetic model.
(ii) A flux distribution which is balanced stoichiometric network is balanced in the kinetic model.
Properties (i) and (ii) imply that the flux distributions shown in Fig. 4.5 may be used in the kinetic model as steady state flux distributions. This gives five different complete sets of fluxes for the kinetic model, corresponding to the five oxygenation conditions in the chemostat cell cultures.
Flux values from these flux distributions were used in two ways in the kinetic model: Some reactions were set to a constant rate of the corresponding flux value i.e. these numbers were used as parameters of the model with knowna priori values. These reactions possess no real ’kinetic’, nevertheless may influence the behaviour of the model by consuming or producing reactants. For the rest of the model’s reactions, flux values were fitted, i.e. set as nominal quantities to be reproduced by the model. Since this involves a kinetic expression for the corresponding reaction, this approach was connected to one or more additional free parameters in the model. Panel B in Fig. 4.4 shows these two sets of reactions in different colours.
With the above motivations in mind, a kinetic model was constructed.
The reactions of the model are listed in Table 4.2, the corresponding list of kinetic expressions is presented in Eq. B.1 in Appendix B; most of these equations are taken from [87]. The balance equations, defining the species and the stoichiometry of the model, are listed in Table 4.3. A visual representation is shown in Panel B, Fig. 4.4.
The kinetic model leans on Teusink’s work [87], discussed in Section 4.1.3, which provides a list of enzyme kinetic expressions, as well as a set of corre-sponding parameters which were determined, as far as possible, on the basis of experimental data.
As in the Teusink model, the glycolytic intermediates GAP and DHAP are lumped together in a metabolite named Triose (Trio), and the reaction catalysed by triose phosphate isomerase (TPI) is assumed to be in equilibrium. Thus, the concentration ratio of this metabolites is fixed to the equilibrium constant KeqTPI, leading to the explicit formulae
DHAP = Trio/(1 +KeqTPI)
GAP=KeqTPITrio/(1 +KeqTPI) (4.13) Hence, time evolution of GAP and DHAP in this model is determined by the single variable Trio(t) and one model parameterKeqT P I. Formally, these two metabolites were not explicitly part of the model equations, Eq. 4.13 was simply used in kinetic expressions requiring GAP or DHAP concentrations.
The three adenosine species (ATP, ADP, AMP) are determined by a single dynamic variable, P(t), and the parameters ΣP andKeqAK. This is implied by two assumptions of the model: the conservation of the sum (ΣP) of these species, and the equilibrium for adenylate kinase reaction with (equillibrum constant KeqAK). These assumptions translate into algebraic constraints (see [87]) which allow to reduce the number of independent variables and calculate concentrations of all three species from a single dynamic variable of the ODE system, P(t), which can be interpreted as the concentration of high energy phosphates. The resulting formulae are
AT P(P) = (P−ADP)/2 (4.14)
ADP(P) =ΣP −p
P2(1−4KeqAK) + 2 ΣP P (4KeqAK−1) + Σ2P (1−4KeqAK)
AM P(P) = ΣP −AT P−ADP
Feasible range of the variable P is limited, as shown in Fig. 4.6.
Figure 4.6: The three adenosine cofactor species and the ATP/ADP ratio as functions
of the single variable P, according to Eq. 4.14. Parameters values are ΣP = 4.1, and KeqAK= 0.45, as used in [87]. Feasible range ofPis limited to the range left from the black slash-dotted line, in which all concentrations have positive value. The plot also highlights the role of adenylate kinase in ’signalling’ low energy states via high AMP concentrations.
points:
• A representation of the TCA cycle and respiration was added. In compari-son to the kinetic network defined by Table 4.1, metabolic processes were further simplified or omitted.
ReactionvTCA represents pyruvate transport into the mitochondria and the TCA cycle up to succinate, hence, it corresponds to the three reactions v10, v11, andv12 of the stoichiometric network.
The metabolite TCAint represents a ’generic intermediate metabolite’
within the TCA cycle. It is an inhibitor of the reaction vTCA. In this work, the nominal concentration of TCAint was set to that of citrate.
TCAint is also consumed by an anabolic reaction corresponding tov38 in the stoichiometric network, the rate of which was treated as an external reaction.
The aerobic reaction vRESP2 corresponds to Reactionv13 in the stoi-chiometric network and hence represents the reaction catalysed by the succinate dehydrogenase/complex II, as well as the corresponding proton flux through ATP synthase, resulting in ADP regeneration, and hence, creating the important stoichiometric link between these processes.
Reaction vRESPrepresents the respiratory chain. The reaction kinetics of ReactionsvRESP and vRESP2 were simplified into mass-action kinetics with a single respective parameter.
The speciesO2 was introduced as an external metabolite. Sincein vivo mitochondrial oxygen concentrations were unknown, its value was set to the value of corresponding oxygen solubility concentration given in Table 1 in [97]. It appears as a substrate in the mass-action kinetics of the Reactions vRESP andvRESP2, implying that zeroO2 concentration will halt these reactions, as it is known to be the case for the corresponding processes.
• In addition to the adenosine cofactors (ATP, ADP, AMP), the nicotinamide species NAD and NADH were introduced as dynamic variables. The sum of the these is preserved by the stoichiometry of the model, which makes the sum a conserved quantity of the model (in the sense of dynamical systems theory) and would allow to reduce the number of associated free variables to one. However, in contrast to the adenosine species, NAD and NADH concentrations were independently calculated; conservation of their sum during time evolution was observed.
• In order create a kinetic model consistent with the presented stoichiometric network, the following were added as external reactions: v29, v32, v33, v38, andv40 (associated with biomass generation), as well asv3 v5, v6, v7, v15, andv16 (associated with processes in the pentose phosphate pathway and gluconeogenesis).
• Within the branch towards ethanol, the model possesses a reaction vPDbp, representing the cytosolic pyruvate dehydrogenase bypass. (It resembles the branch towards succinate in [87], c.f. Fig. ??) The mitochondrial part
of this pathway was concluded to be inactive during the experiments, and is not represented in the model, c.f. Section 2.4, and [47].
Lumping together a number of metabolites into the species TCAint seems an acceptable approximation in this scenario, since we did not aim to explore their individual dynamics, and the change pattern of the metabolites represented in the dataset (citrate, malate, fumarate, succinate) is similar across the five conditions, even if their concentration values differ.
The metaboliteTCAinthas two roles in this model: Because of stoichiometric constrains, at least one intermediate in the TCA cycle is necessary to enable the model to have an active TCA cycle under anaerobic conditions. As reported in the flux dataset, reaction v38 consumes all of the TCA carbon flux under anaerobic condition when flux through the final part of the TCA cycle (starting with succinate dehydrogenase) is reported as zero.
Further,TCAint also plays a role in the kinetics of the model, since it acts as an inhibitor of the reactionvTCA (This inhibitory link is parametrised by vTCA_KiTCAint). This was necessary, since vTCA, being irreversible, is not effected by product concentration, hence, prevents the model from reaching a steady state by accumulating TCAintad infinitum in a number of scenarios, if not ’explicitly’ inhibited by it.
This simple regulatory loop is intended to mimic the complex regulatory mechanism of the pyruvate dehydrogenase complex which possesses a dedicated kinase-phosphatase pair, regulated by the product of the enzyme complex (c.f.
Section 3.3). While this mechanism certainly allows for much higher degree of versatility than the inhibition loop in our model, one of its functions is presumably to simply prevent harmful product accumulation e.g. when oxygen provision is cut off in the environment, causing the TCA cycle reaction succinate dehydrogenase/complex II to halt.
The complete list of kinetic equations is presented in Appendix B, Eq. B.1.
The the following list shows the three phenomenological kinetic expressions, used in reactions added to the model to represent aerobic energy metabolism.
vT CA= vTCA_Vmax
(1 + vTCA_KiTCAintT CAint )· 1
(1 +vTCA_KmPyrP yr )·(1 + vTCA_KmNADN AD ) vRESP =vRESP_K·O2·N ADH·ADP
vRESP2 =vRESP2_K·O2·T CAint·ADP·N AD vAT P ase=vATPase_K·AT P
(4.16) The parameters vTCA_Vmax, vRESP_K, and vRESP2_K were estimated, as discussed later. As discussed below, the inhibiton parametervTCA_KiTCAint could not be estimated using steady state datasets, hence its value was set to 5 mM. This value was chosen fall between concentration values for citrate (which were set as nominal values for the species TCAint) measured in the anaerobic and aerobic experiments.
name in
vPPP main PPP flux v3
v5, v6, v7 fluxes to and from PPP v5, v6, v7
vTCA pyruvate transport to mitochondria and TCA cycle up to OGA
v10, v11, v12 vRESP NADH dehydrogenase, electron transport chain,
ATP synthase
v53 vRESP2 TCA cycle from OGA, complex II, electron
trans-port chain, ATP synthase
v13
vPDC pyruvate decarboxylase v24
vADH alcohol dehydrogenase v26
vPDbp pyruvate dehydrogenase bypass v17, v18
vATPase ATP consuming processes v54
v29, v40, v32, v33, v38
flux towards anabolic processes from G6P, F6P, Triose, PEP, and TCAint, resp.
v29, v40, v32, v33, v38 Table 4.2: List of reactions in the kinetic model. In addition, the reactions catalysed by adenosine kinase (AK) and triose phosphate isomerase (TPI) are implicitly represented via parametrised algebraic constraints, as they are assumed to be in equilibrium, s.text.
PPP denotes the pentose phosphate pathway.