The inverse problem of brain energetics: ketone bodies as alternative substrates

Little is known about brain energy metabolism under ketosis, although there is evidence that ketone bodies have a neuroprotective role in several neurological disorders. We investigate the inverse problem of estimating reaction fluxes and transport rates in the different cellular compartments of the brain, when the data amounts to a few measured arterial venous concentration differences. By using a recently developed methodology to perform Bayesian Flux Balance Analysis and a new five compartment model of the astrocyte-glutamatergic neuron cellular complex, we are able to identify the preferred biochemical pathways during shortage of glucose and in the presence of ketone bodies in the arterial blood. The analysis is performed in a minimally biased way, therefore revealing the potential of this methodology for hypothesis testing.


Introduction
The brain, a highly metabolic organ requiring a continuous supply of nutrients to satisfy the energetic demand of its constituent cells, under normal physiologic conditions derives the energy for performing and maintaining its functions exclusively from the oxidation of glucose. During prolonged starvation, when the availability of glucose to the brain is limited, ketone bodies have been shown to be an alternate source of energy [14]. In this paper we apply a recently developed methodology to perform Bayesian Flux Balance Analysis (BFBA) [9,3,12] to a new five compartment model of the astrocyte-neuron cellular complex to determine the distribution of brain energetics during shortage of glucose and in presence of ketone bodies in the arterial blood. The main goal of flux balance analysis is the estimation of the reaction fluxes and transmembrane transport rates in and across the different compartments at steady state. Determining the parameters identifying the steady state model for brain metabolism is a very challenging underdetermined inverse problem. In vivo and in situ data is in fact difficult to gather, and in addition, little is known about the distribution of brain energetics between the two cells, as indicated by the extensive debate carried out over the last twenty years [4,11]. The unavailability of direct measurements for many of the reaction fluxes and transport rates reinforces the need for computational models which permit to test in silico different hypotheses. The importance of investigating ketone bodies metabolism, which is poorly understood, is emphasized by the evidence of their neuroprotective role in several neurological disorders and will be the focus of this paper.
Computational models proposed in the literature [7,1] lump together several key metabolic reactions, in an effort to reduce their complexity, and are typically biased by adhering to a particular hypothesis. In silico experiments and their underlying mathematical models should be maximally unbiased and the methodology should make it possible for the system itself to choose which fluxes and transports are activated and what their rates are under different physiological conditions. Our model has been constructed following this paradigm. The BFBA allows the inclusion of additional information in the form of prior constraints dictated by the physiology of the problem, without imposing them.
We explore the posterior probability density, which is the solution of the BFBA, using a hybrid of full scan Gibbs sampler [5] and Hit and Run algorithm [15]. The predictions of our model with regards to brain energetics during ketosis are presented and a few conclusions are drawn.

The five compartment computational model
Our five compartment computational model of cellular brain metabolism is based on the biochemical pathway chart shown in Figure 1 and is a modification of the one proposed in [12] and [13]. The model describes a spatially homogenized astrocyte-neuron cellular system, which differentiates between blood and extracellular space on one hand, and two cell types, astrocyte and glutamatergic neuron, on the other. Each cell type is further sub-compartmentalized into cytosol and mitochondria, and each compartment is characterized by the biochemical reactions taking place in it. Both astrocyte and neuron are equipped with a detailed glycolytic chain, including a feedback control of the phosphofructokinase enzyme activity, as described in [13]. Since we are concerned with steady state analysis, it is reasonable to assume that the blood and extracellular space are lumped together, constituting what we will refer to as "blood compartment". The compartments are linked together through transports of some metabolites between them. The glutamate-glutamine cycle, which describes the exchange and clearance of neurotransmitters in the synaptic cleft, is included in the model, and is referred to as the Vcycle. We emphasize that the shuttling of lactate from astrocyte to neuron hypothesized in [11] is not forced in either direction, but it is modelled through a mechanism of efflux/influx from the cells to the blood domain. Since in the present study we are interested in the effects of the ketone bodies (beta-hydroxybutyrate (BHB) and acetoacetate (AcAc)) on brain metabolism, we have included in each cell type the two key biochemical reactions of ketone bodies metabolism, BHB+NAD + −→ AcAc+NADH and AcAc+SCoA −→ ACoA+SUC, see Figure 1, shaded area. Ketone bodies are carried into the cells by the arterial blood, therefore two additional transports, one for AcAc and one for BHB, have been modelled in each cell type.
The dynamics of cellular brain metabolism is described by a system of ordinary differential equations based on mass balance equations, [2]. The system of ordinary differential equations describing the evolution in time of the concentrations in the blood domain is where C b is the vector containing the concentrations of the biochemical species present in the blood compartment, V b is the volume of the blood compartment, Q = Q(t) is the blood flow, and C a − C v is the vector of the differences between the arterial and venous concentrations. The vector J j c,b contains the net transport rates of the substrates between the blood and the cytosol domains of the neurons (j = 1) and astrocytes (j = 2). The dynamics in the cytosol domain of each cell type j is described by the following mass balance equation where C j c is the vector of the cytosolic concentrations, V j c is the compartment domain volume, Φ j c is the vector of reaction fluxes in the domain and S j c is the corresponding stoichiometric matrix. The vectors J j c,b and J j c,m contain the net transport rates between cytosol and blood, and cytosol and mitochondria in each cell type. The vector T j accounts for the neuronal energy demand needed to sustain the glutamate-glutamine cycle between astrocyte and neuron and only affects the cytosolic mass balance equations of ATP, ADP and Pi (T 1 ADP = T 1 P i = −T 1 ATP > 0) of the neuron. The cost of neuronal activity is coupled with the transport of the neurotransmitter glutamate, J GLU,N→A , through the equation T 1 ATP = α GLU J GLU,N→A . The value of the coefficient α GLU = .68 × 38 = 25.84 encodes the experimental information that 68% of the 38 moles of ATP produced per mole of glucose is required to transport one unit of glutamate, ([10], [6]).
Finally, the mass balance equations for the concentration vector in the jth mitochondria are Since at steady state the concentrations of the biochemical species do not change, the time derivatives must vanish. Setting the left hand-side of equations (1)-(3) to zero, we obtain a set of homogeneous linear equations with respect to the reaction fluxes, the transport rates and the arterial-venous differences. Assuming that we have access to a set of measured arterial venous difference, we consider the vector r = Q(C a −C v ) as data and we write the steady state equations as where A is a sparse matrix having the form [S M] with S the stoichiometric matrix and M the transport matrix, u = [Φ; J] is the vector containing the unknowns of primary interest and b is a vector with zeroes in correspondence to the cytosolic and mitochondrial mass balance equations, and the differences between the arterial and venous concentrations in correspondence of the blood equations. We remark that T 1 does not explicit appear among the unknowns of interest because it is expressed in terms of transport rates. The linear system (4) is underdetermined, because the nullspace of A is non trivial.
To ensure that the estimated reaction fluxes and transport rates are physiological meaningful, the unknown vector u needs to satisfy a set of inequality constraints dictated by the physiology of the problem. Hence, in addition to the linear system (4), the unknown vector u must satisfy the vectorial inequality constraints where c(u) indicates that some of the bounds might be non linear.

Description of the data and the stochastic inverse problem
The main obstacle for the estimation problem in cellular brain metabolism in humans, namely the sparsity of data due to the obvious difficulties in obtaining in vivo and in situ measurements, applies also to the study of ketone bodies metabolism. While there are measurements collected through experiments performed on rodents, it is known that ketone bodies metabolism in rats differs from that of humans, hence increasing the difficulties in using such data. In our estimation we use six arterial venous differences whose values were reported by Owen et al. [14]. This data was collected from a population of three obese patients undergoing 38-41 days of fasting.   The statistical setting is the natural environment for encoding the lack of information and uncertainties in the data, therefore we replace the deterministic linear system (4) with the stochastic additive noise model where the term e is a random noise vector that takes into account the uncertainties in the data and in the pathway network model. This information is contained in the covariance matrix Γ, where the variance δ 2 expresses our confidence that the system is at steady state and σ 2 j are the variances of the measured arterial-venous differences. If π noise is the probability distribution of the noise vector e, we can write the likelihood density of b given u as i.e., Since our only assumption about u is that it satisfies (5), the prior is constructed from the bound constraints (5), that is π prior (u) ∝ π + (Cu − c(u)), where π + is a vectorial Heaviside function taking on the value one if all the components of its argument are positive and zero otherwise. By Bayes' formula, the posterior probability density is The inverse problem is to infer on the posterior distribution, π post (u), from the likelihood π(b | u) and the prior, π prior (u), probability densities.

Exploring the posterior density
To explore the posterior probability density we generate a large sample {u 1 , u 2 , . . . , u N } of vectors distributed according to the posterior π post , using an appropriately tuned Gibbs sampling algorithm. The Markov transition rule that defines the move from the previous sample point u j−1 to the next one u j is defined as follows. We write the vector u as the sum of two mutually orthogonal vectors, u = v+w, where w is in the null space of A, i.e., Aw = 0. Since the likelihood is independent of the component w in the nullspace, we can write Similarly, the null space component is independent of the likelihood, and the new update w j is drawn from the uniform distribution over the polyhedral domain using the Hit-and-Run algorithm. The separation of the space into the null space of A and its orthogonal complement can be obtained in different ways; in this case we use the singular value decomposition of the matrix A.
Diagnostics of the convergence of the sample points is done by an output analysis of the sample histories and the normalized autocorrelation function (ACF), as illustrated in Figure 2.

Prior inequality
The inverse problem of estimating reaction fluxes and transport rates at steady state (6) is severely underdetermined and ill-conditioned, because the matrix A is rank-deficient. In the case of the steady state analysis of brain energetics in the presence of ketone bodies, the information encoded in the prior which, in the Bayesian setting, can be used to compensate for the lack of data, is in the form of loose upper bounds for the absolute values of the variables. Since we know the preferred directions of most reversible reactions and transport rates, we include this information in the prior in the form of inequality constraints. The reversible reactions lactate dehydrogenase, malic enzyme, glutamate dehydrogenase , and the bidirectional transports J b c,LAC , for which we have no a priori belief, are estimated as net fluxes which can take on positive or negative values. By this procedure we minimize the bias in the estimation process.
The feedback control mechanism exerted by some key metabolic enzymes can be incorporated in the estimation problem also in the form of inequality constraints, as proposed in [13] for the implementation of the feedback mechanism exerted by the phosphofructokinase glycolytic enzyme. We simulate the neuronal activity by setting a lower bound for the transport of glutamate. In our computed experiments we consider two activity levels by setting: J GLU,N→A ≥ 0.8 mmol/min for high activity and J GLU,N→A ≥ 0.2 mmol/min for low activity. These estimates for the lower bounds are based on the experimental results reported by Hyder et al. [10].

Results
The data consists of the arterial-venous concentration differences for glucose, oxygen, carbon dioxide, lactate, betahydroxybutyrate and acetoacetate, whose values are respectively, C a,GLC − 90mM, C a,AcAc − C v,AcAc = 0.06mM and C a,BHB − C v,BHB = 0.34mM, [14]. The standard deviations δ in the covariance matrix of the likelihood is set to δ = 0.001, and the standard deviation σ j of the observed arterial-venous differences is equal to the maximum between the 5% of the above mean values and 0.1 mM. The numerical values of reaction fluxes and transport rates are expressed in mmol/min, and are scaled to be representative of the whole human brain, whose volume is estimated at 1.5 liters.
The standard deviation for the arterial-venous differences of glucose and BHB was set equal to 0.0001, meaning that we limit the glucose supply and we force the cells to use ketone bodies. This tighter standard deviation was chosen after analyzing the output of the samples generated using a standard deviation of 0.1 for GLC and BHB. In that case the only metabolic fuel used by both cells was glucose, and no uptake of ketone bodies was observed, in obvious contradiction with what is understood about ketone bodies metabolism in the brain. In fact, when the availability of glucose to the brain is reduced and the arterial concentration of ketone bodies increases, as in our experiments, Owen et al. [14] showed that the ketone bodies could supply up to 60% of brain energy requirement, their uptake being proportional to their arterial concentration, [8].
Two samples of size 100,000, one for high activity and one for low activity, were generated saving every tenth vector from a Markov chain of size one million to obtain a sample of more independent realizations.
The preferred biochemical pathways suggested by our computational model are illustrated in Figure 3, where the posterior conditional means of the relevant reaction fluxes and transport rates are reported to summarize our results. The ticker arrows indicate the most active pathways and biochemical reactions. At high neuronal activity we predict a higher glycolytic activity in the neuron. The neuronal glucose uptake is in fact twice the astrocytic glucose influx and is equal to the intake of ketone bodies in both cells. The situation is reversed at low neuronal activity. Interestingly, the uptake of ketone bodies, which is the same at low and high activity, seems not to depend on the activity level. The rate of production of lactate by the astrocyte is almost unchanged in both scenarios, while the neuronal uptake of lactate from the blood compartment is five times faster during high activity. The oxidative phosphorylation is always more active in  the neuron, as expected because of the higher energetic demand of the glutamatergic neuron to sustain the V-cycle. Finally, in astrocyte the fluxes of malic enzyme and pyruvate carboxylation, which are always present, are more pronounced during low neuronal activity.

Conclusions
The methodology developed for performing Bayesian flux balance analysis has allowed for the first time the study of the metabolic interactions between astrocyte and glutamatergic neuron during ketosis, using a complex and detailed computational model. The preliminary results presented here are promising and might help understanding the neuroprotective potentials of ketone bodies for neurodegenerative conditions, epilepsy, hypoxia, ischemia and traumatic brain injury.