next up previous
Next: Metabolic rates and accounting Up: All model parameters Previous: Basics

Subsections



Estimating rate constants for reactions

We now describe some methodologies that we use frequently for estimating rate constants in reactions.


Estimating rate constants from equilibrium and turnover data for reactions not normally at equilibrium

Consider a generalised reaction of the form:
\begin{displaymath}
\sum_{i=1}^{z_1}\ensuremath{\mathrm{n_iA_i}}{{\ensuremath{\m...
...mathrm{k_{-1}}}}}} \sum_{i=1}^{z_2}\ensuremath{\mathrm{m_iB_i}}\end{displaymath} (6)

with, for maximum generality, substrate $\ensuremath{\mathrm{A_i}}$ residing in a compartment with volume $V_i$ and product $\ensuremath{\mathrm{B_i}}$ residing in a compartment with volume $\tilde V_i$. Let $a_i \equiv [A_i]$, etc. What we mean by the above representation is that we have the dynamics:
\begin{displaymath}
\ensuremath{\mathrm{\dot{a_i}}}= -\frac{\ensuremath{\mathrm{...
...) - \ensuremath{\mathrm{k_{-1}}}\mathcal{G}_2(\{b_i\})) \cdots
\end{displaymath} (7)

etc. for some functions $\mathcal{G}_1$ and $\mathcal{G}_2$. $\ensuremath{\mathrm{k_1}}$ and $\ensuremath{\mathrm{k_{-1}}}$ could represent either mass action rates, or Vmax values if the reaction actually consists of two opposing saturable processes. Quantities we might have access to are the normal turnover rate:
\begin{displaymath}
\mathcal{R} \equiv \ensuremath{\mathrm{k_1}}\mathcal{G}_1(\{...
...\}) - \ensuremath{\mathrm{k_{-1}}}\mathcal{G}_2(\{b_{i, n}\}))
\end{displaymath} (8)

and the ratio $\ensuremath{\mathrm{K_{eq}}}\equiv \ensuremath{\mathrm{k_1}}/\ensuremath{\mathrm{k_{-1}}}$ (calculable from the equilibrium concentrations of the reactants and products if it is a mass-action reaction). Then we can calculate:
\begin{displaymath}
\ensuremath{\mathrm{k_{-1}}}= \frac{\mathcal{R}}{\ensuremath...
...q}}}\mathcal{G}_1(\{a_{i, n}\}) - \mathcal{G}_2(\{b_{i, n}\})}
\end{displaymath} (9)

and similarly:
\begin{displaymath}
\ensuremath{\mathrm{k_{1}}}= \frac{\mathcal{R}}{\mathcal{G}_...
...c{1}{\ensuremath{\mathrm{K_{eq}}}}\mathcal{G}_2(\{b_{i, n}\})}
\end{displaymath} (10)

Alternatively we might have an estimate for the normal ratio of backward to forward reaction rates (i.e. of the degree of bi-directionality of the reaction at normal substrate concentrations). In such cases we often define a quantity:
\begin{displaymath}
\ensuremath{\mathrm{Rbf}}\equiv \frac{\ensuremath{\mathrm{k_...
..., n}\})}{\ensuremath{\mathrm{k_1}}\mathcal{G}_1(\{a_{i, n}\})}
\end{displaymath} (11)

In some situations - for example in the conversion of lactate to pyruvate in the brain - we know that a reaction is near equilibrium in normal conditions, but cannot actually be at equilibrium. We can use such heuristic knowledge by setting the value of $\ensuremath{\mathrm{Rbf}}$ to be near one.

As an example of the application of these methods, consider the transport of potassium ions across the blood brain barrier. This processes is known to be active, bi-directional, and possibly involving multiple mechanisms (2) . To avoid extensive complications, in the current version of the model we describe the sum of these processes very simply via two Michaelis Menten reactions:

\begin{displaymath}
\ensuremath{\mathrm{K^+_{c}}}{{\ensuremath{\mathrm{\ensurema...
...KK}}}}}}}} \ensuremath{\mathrm{\ensuremath{\mathrm{K^+_{ec}}}}}\end{displaymath} (12)

Further we assume for simplicity the same Km value in each direction. (i.e. a single kind transporter, but perhaps expressed in different quantities on each side of the barrier). We now calculate the forward and backward Vmax values for the transport of potassium from blood to brain, i.e. $\ensuremath{\mathrm{Vmax1_{KK}}}$ and $\ensuremath{\mathrm{Vmax2_{KK}}}$. Let $\ensuremath{\mathrm{T_{fKK}}}$ and $\ensuremath{\mathrm{T_{bKK}}}$ be the normal forward and backward transport rates, i.e.
\begin{displaymath}
\ensuremath{\mathrm{T_{fKK}}}= \ensuremath{\mathrm{Vmax1_{KK...
...ax2_{KK}}}\ensuremath{\mathrm{[K^+_{ec}]/(Km_{KK} + K^+_{ec})}}\end{displaymath} (13)

We must have that at steady state the difference between these rates is the same as the rate of transport of potassium out of the extracellular space during CSF drainage:

\begin{displaymath}
\ensuremath{\mathrm{T_{fKK}}}- \ensuremath{\mathrm{T_{bKK}}}= \ensuremath{\mathrm{[K^+_{ec}]}}G_o(P_{ic} - P_{vs})
\end{displaymath} (14)

where $G_o$ is a parameter representing the conductance to CSF outflow, $P_{ic}$ is intracranial pressure, and $P_{vs}$ is venous pressure. If we can estimate:
\begin{displaymath}
\ensuremath{\mathrm{Rbf_{KK}}}\equiv \frac{\ensuremath{\mathrm{T_{bKK}}}}{\ensuremath{\mathrm{T_{fKK}}}}
\end{displaymath} (15)

then we can calculate that
\begin{displaymath}
\ensuremath{\mathrm{T_{fKK}}}= \frac{\ensuremath{\mathrm{[K^...
...]}}G_o(P_{ic} - P_{vs})}{(1 - \ensuremath{\mathrm{Rbf_{KK}}})}
\end{displaymath} (16)

from which we calculate $\ensuremath{\mathrm{Vmax1_{KK}}}$ and $\ensuremath{\mathrm{Vmax2_{KK}}}$.

This methodology is used very frequently in the biochemical part of the model because as we shall see in section 3 we quite often can calculate reasonable estimates of the normal flux through a particular reaction pathway.

A comment on Vmax values in reactions

We note a fact that we use often in estimating parameters for the model: This is that reliable and reproducible estimates of maximum rates of reaction (or alternatively of enzyme concentrations and specific activities) are hard to come by for the human brain in vivo. Because we are much more likely to have reasonable estimates of the normal flux through a reaction, and because other constants such as Km values can be more reliably estimated from in vitro experimental data, we almost never use experimental Vmax values, even where they are available, and instead using the techniques above to estimate them.


Estimating rate constants for reactions from equilibrium, concentration and half-time data


One compartment

Sometimes we are ignorant about the normal flux through a reaction pathway. Alternatively, sometimes the very structure of a model implies that at model steady state, a reaction must be at equilibrium, i.e. the flux must be zero. For example, given a reaction $R$ with reactants $\ensuremath{\mathrm{A_i}}$, if for any one of the reactants, $R$ is the only reaction in which it partakes, then at steady state (if the model achieves a steady state!) $R$ must be at equilibrium.

In both of these situations, we might have knowledge of the equilibrium constant but this is not sufficient to determine the actual values of the reaction rates. For this we need further information. The quantity that we most frequently use to determine the rate constants in such situations is a putative half-time to reach equilibrium for the reaction.

To make this concept more precise, consider a generalised reaction:

\begin{displaymath}
\sum_{i=1}^{z_1}\ensuremath{\mathrm{n_iA_i}}{{\ensuremath{\m...
...mathrm{k_{-1}}}}}} \sum_{i=1}^{z_2}\ensuremath{\mathrm{m_iB_i}}\end{displaymath} (17)

For the time being we assume that all reactants and products are in one compartment, though sketch the (simple) generalisation at the end. Let $z \equiv z_1 + z_2$, $K \equiv \frac{k_1}{k_{-1}}$, and again $a_i \equiv [A_i]$, etc. Considered as a dynamical system, evolution is in the $z$-dimensional space of chemical concentrations. We are of course only interested in the area where all $a_i \geq 0$. We call this quadrant $X$, and a given point in it $x$. The equilibria of this reaction form a co-dimension 1 surface in this space $\mathcal{S} \subset X$, defined as the set of points satisfying the equation:
\begin{displaymath}
\frac{\mathcal{G}_1(\{a_i\})}{\mathcal{G}_2(\{b_i\})} = K
\end{displaymath} (18)

Any disturbance away from the surface of equilibria relaxes to a unique point on $\mathcal{S}$.

We get the dynamics:

\begin{displaymath}
\dot{a_i} = \ensuremath{\mathrm{-n_i}}(\ensuremath{\mathrm{k...
...) - \ensuremath{\mathrm{k_{-1}}}\mathcal{G}_2(\{b_j\})) \cdots
\end{displaymath} (19)

Because $\dot a_i/\ensuremath{\mathrm{n_i}}= \dot a_j/\ensuremath{\mathrm{n_j}}$ for all $i, j \leq z_1$ and similarly $\dot b_i/\ensuremath{\mathrm{m_i}}= -\dot a_j/\ensuremath{\mathrm{n_j}}$ for all $i \leq z_2, j\leq z_1$, given an initial condition we can determine the unique path by which an initial point $x_{init} = (a_{1, init}, a_{2, init}, \cdots, b_{1, init}, b_{2, init}, \cdots)$ relaxes to the final point $x_{fin} \in \mathcal{S}$. This is done simply by integrating the above relations to get:


\begin{displaymath}
\frac{\Delta a_i}{\ensuremath{\mathrm{n_i}}} = \frac{\Delta ...
...c{\Delta b_k}{\ensuremath{\mathrm{m_k}}} \quad \forall i, j, k
\end{displaymath} (20)

where $\Delta a_i \equiv a_{i, fin} - a_{i, init}$, etc. Given an initial condition $x_{init}$, the above equations define the line in $X$ along which the system evolves to reach a point on $\mathcal{S}$. The point on $\mathcal{S}$ reached is derived by considering in addition equation 18.

This geometry means that there is a unique eigenvalue associated with any given point on the surface and hence any given initial condition.

Perturbing equations 19 by letting $a_i \to a_i + \alpha_i$ and $b_i \to b_i + \beta_i$, we get the linearised equations:

$\displaystyle \dot \alpha_k$ $\textstyle =$ $\displaystyle -\left(\ensuremath{\mathrm{k_1}}\sum_{i=1}^{z_1}\left[ \ensuremat...
...{\partial \mathcal{G}_2}{\partial b_i}\right]\right)\alpha_k \quad \mbox{any k}$  
$\displaystyle \dot \beta_k$ $\textstyle =$ $\displaystyle -\left(\ensuremath{\mathrm{k_1}}\sum_{i=1}^{z_1}\left[ \ensuremat...
...c{\partial \mathcal{G}_2}{\partial b_i}\right]\right)\beta_k \quad \mbox{any k}$ (21)

That the same equation holds for all $\alpha_k$ and $\beta_k$ is not surprising considering the topology of the system which requires that any perturbation which is to approach a particular equilibrium solution must do so along a particular path.

We can work out the halftime for relaxation to equilibrium to be:

\begin{displaymath}
t_{1/2} = \frac{\ln 2}{\left(\ensuremath{\mathrm{k_1}}\sum_{...
...\mathrm{m_i}}\frac{\partial \mathcal{G}_2}{b_i}\right]\right)}
\end{displaymath} (22)

Let:
\begin{displaymath}
P \equiv \sum_{i=1}^{z_1}\left[ \ensuremath{\mathrm{n_i}}\fr...
...h{\mathrm{K_{eq}}}\equiv \ensuremath{\mathrm{\frac{k_1}{k_-1}}}\end{displaymath} (23)

Note that this definition of $\ensuremath{\mathrm{K_{eq}}}$ is not equivalent to the usual definition of the equilibrium constant unless the dynamics are mass action, but is more convenient for our purposes. So we get:
\begin{displaymath}
t_{1/2} = \frac{\ln 2}{\left(\ensuremath{\mathrm{k_1}}P + \ensuremath{\mathrm{k_{-1}}}Q\right)}
\end{displaymath} (24)

We can then solve for $\ensuremath{\mathrm{k_1}}$ and $\ensuremath{\mathrm{k_{-1}}}$ in terms of $t_{1/2}$ and $\ensuremath{\mathrm{K_{eq}}}$, i.e.
\begin{displaymath}
\ensuremath{\mathrm{k_{-1}}}= \frac{\ln 2}{(\ensuremath{\mat...
...}}}= \frac{\ln 2}{(P + Q/\ensuremath{\mathrm{K_{eq}}})t_{1/2}}
\end{displaymath} (25)

Thus if we know $\ensuremath{\mathrm{K_{eq}}}$ and have an estimate for $t_{1/2}$ at given equilibrium concentrations of reactants, then we can calculate estimates of $\ensuremath{\mathrm{k_1}}$ and $\ensuremath{\mathrm{k_{-1}}}$.

Take the particular case of a mass action reaction. In this case $\mathcal{G}_1 = \prod_{i=1}^{z_1}a_i^{\ensuremath{\mathrm{n_i}}}$ and $\mathcal{G}_2 = \prod_{i=1}^{z_2}b_i^{\ensuremath{\mathrm{m_i}}}$. So that:

\begin{displaymath}
\ensuremath{\mathrm{n_i}}\frac{\partial \mathcal{G}_1}{\part...
...od_{j=1, j\not= i}^{z_2}b_j^{\ensuremath{\mathrm{m_j}}}\right]
\end{displaymath} (26)

giving:
\begin{displaymath}
P \equiv \sum_{i=1}^{z_1}\left[\ensuremath{\mathrm{n_i}}^2 a...
...od_{j=1, j\not= i}^{z_2}b_j^{\ensuremath{\mathrm{m_j}}}\right]
\end{displaymath} (27)

As a somewhat non-trivial example we consider the generalised reaction:

\begin{displaymath}
\ensuremath{\mathrm{Hb}}+ 4\ensuremath{\mathrm{O_2}}{{\ensur...
...emath{\mathrm{k_{nXO4}}}}}} \ensuremath{\mathrm{Hb(O_2)_4}}\\
\end{displaymath} (28)

with the forward rate of reaction depending on $\ensuremath{\mathrm{[O_2]^{2.5}}}$. i.e.:
\begin{displaymath}
\dot{\ensuremath{\mathrm{Hb}}} = - \ensuremath{\mathrm{k_{XO...
...\ensuremath{\mathrm{k_{nXO4}}}\ensuremath{\mathrm{[Hb(O_2)_4]}}\end{displaymath} (29)

In this case:
\begin{displaymath}
\mathcal{G}_1 = \ensuremath{\mathrm{[Hb][O_2]^{2.5}}}\quad \mbox{and} \quad \mathcal{G}_2 = \ensuremath{\mathrm{[Hb(O_2)_4]}}\end{displaymath} (30)

and so:
\begin{displaymath}
P = \ensuremath{\mathrm{\frac{\partial ([Hb][O_2]^{2.5})}{\p...
...ial [O_2]}}}= \ensuremath{\mathrm{[O_2]^{1.5}(10[Hb] + [O_2])}}\end{displaymath} (31)

and
\begin{displaymath}
Q = \ensuremath{\mathrm{\frac{\partial [Hb(O_2)_4]}{\partial [Hb(O_2)_4]}}}= 1
\end{displaymath} (32)

It is very important to note that to do these calculations one needs to know the equilibrium constants and halftimes for the reaction occuring in isolation, rather than the steady state values when the reactants are part of a bigger scheme of things (as is usually the case with in vivo data).


Different compartments

We only sketch how the above argument generalises to the case of multiple compartments. When there are a number of compartments involved, instead of equations 19, we get:
\begin{displaymath}
\dot{a_i} = -\frac{\ensuremath{\mathrm{n_i}}}{V_i}(\ensurema...
...) - \ensuremath{\mathrm{k_{-1}}}\mathcal{G}_2(\{b_j\})) \cdots
\end{displaymath} (33)

(notation as in the previous section). We then get instead of equations 23:
\begin{displaymath}
P \equiv \sum_{i=1}^{z_1}\left[ \frac{\ensuremath{\mathrm{n_...
...\tilde V_i} \frac{\partial \mathcal{G}_2}{\partial b_i}\right]
\end{displaymath} (34)

with
\begin{displaymath}
\ensuremath{\mathrm{k_{-1}}}= \frac{\ln 2}{(\ensuremath{\mat...
...}}}= \frac{\ln 2}{(P + Q/\ensuremath{\mathrm{K_{eq}}})t_{1/2}}
\end{displaymath} (35)

as before.

Unavailability of half-time data

While equilibrium data is often available, half-time data is often not. In particular authors often simply mention in passing that a reaction is ``rapid''. If the model structure implies that at steady state the reaction must be at equilibrium (see comments above), then one possibility is simply to say that the reaction is always at equilibrium, and use the equilibrium equation to determine the concentration of one of the reactants, rather than simulating the reaction itself. An alternative possibility, which is equivalent, but computationally more robust, and works even when the reaction is not forced by the model structure to be at equilibrium during model steady state, is to set a very small half-time.


Reactions involving protons

When a reaction involves protons, the halftime as defined above - i.e. the halftime for the reaction to reach equilibrium under closed container conditions - is usually of little physiological relevance. Since pH is usually buffered by several mechanisms, a more useful quantity is probably the halftime to reach equilibrium in the situation where the pH is held constant. If for example we have the reaction
\begin{displaymath}
\ensuremath{\mathrm{X}}{{\ensuremath{\mathrm{k_1}}\atop{\lon...
...{k_{-1}}}}}} \ensuremath{\mathrm{Y}}+ \ensuremath{\mathrm{H^+}}\end{displaymath} (36)

the halftime to reach equilibrium under conditions of constant pH is:
\begin{displaymath}
t_{1/2} = \frac{\ln(2)}{\ensuremath{\mathrm{k_1 + k_{-1}[H^+]}}}
\end{displaymath} (37)

giving
\begin{displaymath}
\ensuremath{\mathrm{k_1}}= \frac{\ln(2)}{t_{1/2}(\ensuremath...
...}}}= \frac{\ln(2)}{t_{1/2}(\ensuremath{\mathrm{Keq + [H^+]}})}
\end{displaymath} (38)


Carrier mediated transport

Carrier mediated transport can be modelled using the methodology described in Section 2.4.1 of (3). If $\ensuremath{\mathrm{X}}$ exists in two compartments $A$ and $B$ with concentrations \ensuremath{\mathrm{X_a}} and \ensuremath{\mathrm{X_a}}, then the theoretical reaction scheme is:
\begin{displaymath}
\ensuremath{\mathrm{X_a}}+ \ensuremath{\mathrm{E_a}}{{\ensur...
...rrow}\atop{\ensuremath{\mathrm{k}}}}} \ensuremath{\mathrm{E_b}}\end{displaymath} (39)

where $\ensuremath{\mathrm{E}}$ is the carrier and $\ensuremath{\mathrm{P}}$ is the complex of carrier and $\ensuremath{\mathrm{X}}$. The steady state assumption leads to evolution of the form:
\begin{displaymath}
\dot [\ensuremath{\mathrm{X_a}}] = \frac{-\ensuremath{\mathr...
...rm{[X_b] + K_X(1 + p_X)}}) - \ensuremath{\mathrm{K_X^2p_X^2}}}
\end{displaymath} (40)

where $\ensuremath{\mathrm{K_X}}$ is a constant reflecting the equilibrium in the binding of $\ensuremath{\mathrm{X}}$ to the carrier, and $\ensuremath{\mathrm{p_X}}$ is the ratio of the rates of conversion between the different states of the carrier, and the rate of binding of $\ensuremath{\mathrm{X}}$ to the carrier. $\ensuremath{\mathrm{Vmax}}$ is not strictly a maximum rate - the maximum rate of transport for a given $[\ensuremath{\mathrm{X_b}}]$ is:


\begin{displaymath}
\frac{-\ensuremath{\mathrm{Vmax}}}{Vol_a(\ensuremath{\mathrm{[X_b] + K_X(1 + p_X)}})}
\end{displaymath} (41)

It is usually $\ensuremath{\mathrm{Km}}$ values which are quoted in textbooks. If we assume that $\ensuremath{\mathrm{Km_X}}$ reprents the concentration of $\ensuremath{\mathrm{X}}$ across the membrane at which we have half-maximal transport, when concentration on the other side is zero, we can calculate that:

\begin{displaymath}
\ensuremath{\mathrm{K_X}}= \ensuremath{\mathrm{Km_X}}\ensuremath{\mathrm{\frac{1 + p_X}{1 + 2p_X}}}\end{displaymath} (42)

In the case of symport, we again follow (3) and assume that all intermediate steps can be ignored. Thus if substrates X and Y are transported together we assume the reaction scheme:

\begin{displaymath}
\ensuremath{\mathrm{mX_a}}+ \ensuremath{\mathrm{nY_a}}+ \ens...
...rrow}\atop{\ensuremath{\mathrm{k}}}}} \ensuremath{\mathrm{E_b}}\end{displaymath} (43)

giving rise to dynamics of the form:
\begin{displaymath}
\dot \ensuremath{\mathrm{X_a}}= \frac{-\ensuremath{\mathrm{m...
..._{XY}(1 + p_{XY})}}) - \ensuremath{\mathrm{K_{XY}^2p_{XY}^2}}}
\end{displaymath} (44)


next up previous
Next: Metabolic rates and accounting Up: All model parameters Previous: Basics
Murad Banaji
2004-07-08