A spatially distributed computational model of brain cellular metabolism.

Journal of Theoretical Biology 376 (2015) 48–65

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

A spatially distributed computational model of brain cellular metabolism Daniela Calvetti n, Yougan Cheng, Erkki Somersalo Department of Mathematics, Applied Mathematics and Statistics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, United States of America

H I G H L I G H T S

A new spatially distributed model for brain energy metabolism is developed. A standard lumped model is recovered by spatial averaging. The sensitivity of parameter estimation to scale change is demonstrated by computed examples. The question of glucose vs. lactate as metabolites is revisited.

art ic l e i nf o

a b s t r a c t

Article history: Received 26 June 2014 Received in revised form 6 March 2015 Accepted 31 March 2015 Available online 8 April 2015

This paper develops a three-dimensional spatially distributed model of brain cellular metabolism and investigates how the locus of the synaptic activity in reference to the capillaries and diffusion affects the behavior of the model, a type of analysis which is impossible to carry out in spatially lumped models, which are shown to be consistent spatially averaged approximations of the distributed model. To avoid a geometrically detailed modeling of the complex structure of the tissue consisting of different cell types and the extracellular space, the distributed model is based on a novel multi-domain formulation of reaction–diffusion equations, accounting also for separate mitochondria. The model reduction relating the spatially distributed model and lower dimensional reduced models, including the well-mixed spatially lumped compartment model, is carefully explained. We illustrate the effects of losing the spatial resolution with a computed example which is based on a reduced one-dimensional distributed radial model, and look into how the model behaves when the locus of the synaptic activity in reference to the capillaries is changed. By averaging the fluxes and concentrations in the distributed radial model to correspond to quantities in a lumped model, and further by estimating the parameters in the lumped, we conclude that varying the locus of the synaptic activity in reference to the capillaries alters significantly the lumped model parameters. This observation seems to be consequential for parameter estimation procedures from data when the spatial resolution is insufficient to determine the locus of the activity, indicating that the model uncertainty is an inherent feature of lumped models. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Reaction–diffusion equation Multi-domain model Reduction of dimensionality Well-mixed lumped model Krogh cylinder

1. Introduction Brain energy metabolism continues to be an active research topic, providing the key to understanding the neurovascular connection, and playing a central role in the interpretation of data arising from different imaging modalities. The complexity of the brain biochemistry, and the difficulty in performing direct measurements of the

n

Corresponding author. Tel.: þ 1 216 368 2884; fax: þ1 216 368 5163. E-mail addresses: [email protected] (D. Calvetti), [email protected] (Y. Cheng), [email protected] (E. Somersalo). http://dx.doi.org/10.1016/j.jtbi.2015.03.037 0022-5193/& 2015 Elsevier Ltd. All rights reserved.

relevant quantities without disrupting its activity are some of the factors that make direct observations of the metabolic state almost impossible and therefore, to interpret the scarce indirect observations and to assimilate them into a coherent picture, mathematical metabolism models need to be employed. Mathematical models of increasing complexity and sophistication have been proposed in the past two decades to provide a narrative context in which to interpret the experimental results, and to test the feasibility of proposed mechanistic interpretation of brain metabolism. Over the years, mathematical models of brain energy metabolism have gradually evolved into complex systems with fine compartmentation and detailed metabolic pathways (Aubert and Costalat, 2005; Gruetter

D. Calvetti et al. / Journal of Theoretical Biology 376 (2015) 48–65

et al., 2001; Calvetti and Somersalo, 2011; Occhipinti et al., 2008, 2010; Calvetti and Somersalo, 2012, 2013; Cloutier et al., 2009; DiNuzzo et al., 2010; Simpson et al., 2007). The evolution of the modeling paradigm, and the rationale for the differences in the published brain metabolism models are discussed in a recent overview (Somersalo et al., 2012). The existing mathematical models of brain energy metabolism can be broadly divided into those which look at the various reaction fluxes and cross membrane exchange rates supporting a given steady state, and kinetic models which track the changes of the metabolites and intermediates over time. A feature common to all mathematical models of brain energy metabolism is the spatially lumped framework, in which a physical region of interest is represented in terms of well-mixed compartments representing, e.g., different cell types, extracellular space, and capillary blood (Aubert and Costalat, 2002; Gruetter et al., 2001; Calvetti and Somersalo, 2011; Cloutier et al., 2009; DiNuzzo et al., 2010; Simpson et al., 2007). While these models shed some light on details of the brain metabolism, they overlook some potentially important factors, including the effect of the locus of the synaptic activity in reference to capillaries, the effect of diffusion, the different roles of pre- and postsynaptic neurons, and possible variations of mitochondrial density within the cells. In principle, the different availability of glucose and oxygen in particular could affect the role of aerobic and anaerobic metabolism and use of lactate, pointing towards an important role of the diffusion. In this work, we develop a spatially distributed metabolism model that better takes into account the aforementioned factors. This paper proposes the first compartmentalized, spatially distributed mathematical model of brain energy metabolism. The compartmentation as well as the description of the reaction rates and cross-membrane exchange rates is locally in agreement with the paradigms followed by the spatially lumped models proposed in the literature, which makes it possible to reduce the distributed model into a lumped model by integrating out the spatial dimension. The new, spatially distributed model that we propose is based on a multi-domain reaction–diffusion system of partial differential equations, akin to the bi-domain models of myocardium (Henriquez, 1993; Colli Franzone et al., 2005; Roth, 1994). The mathematical framework, which allows a finer or coarser space resolution that can accommodate a more or less detailed local compartmentation at the metabolic level, is naturally suited to model regions where, for example, the density of mitochondria may change in space. Furthermore, the model allows geometric model reduction in particular geometries, making it possible to investigate the effects of spatial distribution without the need to implement the computationally challenging full three-dimensional model. To illustrate the features of the model we present results of computed examples in a cylindrical geometry modeling a Krogh cylinder around a single blood vessel (Krogh, 1919), showing that the predictions about oxygen concentrations as a function of the distance from the blood vessel are in agreement with experimental findings (Kasischke et al., 2011). Moreover, our model suggests that the different oxygen availability in the watershed of the Krogh cylinder may cause a change in glycolytic activity at different spatial locations. To understand how the lack of fine-scale spatial details may affect the interpretation of metabolic processes when a lumped model is used, we derive the latter as a spatial mean model of the fine-scale distributed model, and investigate the effect of different activation scenarios on the parameter distribution of the lumped model. One of the main findings of this analysis is that the lack of resolution in spatially lumped models seems to induce significant uncertainties to model parameter estimates. The metabolic network in the computed example is strongly simplified; however, this simplification is not a limitation of our general approach, and the example is intended only to demonstrate the effect of the spatial distribution of the model as compared to a lumped model with the

49

same biochemical complexity. In particular, the model demonstrates that ignoring the effect of diffusion may lead to parameter estimates that disagree with the underlying microscopic parameters, underlining the multi-scale nature of the metabolic models.

2. Multi-domain reaction–diffusion model Consider a bounded three dimensional domain D R3 representing a sample of the brain, with brain tissue occupying a subset Ω D, and L capillaries running through the sample. Brain tissue consists of various neurons and glial cells, surrounded by a diffusive extracellular space. To avoid having to deal with the complex details of the underlying geometry, we define a multi-domain structure as follows: let Ωk be a copy of the domain Ω with a corresponding volume fraction ηk ; 0 r k r K. We assume that

η0 þ … þ ηK ¼ 1: The domain Ωk represents a homogenized compartment that occupies a fraction ηk of the domain; a point x A Ω belongs simultaneously to all subdomains, each one of them contributing according to a weight reflecting the corresponding volume fraction. We formalize this idea by defining the kth compartment as Ωk ¼ Ω; ηk dx ; with Ωk being an identical copy of the physical domain Ω equipped with the Lebesgue measure weighted by ηk, i.e., for any integrable function f : Ω-R, we define the integral over Ωk by Z Z f ðxÞ dx ¼ ηk f ðxÞ dx: Ωk

Ω

Hence, we may interpret the volume fractions as

ηk ¼

jΩ j ; j Ωj k

with j j denoting the volume of the set. The multi-domain is then defined as

Ω ¼ Ω0 Ω1 … ΩK ; where the notation indicates that an element of Ω has Kþ1 components, one for each subdomain. The construction bears some similarity with the bi-domain theory in a different context, see Henriquez (1993), Colli Franzone et al. (2005), and Roth (1994). The model problem with K¼2 consists of three subdomains, the extracellular space (ECS), (k¼ 0), neuron (k¼ 1), and astrocyte (k¼2). However, the formalism extends to more detailed models by further subdividing these subdomains to comprise, e.g., pre- and postsynaptic neurons, and glutamatergic and GABAergic neurons. In general, the tissue domain Ω and the blood in the capillaries exchange certain metabolites, which diffuse within each subdomain. In turn, the subdomains within the tissue may exchange some of the metabolites: these processes are described in terms of a reaction–diffusion model. Consider a single metabolite A, for example oxygen, glucose, or lactate. Suppressing the reference to the particular metabolite, we denote by C k ðt; xÞ the concentration in the subdomain Ωk at time t at x A Ω. The metabolite may diffuse within the subdomain Ωk with a diffusion constant Dk which is characteristic to the subdomain and to the metabolite. In the general model the diffusion coefficient Dk can be space and time dependent, and anisotropic. If there were no exchanges between subdomains and no biochemical reactions involving the metabolite took place, the process which changes the concentration of Ck over time could be described by a diffusion

50


produces si 4 0 molecules of A, or consumes si 4 0 molecules of A with a reaction rate ψ i ¼ ψ i ðt; xÞ, we add to Eq. (5) a reaction term

equation:

ηk

∂C k ¼ ∇ Dk ∇C k : ∂t

ð1Þ

To assign the boundary conditions for (1), we need to decide how to model (a) the interface between the blood vessels, and (b) the fictitious boundary enclosing the computational domain, denoted by S0. Starting with the latter question, assuming that the domain of interest Ω is independently nourished through the capillaries traversing it, it is appropriate to impose a homogeneous Neumann condition (no diffusion in or out of Ω): k ! k k ∂C ð2Þ D ¼ n ðDk ∇C Þ 0 ¼ 0; ∂n 0 S S

! where n is the exterior unit normal to S0. The boundaries between Ω and the L capillaries correspond to the blood–brain– barrier (BBB), and hence must be carefully modeled. Our preliminary model assumes that only the ECS exchanges substrates with the capillary blood. This simplifying assumption is not necessary, and its validity is briefly discussed in Section 5. Denoting by Sℓ the boundary between the ℓth capillary and Ω, we set a no-flux boundary condition for concentrations in the tissue compartments: k ! k k ∂C D ð3Þ ¼ n ðDk ∇C Þ ℓ ¼ 0; 1 r k r K; 1 r ℓ r L; ∂n ℓ S S

while allowing exchanges with the ECS: ∂C 0 ℓ D0 ¼Φ ; ∂n ℓ

ð4Þ

1 X 1X ∂C k 1 ¼ k ∇ Dk ∇C k k φk-j φj-k þ k si ψ i ∂t η η j η i

scaled by the volume fraction to make the same biochemical reaction in different subdomains comparable and independent of the volume fraction. 2.1. Adding mitochondria The assumption that the subdomains Ωk may communicate with each other and that diffusion may take place within the subdomains are reasonable when the domain Ω represents a volume comprising few cells, with one of the subdomains representing the ECS and the others the cytoplasms of neurons and astrocytes. At this scale, mitochondria can be thought of as inclusions in the cytoplasm that do not communicate with each other, but exchange substrates only with the cytoplasm surrounding them. This is what motivates the following modification of the model when the mitochondrial subcompartments are added. The ECS does not contain mitochondria, while the other subdomains do. We denote the mitochondrial domains by

ωk Ωk ;

ℓ

1 X ∂C k 1 ¼ k ∇ Dk ∇C k k φk-j φj-k ; ∂t η η j

γk ¼

φ

¼φ

k-j

j ωk j jΩ j k

:

We remark that when separate mitochondria are added to the subdomain Ωk, the volume occupied by cytoplasm must be reduced accordingly, thus the volume fraction of the latter becomes

ηk -ð1 γ k Þηk : This scaling must be implemented in the kinetic models describing diffusion, reactions and exchange. The mitochondrion density in a cell can be approximated in terms of the factor ηk. If vmito represents the volume of a typical mitochondrion, the number of mitochondria in the subdomain Ωk is # of mitochondria in Ω ¼ k

ð5Þ

where the sum extends over all compartments communicating with Ωk. We assume that the exchange terms are local both in time and space, that is k-j

k ¼ 1; …; K:

The volume fractions of mitochondria are typically rather small, comprising only a small percentage of the cell, which we denote by

S

where the flux density Φ is related to the cerebral metabolic rate ! of the species in question. The vector n normal to the boundary ℓ surface S points from tissue into the capillary, thus a positive flux ℓ density Φ corresponds to a flux from blood to tissue. The dependency of the flux density on the concentrations in blood and ECS will be discussed later. Consider the fluxes between the subdomains. For each metabolite which can be exchanged between compartments k and ja k, the corresponding diffusion equation (1) must be modified by adding at each point x A Ω a pair of exchange terms, yielding

ðt; xÞ:

The scaling by volume fractions is necessary for mass balance: if no diffusion would take place, the rate of change of mass in the kth compartment occupying a volume B Ω due to solely the transport from compartment k to compartment j would be Z Z Z d k ∂C k 1 mB ¼ ðt; xÞ dx ¼ k φk-j ðt; xÞηk dx ¼ φk-j ðt; xÞ dx; k dt η B B \ Ω ∂t B and an equal expression but with opposite sign would appear in compartment j: Z d j mB ¼ φk-j ðt; xÞ dx: dt B In the next step we adjust the model to account for the effect of biochemical reactions on concentrations. Assume that at x A Ω in compartment k; 1 r k r K, the metabolite A participates in reactions which either depletes it, or produces it. To indicate that the ith reaction

ð6Þ

j ωk j ; vmito

which, divided by the volume of the subdomain, yields the mitochondrion density

ρk ¼

j ωk j vmito j Ω j k

¼

γk vmito

:

In general, if the mitochondrion density is not assumed to be constant over the cytoplasm, it is denoted by γ k ¼ γ k ðxÞ. In this paper, we will assume that γk is constant over Ω. Assume that substrate A is exchanged between cytoplasm and mitochondria of subdomain Ωk. The kinetic model (6) of A is augmented accordingly: X ∂C k 1 1 ¼ ∇ Dk ∇C k φk-j φj-k ∂t ð1 γ k Þηk ð1 γ k Þηk j X 1 1 ϑk þ ϑk þ si ψ i ; k k k k ð1 γ Þη ð1 γ Þη i ∂ck 1 kþ 1 X r r k ¼ k k ϑ ϑ þ k k sψ ; ∂t γ η γ η r

ð7Þ

ð8Þ

where Ck denote the cytosolic concentration and ck denote the k; 7 mitochondrial concentration of A; ϑ is the exchange term into


(þ ) and out of ( ) the mitochondria, and the summing index r runs over the indices referring to mitochondrial reactions. 2.2. Capillary dynamics Consider the ℓth capillary, denoted by Rℓ D, 1 r ℓ rL, whose boundary consists of the capillary wall Sℓ and the capillary ends, Sℓa at the arterial end and Sℓv at the venous end. The concentration of a given species in Rℓ is denoted by C ℓb , and is assumed to satisfy the advection–diffusion equation: ∂C ℓb ! ¼ Db ΔC ℓb v ∇C ℓb ; ∂t

! where the velocity of the blood through the capillary v satisfies the boundary condition: ! ! ð9Þ n v ℓ ¼ 0:

51

Table 1 Reactions in neuron and astrocyte. The substrate of each reaction controlling the reaction rate in Eq. (15) is listed first on the right. Observe that in the last reaction of ATP dephosphorylation, ATP is considered as substrate. Reaction

Compartment

Glc þ 2 ADP þ 2 NAD þ þ 2 Pi-2 Pyr þ 2 ATPþ 2 NADH Pyr þ NADH-Lac þ NAD þ Lac þ NAD þ -Pyr þ NADH Pyr þ 5 NAD þ þ ADP þ Pi-3 CO2 þ 5 NADH þ ATP O2 þ 2 NADH þ 5 Pi þ 5 ADP-2 NAD þ þ 5 ATP þ 2 H2 O Cr þ ATP-PCr þ ADP PCr þ ADP-Crþ ATP Gln-Glu Glu þ ATP-Gln þ ADP þ Pi ATP-ADP þ Pi

Neuron and astrocyte Neuron and astrocyte Neuron and astrocyte Neuron and astrocyte Neuron and astrocyte Neuron and astrocyte Neuron and astrocyte Neuron Astrocyte Neuron and astrocyte

S

This condition is equivalent to stating that the component of the velocity of the blood in the direction orthogonal to the capillary wall vanishes, hence the blood velocity is tangential to the capillary wall. The boundary conditions for the concentration at the end of the capillaries are C ℓb Sℓ ¼ C ℓa ; C ℓb Sℓ ¼ C ℓv ; a

v

where C ℓa and C ℓv are the arterial and venous concentrations, respectively. Because diffusion across the capillary wall into ECS needs to match the uptake flux (4): Db

∂C ℓb ∂n

ℓ

¼Φ ;

where the normal vector is exterior to capillary.

ð10Þ

Ω0 and points into the

2.3. Exchange rates, reaction rates Consider the exchange terms φk-j and φj-k in more detail. If the transfer of the species A occurs by passive diffusion, we may use a basic Fick's law: φk-j φj-k ¼ λ C k C j ; ð11Þ where λ is specific for the species, and could possibly depend on x. For the moment, we assume a homogeneous model in which λ is constant. Observe that the reciprocal of λ is a time-like quantity, the characteristic time of the transmembrane diffusion between the compartments. Most of the species of interest do not diffuse across the membranes, but require the presence of particular membranebound transporters. Active, or facilitated, exchange is modeled by the Michaelis–Menten type formula: ! Ck Cj φk-j φj-k ¼ T max ; ð12Þ M þ C k M þC j where T max is the maximum exchange velocity that depends, e.g., on the density of the transporters, and M is the affinity constant. The formula above assumes symmetry of the exchange process. The characteristic time constant for the facilitated transfer is M=T max . All the kinetic parameters are species and subdomain specific. The fluxes between the blood and ECS compartment are treated similarly. If the flux is diffusive, as is commonly assumed for gases (Boron and Boulpaep, 2012), we write ð13Þ Φℓ ¼ ϕRℓ -ECS ϕECS-Rℓ ¼ λBBB C ℓb C ECS ; where “BBB” stands for Blood–Brain–Barrier. Similarly, for species such as glucose or lactate requiring transporter mechanisms that

are subject to saturation, we choose a Michaelis–Menten type condition: ! C ℓb C ECS Rℓ -ECS ECS-Rℓ ℓ BBB Φ ¼ϕ ϕ ¼ T max : ð14Þ M BBB þ C ℓb M BBB þ C ECS Consider a simple enzymatic reaction R in which a substrate A with concentration [A] is transformed into species B. In our model, the flux ψ of such reaction is expressed in terms of a simple Michaelis–Menten form: R : A-B;

ψ ¼ V max

½A ; K þ ½A

ð15Þ

where V max is the reaction specific maximum flux with unit ½V max ¼ mM=min and K is the affinity coefficient with unit mM. The reactions included in the current model are listed in Table 1. A modification is introduced for reactions that change the phosphorylation state r phos or the redox state r redox of the system, defined as r phos ¼

½ATP ; ½ADP

r redox ¼

½NADH : ½NAD þ

In the reaction flux ψ for every reductive reaction (involving the cofactor reaction NADH-NAD þ ), we modify formula (15) by multiplicative factor: r

ψ - redox ψ ; ν þ rredox while for oxidative reactions (involving the converse cofactor reaction NAD þ -NADH), we modify the flux as r1

ψ - redox ψ: 1 ν þ rredox The dimensionless affinity parameters ν are reaction specific. Similarly, in reactions involving dephosphorylation (ATP-ADP), or the converse, phosphorylation (ADP-ATP), as a cofactor reaction, the corresponding multiplicative factors are r phos μ þ rphos

or

1 r phos

1 μ þ rphos

;

respectively, with reaction specific dimensionless affinity μ. Observe that some reactions may involve more than one cofactor reaction, requiring therefore more than one cofactor coefficient in the corresponding reaction flux formula. We point out that although not indicated in the formulas above, all kinetic parameters are species and subdomain specific. Although they can depend spatial location, e.g., through different densities of transporters in different areas of the tissue, we treat them as constants in space.

52


2.4. Dimensions and volume scaling The volume scaling in multiscale models sometimes becomes a source of confusion, so it is worthwhile to clarify our assumptions. In brain metabolism research, it is common to scale the quantities so that they correspond to one gram brain tissue. The density of brain tissue is approximately ρ ¼ 1:04 g=cm3 , which means that one gram of tissue occupies a volume of approximately V ¼ 0:96 cm3 ¼ 0:96 mL. In our discussion, we assume that the exchange and reaction fluxes are all scaled to correspond to this volume. The times, following the conventions on the field, are usually expressed in minutes (min, and concentrations are expressed in millimolars (mM), that is, 1 mM ¼ 10 3 mol=L ¼ 1 μmol=mL. The Michaelis–Menten affinity constants, regardless of whether we refer to a distributed or lumped model, are given in volume-independent units: ½K ¼ ½M ¼ mM;

½T max ¼ ½V max ¼ mM=min;

and similarly, the unit of the membrane permeability that

λ is scaled so

½λ ¼ 1=min: Diffusion coefficients are given in units ½D ¼ cm2 =min:

3. Lumping as a way of reducing spatially distributed models The importance of the spatial structure withstanding, there are two main reasons why spatially lumped models continue to be favored over more realistic, spatially distributed ones: one is the limited availability of parameter values, which makes the identification and validation of complex models particularly challenging, and the other one is the great increase in computing resources needed to run more realistic simulations of the brain. A question which arises naturally in this context is how changes in the locality of neural activity should be taken into account in spatially lumped mathematical models of brain energy metabolism which can be found in the literature: see, e.g., Aubert and Costalat (2005), Cloutier et al. (2009), DiNuzzo et al. (2010), and Simpson et al. (2007). To understand if, and in which way, a spatially lumped model is sensitive to changes in the unresolved spatial dimension, we now derive the standard lumped model by performing model reduction on a spatially distributed one through averaging. A schematics of this process is indicated by the diagonal arrow in Fig. 1. To further understand how the collapse of the spatial resolution may affect the spatially lumped model, we construct an intermediate, one dimensional model which allows us to control the location of the neural activation while keeping the computational complexity within reasonable limits. In the reduction from a threedimensional to a one-dimensional model, indicated by the top arrows in Fig. 1, the number of unknown parameters is much smaller, reducing significantly the run time of simulations. While many of the essential properties of the model are retained in the dimension reduction, some will be inevitably lost. A more realistic intermediate reduction from three to two spatial dimensions, bottom arrows in Fig. 1, which relies on finite elements to solve the associated governing equations, will be studied in a later work. 3.1. Spatially lumped model as a mean approximation of a three dimensional model In this section, we show that the lumped model can be seen as an average over the multi-domain model developed above. The derivation provides means to interpret multi-domain simulations in the context of the lumped model, and it also shows that the new

Fig. 1. The model reduction chart. The spatially distributed model and the spatially lumped model are connected through proper reduction (diagonal path). Two ways of reduction can be performed to lower the dimension of the problems: one is indicated in the top path, which reduces the 3-D system into 1-D system. Then the 1-D system can be solved using ODE solvers; the other is indicated in the bottom path, which reduces the 3-D system into 2-D system. Then the 2-D system can be solved using finite element method.

model is not in conflict with well-established lumped models developed previously in the literature. We consider, for simplicity, the reaction–diffusion model in Ωk for a sub-compartmentation which does not account for separate mitochondria. Adding mitochondria is a straightforward extension. k Define the mean concentration C of a given species in k subdomain Ω to be Z Z Z 1 ηk 1 k k k C ðtÞ ¼ C ðt; xÞ dx ¼ C ðt; xÞ dx ¼ C k ðt; xÞ dx: k k j Ωj Ω j Ω j Ωk jΩ j Ω Integrating both sides of the differential equation (6), we obtain, by the divergence theorem, Z k dC 1 ∂C k ¼ ðt; xÞ dx dt j Ω j Ω ∂t Z Z X 1 1 Dk ΔCðt; xÞ dx þ k φj-k ðt; xÞ φk-j ðt; xÞ dx ¼ k η j Ωj Ω η j Ωj Ω j

Z X 1 si ψ i ðt; xÞ dx ηk j Ω j Ω i ! Z XZ 1 ∂C 1 X j-k þ φ ðtÞ φ k-j ðtÞ Dk dS þ k ¼ k ∂n η j Ω j S0 ℓ Sℓ η j þ

þ

1X i i s ψ ðtÞ

ηk

ð16Þ

i

where the average exchange rates and reaction fluxes are defined as Z Z 1 1 φ k-j ðtÞ ¼ φk-j ðt; xÞ dx; ψ i ðtÞ ¼ ψ i ðt; xÞ dx: ð17Þ j Ωj Ω j Ωj Ω We define the total flux from ℓth capillary, and the total flux from the blood domain to ECS by the formulas Z Z ℓ ∂C Φ ¼ ℓ D0 dS ¼ ℓ Φℓ dS; ∂n S S X ℓ Φ¼ Φ ; ð18Þ ℓ

By using the no-flux boundary condition at S0, we then obtain the familiar lumped model equation: 1X dC 1 1 X j-k ¼ δ0;k Φþ k φ ðtÞ φ k-j ðtÞ þ k si ψ i ðtÞ: 0 dt η η i jΩ j j k


53

Consider next a single capillary Rℓ . We define Z 1 ℓ C b ðtÞ ¼ ℓ C ℓb ðt; xÞ dx: j R j Rℓ

microscopic and macroscopic parameters is no longer straightforward. More precisely, consider the exchange model (11). The lumped equivalent of it is

The differential equation for this average concentration is Z ℓ dC b 1 ! ¼ ℓ ∇ Db ∇C ℓb v C ℓb dx dt j R j Rℓ ! Z Z ! Z 1 ∂C ℓb ! ! ℓ n v C a dS; þ þ ¼ ℓ Db ∂n jR j Sℓ Sℓa Sℓv

φ k-j φ j-k ¼ λ C C ;

To simplify the integrals over the capillary ends, we assume that the blood transport is advection-dominated, so we can neglect the diffusion part and write Z !

Z

j Rℓ j

Sℓa

þ

Sℓv

Db

j

which is readily obtained by averaging (11) over Ω. In contrast, for the facilitated transfer such as (12), the averaged version no longer follows by linearity. We may write a saturation model

where the minus-sign in front of the integral over Sℓ is due to the convention that the normal vector along the capillary wall points into the capillary. By imposing the boundary conditions (9) and (10), we obtain ! Z ! Z ℓ ℓ dC b 1 1 ∂C ℓb ! ! ℓ ¼ ℓ Φ þ ℓ n v C a dS: þ Db dt ∂n jR j jR j Sℓa Sℓv

1

k

! ℓ ∂C ℓb ! ! ℓ ! A n v C a dS j v j ℓ C ℓa C ℓv ¼ qℓ C ℓa C ℓv ; ∂n jR j

φ k-j φ j-k ¼ T max

C

k

M þC

k

C

!

k

M þC

k

;

but the macroscopic parameters T max and M need to be no longer equal to the microscopic ones, T max and M. This observation, although quite obvious, has a significant consequence: Michaelis–Menten parameters estimated from cell level laboratory experiments need not be good values for lumped models, making the use of literature values significantly challenging. This point will be elucidated in computed examples.

3.2. Reduction of dimensionality: a radial model

ℓ

where A is the area of the lumen of the capillary. Here we assume that the velocity field over the capillary is approximately constant and normal to the capillary ends. The quantity qℓ is given in units 1/min. The mean concentration in the capillary therefore satisfies ℓ dC b 1 ¼ ℓ Φ þ qℓ C ℓa C ℓv : dt jR j ℓ

Finally, as is customary in the lumped models, we combine the capillaries into a single blood domain. Assuming that each capillary is identical, that is, qℓ ¼ q, and j Rℓ j ¼ V b =L, where V b is the total volume of the blood compartment, the average concentration over capillaries satisfies dC b 1 Q ¼ Φ þ ðC a C v Þ; Vb Vb dt where Cb ¼

L 1X ℓ C ; Lℓ¼1 b

Φ¼

L X

ℓ

Φ ;

Q ¼ qV b :

ℓ¼1

The quantity Q is the total blood flow through the tissue volume, with the dimension ½Q ¼ mL=min. Observe that the total flux Φ out of the blood compartment matches the total flux into the ECS domain given by (18). It is customary in lumped models to specify only the arterial or the venous concentration, otherwise the species uptake will be determined by the arterial–venous difference rather than by the energy demand of the system. Therefore, the concentration in blood is expressed as a convex combination of the arterial and venous values: C b ¼ ð1 FÞC a þ FC v ;

0 r F r 1;

with the coefficient F being the mixing ratio. With this notation, we may write C a C b ¼ FðC a C v Þ;

C b C v ¼ ð1 FÞðC a C v Þ;

leading to the alternative formulas: dC b 1 Q 1 Q ¼ Φþ Ca C b ¼ Φ þ C C v ; Vb FV b Vb ð1 FÞV b b dt

ð19Þ

Kinetic lumped models are commonly equipped with reaction and transport laws similar to the ones we use in the distributed model (Section 2.3), while the connection between the

We consider a computed example in a simplified geometry: the tissue occupies a cylindrical domain around a single capillary. The capillary of length h is assumed to run parallel to the z-axis and to have a constant radius δ 4 0. The outer radius of the tissue cylinder is R 4 δ, and it is assumed that the tissue in the cylinder draws all its nutrients from the central capillary. The model corresponds physiologically to a Krogh (1919) cylinder. Two-photon nicotinamide adenine dinucleotide (NADH) imaging enables direct visualization of the Krogh tissue cylinder and the typical radius of a Krogh cylinder is R 50 μm (Kasischke et al., 2011). This value is used in our computations. The capillary length of the cylinder is denoted by H, the value of which turns out to be of no importance here. It is pointed out in Kreft et al. (2012) that “the astrocyte may be regarded as a crowded place; metabolites and organelles are moving around in a morphologically and functionally complex manner and diverse cell in a semi-aqueous cytosolic environment with internal physical barriers (i.e. organellar membranes) to isotropic diffusion”, which makes the diffusion process much slower compared to the diffusion in the ECS compartment. Furthermore, we point out that the cellular compartments do not represent single cells, but rather the agglomerate of possibly numerous cells separated by cell membranes in the unit, complicating the diffusion modeling. On the other hand, extracellular space is commonly agreed to play the role of diffusion channel in the brain (Nicholson and Syková, 2008). Motivated by these observations, our simplified model neglects the diffusion process in the neuron and astrocyte compartment. This simplification of the model can be easily removed, and it is further commented in Section 5. The general three-dimensional model can be reduced into a one-dimensional radial model by averaging over the vertical and angular variables in cylindrical coordinates, i.e., by defining C ðr; tÞ ¼

1 2π H

Z 0

H

Z 2π 0

Cðr; θ; z; tÞ dθ dz;

ð20Þ

and derive the one-dimensional reduced model similar to the lumped model. However, rather than going through the radial diffusion model, we derive the discrete computational model directly for the averaged concentrations (20), leading to a timedependent system similar to that produced by the method of lines. For notational simplicity, we drop out the overline of the averaged concentrations (20).

54


3.2.1. Subunit geometry We subdivide a Krogh cylinder into N subunits along r, each of which occupies the same volume. These subunits interact with each other through diffusion and each subunit is compartmentalized into neuron, astrocyte, synaptic cleft and extracellular space. The blood compartment is taken into account only in the innermost subunit adjacent to the capillary. Furthermore, we assume that the concentrations over a compartment in the specified subunit are constant. At each discretized unit, neuron, synaptic cleft and astrocyte communicate with its associated extracellular space. With no diffusion in neuron and astrocyte, the coupling between units is mediated through diffusion in ECS. A lumped model for the blood compartment will be used, thus neglecting the capillary geometry, and we assume a line-like capillary with radius δ ¼ 0. We shall use the lumped model (19) for the blood concentrations, however, regardless of the geometric simplification, the blood volume V b , is non-zero. The total flux, which is volume specific, is denoted by Φ, and we have Vb

dC b Q Ca C b : ¼ Φþ F dt

ð21Þ

Observe that the flux Φ goes from blood to ECS only. After subdividing the cylinder into N annular units, we make sure that all units have the same volume by setting

π ðr2j r 2j 1 ÞH ¼ π ðr2j þ 1 r 2j ÞH;

r 0 ¼ 0; r N ¼ R;

where rj is the outer radius of the jth subunit. It follows that, for 1 r j r N, rffiffiffiffi j R; rj ¼ N and if we let

be the thickness of the subunits, it follows that rffiffiffiffi rffiffiffiffiffiffiffiffiffi j j1 : hj ¼ αj R; αj ¼ N N

V π R2 H : Vu ¼ ¼ N N

DSj ECS 4Dπ r j H ECS ðC C ECS Þ¼ ðC C ECS Þ; j j hj þ hj þ 1 j þ 1 ℓj j þ 1

1 r jr N 1:

Here Sj and ℓj are the cross sectional area and the average distance between the jth subunit and the (jþ 1)th subunit, respectively, given by hj þ hj þ 1 : 2

In the first and the last subunit, we set

Φ0 ¼ Φ;

Φ0

S2

2 Φ0 3 0 Φ4

S3 S

4

Fig. 2. An axial view of the flux representations in the discretized Krogh cylinder geometry when N ¼ 5. In this case, the thickness of each subunit is h1 0:45R; h2 0:19R; h3 0:14R; h4 0:12R; h5 0:11R. Here Φj ð1r jr 4Þ is the diffusion flux in ECS between the jth subunit and the ðj þ 1Þth subunit and Φ0 is the exchange rate through the blood brain barrier.

nutrients from blood, and that no diffusion takes place at the outer boundary of the Krogh cylinder. After accounting for the exchange rates from ECS to neuron (n) and astrocyte (a) and vice versa, from neuron and astrocyte to ECS, the total P volume-specific flux in and out of ECS in subunit j are x ¼ n;a J x-ECS and j P ECS-x respectively, and we have the mass balance equation x ¼ n;a J j ECS

ηECS V u

dC j

¼ Φj 1 Φj

dt

X ECS-x x-ECS : Jj Jj

x ¼ n;a

after substituting the expressions for jr N 1,

ΦN ¼ 0;

where Φ is the total flux out of the capillary, in agreement with the understanding that the subunit adjacent to the capillary receives the

Φ we have that, for 2 r

ECS

ηECS V u

3.2.2. Extracellular space In the simplified radial model, the extracellular space is the only diffusive compartment. To derive the mass balance equations for the jth subunit of the ECS compartment, we let C ECS be the j constant concentration of the species in the unit, 1 rj o N. We denote by Φj the volume-specific diffusion flux of a given species from the ECS in the jth subunit to the ECS in the (jþ 1)th subunit, and we model it by discrete Fick's law:

ℓj ¼

1

x ¼ n;a

Fig. 2shows the annular domains when N ¼5. If denote by V the volume of the Krogh cylinder, the volume of each subunit is

Sj ¼ 2π r j H;

S1

Φ0

We remark that the exchange rates may vary in time, although here for notational convenience we do not indicate explicitly the time dependence. Denoting the net exchange rate within the jth unit by X ECS-x x-ECS Jj ¼ ; Jj Jj

hj ¼ r j r j 1 ;

Φj ¼

Φ00

dC j DSj 1 ECS DSj ECS ¼ ðC C ECS ðC C ECS Þ Jj: j 1Þ þ j dt ℓj 1 j ℓj j þ 1

After rearranging the terms and substituting the geometric expressions, we arrive at the equation ECS dC rj 1 rj rj 2DN r ηECS j ¼ 2 j 1 C ECS þ þ C ECS ð22Þ C ECS j1 j j þ 1 ϕj ; dt ℓj 1 ℓj 1 ℓj ℓj R where

ϕj ¼

X Jj ¼ φECS-x φx-ECS j V u x ¼ a;n j

is the scaled net flux from the jth ECS compartment into the cells. Eq. (22) corresponds to a finite difference approximation of the radial diffusion model with local uptake by neurons and astrocytes. Similarly, for j¼ 1, we have the equation ECS

dC 1 ¼ Φ0 Φ1 J 1 dt DS1 ECS ¼ Φþ C 2 C ECS J1 ; 1 ℓ1

ηECS V u

and further, denoting ϕ0 ¼ Φ=V the uptake of the Krogh cylinder scaled by the volume

ηECS

ECS

dC 1 dt

¼ N ϕ0 þ

2DN r 1 ECS ϕ1 : C 2 C ECS 1 2 ℓ 1 R

ð23Þ


Finally, in the unit j¼N, we obtain the mass balance equation: ECS ECS dC N

η

dt

2DN r N 1 ECS ¼ 2 C N C ECS N 1 ϕN : R ℓN 1

ð24Þ

Using the notation ξj ¼ r j =ℓj , we introduce the finite difference matrix: 2

ξ1 6 6 ξ1 6 6 6 L¼6 6 6 6 4

3

ξ1 ðξ1 þ ξ2 Þ ξ2 ξ2 ðξ2 þ ξ3 Þ ξ3 ⋱

⋱ ðξN 2 þ ξN 1 Þ

ξN 1

ξN 1

ξN 1

7 7 7 7 7 7 A RNN 7 7 7 5

2

C ECS 1

3

6 7 7 C ECS ¼ 6 4 ⋮ 5; C ECS N

2

ϕ0

3

6 0 7 6 7 p ¼ N6 7; 4 ⋮ 5 0

2

ϕ1

ECS

dC dt

7 ϕ ¼6 4 ⋮ 5; ϕN

¼ κ LC ECS þp ϕ;

ηsc

ð25Þ

where κ ¼ 2DN=R2 . Fig. 2 illustrates the geometric setting when N ¼5. 3.2.3. Neuron and astrocyte In neuron and astrocyte, under the simplifying assumption of negligible diffusion, the concentration of substrates in the jth subunit changes as the results of chemical reactions, and exchange with the ECS. The mass balance equation, based on Eq. (6), is of the form

ηx V u

dC j X ¼ si Ψ i;j þJ ECS-x J x-ECS ; j j dt i

where x ¼n, a and the sum on the right hand side extends over all reactions taking place in the cell, Ψ i;j is the volume-specific reaction rate of the ith reaction and si is the stoichiometric constant. In this very simplified model, the coupling between the TCA cycle and the neurotransmitter cycle occurs only through the energy balance, and the carbon exchanges between the two cycles are not explicitly modeled. Each subunit is assigned the same metabolic pathways following Calvetti and Somersalo (2011), and 14 substances are tracked in each cell compartment in every subunit. The list of reactions considered is listed in Table 1. After scaling by the volume of the subunit, the mass balance equations are of the form

ηx

dC j X ¼ si ψ i;j þ φECS-x φx-ECS ; j j dt i

dC j ¼ φGlu;n-sc φGlu;sc-a ; j j dt

sc C j ¼ Glu j :

ð27Þ

where the fluxes have been scaled by the volume of the unit V u . The scaled mass balance equation for glutamine is of the similar form:

3

the equations governing the concentrations in the ECS subunits can be expressed in matrix form as

ηECS

to neighboring synapses. The removal relies mostly on the active glutamate uptake by the excitatory amino acid transporters (EAAT) in astrocytes. The glutamate uptaken by astrocytes will generate glutamine, which enters the synaptic cleft and eventually enters the glutamatergic neurons. These neurotransmitter reactions and transports are commonly referred to as the V-cycle. Our model only tracks the change in concentration of glutamate (Glu) and glutamine (Gln), and includes the transmembrane glutamate and glutamine transport. For more details about the model, see Calvetti and Somersalo (2011). The mass balance equation of glutamate in synaptic cleft (sc) with volume fraction ηsc in the jth unit is

ηsc

and the vectors:

55

ð26Þ

, and φx-ECS are the fluxes scaled by the volume where ψ i;j , φECS-x j j of the subunit. 3.2.4. Synaptic cleft The synaptic cleft is the small compartment where the neurotransmitters are exchanged. The glutamatergic neuron synthesizes the excitatory neurotransmitter glutamate from glutamine. The vesicular glutamate transporter (VGLUT) collects glutamate in the presynaptic neuron into vesicles which, in response to electropotential signalling, fuse with the presynaptic plasma membrane, leading to the secretion of the neurotransmitter into the synaptic cleft. Glutamate, which signals the neuronal activity to the post-synaptic neuron, needs to be removed rapidly from the cleft to prevent exitotoxicity and spillover

dC j ¼ φGln;a-sc φGln;sc-n ; j j dt

sc C j ¼ Gln j :

ð28Þ

To maintain mass balance, the corresponding fluxes need to be subtracted from the astrocyte and neuron equations (26) in the respective compartments. The V-cycle transports of glutamate and glutamine are assumed unidirectional, and we ignore the reuptake of glutamate by neurons, assuming that all glutamate goes to astrocyte. These simplifications give rise to a unidirectional transport flux expression: we use the model n o Cj sc φsc-x ¼ T max ; ðC j ; xÞ A ð½Glusc j j ; aÞ; ð½Glnj ; nÞ M þ Cj for the uptake fluxes from the cleft, and n o Cj φx-sc ¼ T max ; ðC j ; xÞ A ð½Glunj ; nÞ; ð½Glnaj ; aÞ : j M þ Cj To model the neuronal activity, the glutamate efflux from neuron is subject to a time dependent modification as described in Appendix A. 3.3. Building the system We are now ready to assemble the system of equations that constitute the full kinetic model in the discretized cylindrical geometry. Blood domain:. We follow four species in the blood: glucose (Glc), lactate (Lac), oxygen (O2 ) and carbon dioxide (CO2 ), each concentration satisfying Eq. (21) (4 equations). Extracellular space:. We follow the same four substance in ECS as in blood. Each one is distributed in N subdomains, and we have thus four systems given by Eq. (25), for a total of 4N equations. Neuron and astrocyte:. Both the neuron and astrocyte domains are further subdivided into N subdomains. In each cell type, we track the concentration of 14 species: glucose, lactate, oxygen, carbon dioxide, pyruvate (Pyr), adenosine triphosphate (ATP), adenoside diphosphate (ADP), inorganic phosphate (Pi), nicotinamide adenine dinucleotide (NADH and NAD þ ), creatine (Cr), phosphocreatine (PCr), glutamate, and glutamine. The concentrations are governed by Eq. (26), for a total of 2 14N ¼ 28N equations. Synaptic cleft:. Each subunit is equipped with its synaptic cleft, in which two species (glutamine and glutamate) are tracked. Their concentrations are governed by Eqs. (27) and (28), for a total of 2N equations. The full system consists therefore of a system of 34N þ4 differential equations, which are coupled by the coupling conditions. Further details concerning the coupling equations, reaction fluxes, time dependent excitation model, as well as the inputs of the system are given in Appendix A.

56


U(t) Activation function

Unit 1

Ca2+ +

Extracellular space

Unit 5

(t), M(t) Tmax(t) Glutamate transporter in neuron

Q(t) Blood flow

Neuron

Synaptic Gap Astrocyte

[Glu]Cleft

Na

ATPase-Neuron

V-cycle

Fig. 3. Left: schematics of the activation protocol. Right: a schematics showing the two locations of the activation in relation to the blood vessel, namely the unit closest to the blood vessel (Unit 1) and that farthest away from it (Unit 5). The smaller inset shows the interaction between compartments in each unit.

4. Implementation and simulations The goal of the numerical simulation with the radially distributed system derived above is twofold. First, we want to show that the fine scale dynamic behavior of the system is sensitive to the locus of the neuronal activity, in the sense that if the energy demand due to synaptic activity increases in a unit adjacent to a source of oxygen and glucose, the response of the system is different from when the activation takes place in a unit removed from the capillary. The second point that we want to make is that when fitting a spatially lumped model to average concentration data corresponding to the two different activation protocols, the lumped kinetic parameters may differ significantly from each other. This latter point provides a natural explanation why carefully calculated model parameters based on experiments may differ from experiment to another. Also, it gives a good reason to believe that a proper description of a kinetic parameter in a lumped model is a distribution rather than a single value. In order to perform simulations with the distributed model we need to assign the parameter values. The parameter selection is discussed in detail in Appendix B. 4.1. How the spatial location of the activation affects the dynamics of the metabolism In this section we address the question of whether and how changing the spatial location of neuronal activity with respect to the location of the blood vessel changes the landscape of the metabolic activity. In particular, a question which has received quite a lot of attention in recent years, and conflicting answers, is whether the neuron relies on glucose or lactate to meet the energetic demands associated with sustaining neuronal activation. To find out if the disagreement in the conclusions of different laboratories may be attributed to the lack of spatial resolution of the mathematical models, we perform a computed experiment where we generate two sets of data with the radially distributed model, corresponding to two different loci of the neuronal activity. More specifically, using a discretization of the volume into N ¼ 5 subcompartments, in the first simulation we activate the unit closest to the blood vessel, while in the second simulation the activity occurs in the unit furthest away from it. Fig. 3 indicates schematically how the activation protocol in a selected unit is computationally organized. Details are given in Appendix A. The figure also indicates the spatial arrangement of the two different locations of the activation, and one of the interactions between neuron, astrocyte and synaptic cleft in each

unit. We point out that the activation in a subunit represents the average neuronal activity integrated over the unit, and therefore no interpretation of what cell level activity configuration it represents can be given. This, of course, is true also for spatially lumped models. The general three-dimensional model allows, within the computational resolution, a more detailed description of the underlying cell level activity configurations.

4.1.1. Time course of substrate concentrations In the first numerical simulation, we examine the time course of the concentrations of some key substrates under the two activation patterns. As shown in the left-most panels in Fig. 4, regardless of the location of the activation, the glucose concentration in ECS drops at the start of the neuronal activity, and the decay is more pronounced and longer lasting when the activation is further from the blood vessel than when it is activated near it. More specifically, when the activation occurs in the first unit, the glucose concentration drops by 6% while when the fifth unit is activated it decreases by 16%. The concentration of glucose in neuron follows the same pattern (not shown). The change in lactate concentration, on the other hand, is remarkably different in the two cases. When the activation occurs in the first unit, the lactate concentration in ECS drops right at the beginning of the activity, rebounding up slightly at t¼6.5 min, before dropping again more significantly, reaching its minimum at t¼9 min and subsequently increasing again and reaching the initial concentration at t¼ 20 min, as shown in the second panel from the left in the top row of Fig. 4. When the activity occurs in the fifth unit, on the other hand, a very brief slight drop in the lactate concentration at the start of the activity in all units is followed by an increase of about 11% at t¼6 min at which time the concentration starts decreasing, reaching its minimum, below the steady state value, at t¼13.5 min, when it begins increasing slowly, returning to steady state value around t¼30 min. Summarizing, when the first unit is activated the lactate concentration drops at all locations, while when the activation takes place in the fifth unit, the lactate concentration increases substantially in all units, and it takes about 50% longer to have it back at the initial concentration level. The lactate concentrations in neuron and astrocyte follow the same pattern (data not shown). While the concentrations of glucose and lactate at steady state are close to each other in all units, the O2 concentration in ECS has a stronger radial gradient, and it increases immediately by approximately 11 % at all locations at the time of activation, returning to baseline value after 4–5 min when the unit nearest to blood vessel is activated. Note that the baseline concentration is highest in the unit


ECS

57

ECS

CGLC, AU=1

ECS

CLAC , AU=1

CO2 , AU=1

1 2

1.2

3 4

1.1

5 1

0

10

20

1 1.3

2 3

1.25

4 1.2

5

1.15 1.1

30

0

10

time(min) ECS , GLC

C

concentration(mmol/L)



1.35 1.3

20

2 0.04

3 4

0.03

5 0.02 0.01

30

1

0.05

0

10

ECS , LAC

AU=5

C

20

30

time(min)

time(min)

ECS , O2

AU=5

C

AU=5

1.3 1 2

1.2

3 4

1.1

1

5

0

10

20

30




1.35 1

1.3

2 3

1.25

4 1.2

5

1.15 1.1

0

10

time(min)

20 time(min)

30

1

0.05

2 0.04

3

0.03

4 5

0.02 0.01

0

10

20

30

time(min)

Fig. 4. Time course of glucose, lactate and oxygen concentrations in ECS in the units with the activation is in the first (top row) and the fifth (bottom row) unit. The neuronal activation effectuated by an activation onset at t¼ 3 min and with exponential decay is indicated by the background shading, see Appendix A for further details.

nearest to blood vessel, where it is nearly twice as large as in the unit farthest from it; see Fig. 4 right most top panel. The pattern is analogous when the fifth unit is activated, with the difference that before returning to baseline values there is a slight dip in concentration levels; see Fig. 4 right most bottom panel. The increase in oxygen concentration, in spite of the increased need for it, is due to the concomitant increase in oxygen supply by increased blood flow. The change in glucose concentration in astrocyte (see Fig. 5) at the beginning of neuronal activation is rather remarkable, dropping about 45% and returning to baseline values after t¼ 10 or t¼ 20 min, depending on whether the activity took place in the first or fifth unit. Glucose concentration in non-activated units exhibits a smaller and delayed drop. The concentration of O2 in neuron and astrocyte decreases only in the activated unit, while increasing in all the others. This ensures that diffusion brings O2 to the location which has highest energy demand. The decrease is three times higher when the activity takes place in the fifth unit rather than in the first, and it takes longer to return to baseline values; see Fig. 5. Baseline concentration of ATP is essentially the same in all units for both neuron and astrocyte, exhibiting a slight decrease in the activated unit (data not shown). The concentration of NADH in both neuron and astrocyte increases in the activated unit only, while decreasing in the other units. We remark that the baseline concentration is not uniform in all units, but in the unit farthest away from the blood vessel it is 2.5 times more than in the unit nearest to it, the opposite of what was the case for oxygen concentrations. When the first unit is activated, the concentration of NADH in neuron in the first unit increases only 1% while the concentration of NADH in neuron in the fifth unit increases 38% when the fifth unit is activated. The return to baseline values takes longer when the fifth unit is activated.

4.1.2. Time course of reaction fluxes and exchange rates We consider now the time courses of some of the fluxes in the system under the two different protocols. The left panels of Fig. 6

shows the glutamate efflux rate from the neurons. Since the glutamate efflux from neuron to the synaptic cleft is a proxy for the activation in the activation protocol, the results confirm that we assign the same activation to the first subunit and the fifth subunits, and no leak of synaptic activity to neighboring units takes place. The rate of glycolysis in the first unit when it is the site of the activity is roughly doubled compared to the steady state baseline, while in the non-active units, it is actually slightly less than at steady state. On the other hand, the increase in the glycolysis rate is fourfold in the fifth unit when the activity occurs there. This can be attributed to the fact that in the time it takes for the oxygen to reach the unit farthest away from the blood vessel, the energetic demand is met by the wasteful non-oxidative glycolysis rather than the more efficient oxidative phosphorylation. The increase in glycolysis leads also to a corresponding increase in glucose uptake by the neuron. In astrocyte, the glycolysis flux increases approximately to 2.5 fold in the activated unit, regardless of the location of the activation, see Fig. 7. When the fifth unit is the site of the activity, the time it takes for the flux to return to baseline value is nearly double than when the activity occurs in the first unit. This seems to suggest that with low oxygen levels, the clean-up of the glutamate from the active cleft requires more time. Since the glycolysis is related to the glucose transport from ECS to astrocyte, the transports follow the same patterns (not shown). In the light of the discussion about the primacy of glucose or lactate as a source of activated neurons (see, e.g., Wyss et al., 2011, DiNuzzo et al., 2010 and references therein), it is of interest to follow the lactate dehydrogenase (LDH) flux. With the sign convention in this work, positive LDH flux of this bidirectional reaction corresponds to oxidation of lactate to pyruvate, and negative sign is the pyruvate reduction to lactate. The LDH fluxes are shown in Fig. 8. Observe that in the simulated steady state, the neuronal LDH fluxes are positive, that is, neurons uptake small amounts of lactate. The figure shows that as soon as the unit next to the capillary is activated, the LDH flux of the first unit increases 21% at around t¼6 min, returning to its baseline value at t¼15 minutes, and in the

58


Astrocyte

CGLC

Neuron

, AU=1

CO2

1 2

0.6

3 4

0.4

5 0

10

20

30

1

0.04

2

0.03

3 4

0.02

5 0.01 0

40

0

5

15

2

0.03

3 4

0.02

5 0.01 0

20

0

5

10

15

time(min)

time(min)

CAstrocyte, AU=5

CNeuron, AU=5

CAstrocyte, AU=5

O2

1 2

0.6

3 4

0.4

5 10

20

30

0.05 concentration(mmol/L)


0.8

1

0.04

2

0.03

3 4

0.02

5 0.01 0

40

20

O2

0.05

1 concentration(mmol/L)

10

1

0.04

time(min)

GLC

0

, AU=1

0.05 concentration(mmol/L)



0.8

0.2

CO2

0.05

1

0.2

Astrocyte

, AU=1

0

5

10

15

0.04

1 2

0.03

3 0.02

4 5

0.01 0

20

0

5

time(min)

time(min)

10

15

20

time(min)

Fig. 5. Time courses of glucose concentrations in astrocyte and of oxygen concentration in neuron and astrocyte with the two activation patterns. In top row, the activation takes place in the first unit, in the bottom row in the fifth unit.

0.5

0

0

5

10 time(min)

15

20

0.8 1 2

0.6

3 0.4

4 5

0.2 0

0

5

1 2 3 4 5

1

0.5

0

0

5

10 time(min)

15

10 time(min)

15

20

20

1 2

0.6

3 0.4

4 5

0.2 0

0

2

0.6

3 4

0.4

5 0.2

5

10 time(min)

10 time(min)

15

20

ECS → Neuron

1

0

5

φGLC

0.8

0

, AU=1

0.8

ψNeuron , AU=5 Glycolysis

GLU

Distributed Fluxes(mM/min)


φNeuron → Cleft, AU=5 1.5


1 2 3 4 5

1

φGLC

Glycolysis


1.5

ECS → Neuron

ψNeuron , AU=1

, AU=1 Distributed Fluxes(mM/min)


Neuron → Cleft

φGLU

15

20

, AU=5

0.8 1 2

0.6

3 4

0.4

5 0.2 0

0

5

10

15

20

time(min)

Fig. 6. Left: time course of glutamate efflux of the 5 units with the two activation patterns. Center: the glycolytic activity of each unit. Right: the glucose uptake of neurons from the extracellular space. The upper row corresponds to activation of the innermost unit, the bottom row to activation of the unit 5 farthest from the capillary.

neurons in other units, the fluxes initially increase too. In contrast, when the activation occurs in the fifth unit, its neuronal LDH flux drops from 0.1746 mM/min to 0.5767 mM/min, indicating that the active unit reduces pyruvate to lactate, a sign of anaerobic glycolysis, returning to its baseline value at t¼10 min. Concurrently, the LDH in neuron of other units increases slightly, due to the increased availability of lactate. This behavior suggests that the location of the activation in

reference to the blood vessel has a dramatic effect on the production or consumption of lactate in neurons. Far from capillaries, lack of oxygen at the beginning of the activation requires non-oxidative metabolism and elevated glycolysis which leads to lactate production, while availability of oxygen makes lactate a viable fuel for neuronal activity. The transport rate of lactate from ECS to neuron follows the pattern of the LDH flux (not shown).


Astrocyte

ψGlycolysis, AU=5

Astrocyte



ψGlycolysis, AU=1 1 2

0.3

3 4

0.25

5 0.2 0.15

0

5

10

59

15

20

1 2

0.3

3 4

0.25

5 0.2 0.15

0

5

10

15

20

time(min)

time(min)

Fig. 7. The time course of glycolysis in astrocyte. Left: activation in the first, innermost unit. Right: activation in the unit farthest away from the capillary.

Neuron

Astrocyte

1

0.35

2 3

0.3

4 0.25

5

0.2 5

10

15

20

0.7

2

0.6

3 4

0.5

5

0.4 0.3

0

5

10

15

20

3.5 1 3

2 3

2.5

4 5

2 1.5

0

5

10

time(min)

time(min)

ψNeuron , AU=5 LDH

ψAstrocyte , AU=5 OxPhos

ψNeuron , AU=5 OxPhos

0.2 1

0

2 3

−0.2

4

−0.4

5 5

1

time(min)

0.4

0

0.8

10

15

20

0.8


0

ψOxPhos, AU=1 Distributed Fluxes(mM/min)


0.4

0.15

Neuron

ψOxPhos , AU=1

, AU=1




ψLDH

1

0.7

2 3

0.6

4 0.5

5

0.4 0.3

0

5

time(min)

10 time(min)

15

20

15

20

3.5 3

1 2

2.5

3 4

2 1.5

5

0

5

10

15

20

time(min)

Fig. 8. Left: time course of the LDH in neurons. Center: oxidative phosphorylation fluxes in astrocyte. Right: oxidative phosphorylation fluxes in neuron. The upper row corresponds to activation of the unit closest to capillary, the bottom row to activation of the unit farthest from it.

Summarizing, our spatially distributed model suggests that the Astrocyte-Neuron Lactate Shuttle (ANLS) can be observed when the neuronal activity occurs in proximity of a blood vessel, whereas both neuron and astrocyte produce lactate when the fifth unit is activated, confirming a strong role of the activation location on the preferred substrates sustaining neuronal activity. Regardless of the locus of the activity, the oxidative phosphorylation flux in neuron increases significantly, between 66% and 79% in the unit where the activation takes places, and slightly in all others, returning to baseline in 10–12 min after the beginning of the activation. The flux increase is more pronounced when the activation is near the blood vessel, due to better availability of oxygen, but its returns to the baseline slower, as can be seen in Fig. 8. The same pattern can be seen in astrocyte, where the proximity of the blood vessel seems to have a stronger effect. In fact, when the locus of the activation is near the blood supply, the oxidative phosphorylation flux in the astrocyte of the first unit increases by 90% and returns to its baseline value at t¼12 min, while when the fifth unit is activated, its oxidative phosphorylation flux in astrocyte increases 47% and returns to its baseline value also at t¼12 min.

Since the pyruvate dehydrogenate (PDH) and tricarboxylic acid cycle (TCA) are lumped into a single reaction flux ψ pdh , and the oxidative phosphorylation (OxPhos) flux φOxPhos are related by φpdh ¼ ð1=3ÞφOxPhos , the reaction fluxes of PDH/TCA (not shown) follow the same pattern of the oxidative phosphorylation, with the values multiplied by a factor of three. In neuron, ATP hydrolysis increases significantly in the units where the activation takes place, while remaining constant in all other units. This is in accordance with out activation patterns since only one unit is activated in the activation scenarios. For astrocyte, the ATP hydrolysis is essentially constant in all units. We can summarized the findings in the following main points: 1. The constant need of oxygen at steady state requires a significant oxygen gradient across the tissue, implying that less oxygen is available far from the capillary. 2. If neuronal activation takes place close to the capillary, the energy need is met with increased oxidative phosphorylation, which may use both glucose and lactate as fuel. For neuronal activity far away from the capillary, the oxygen level is so low that non-oxidative

60


metabolism is recruited to meet the energetic need, leading to a concomitant lactate production by neuron. 3. The oxygen replacement due to increased blood flow leads a rapid rise in oxygen levels throughout the ECS during activation, but it is not sufficient to sustain oxidative energy production during the distal neuronal activity.

To elucidate the latter point further, we perform a controlled experiment in which we generate two sets of data with the spatially distributed model with the two activation locations described above. Subsequently, mimicking the procedure followed when experimental data is available, we fit a spatially lumped parametric mean model to the two data sets and compare the outcomes.

4.2. From distributed to lumped: the effect of resolution on parameter estimation

4.2.1. Estimating the parameters of the lumped model For the sake of definiteness, to explain how we estimate the parameters of the lumped model, we consider the lumped glycolysis in neuron:

The problem of establishing a connection between parameters of models at different scales is a difficult and important one. In this section, we study the connection between the spatially lumped and the spatially distributed model, and investigate how the respective parameter can be related. The computational dynamic models of cellular brain energy metabolism which have appeared in the literature in recent years, see, e.g., Aubert and Costalat (2005), Calvetti et al. (2007), Cloutier et al. (2009), DiNuzzo et al. (2010) and Simpson et al. (2007), are spatially lumped. While there is no question about their contribution to advancing our understanding of brain energetics, it is puzzling that different groups have come to discording conclusions when using mathematical models together with experimental data. Possible explanations for the discrepancy include the different details incorporated in the models, differences in values assigned to the model parameters, and, motivational for the present work, the lack of the spatial resolution in lumped models, which does not allow us to take into account the effects of the location of the neuronal activity in reference to the blood vessel.

Table 2 Glycolysis: values of the parameters of the spatially distributed model used for data generation (True); estimates of the parameters for the lumped model using data with activation in the innermost compartment adjacent to the capillary (Near) and far away from it (Far), with the percentage change from the corresponding values for the distributed model. μ

ν

K

Neuron True 19.7561 Near 19.8481 (0%) Far 2.9587 ( 85%)

0.09 0.1856 (106%) 0.0134 ( 85%)

10 16.5276 ( þ 65%) 5.6471 ( 43%)

4.4656 1.1060 ( 75%) 3.7756 ( 15%)

Astrocyte True 19.5895 Near 3.7168 ( 81%) Far 5.9501 ( 70%)

0.09 0.0195 ( 78%) 0.0088 ( 90%)

10 9.6478 ( 3%) 4.2474 ( 58%)

2.4927 1.5831 ( 36%) 7.2667 ( þ191%)

V max

Glc þ 2 ADP þ 2 NAD þ þ 2 Pi-2 Pyr þ 2 ATP þ 2 NADH: The reaction is oxidative, involving the phosphorylation of ADP, and therefore the Michaelis–Menten form for the reaction flux in this case is

ψ ¼ V max

1 1 r phos r redox ½Glc ; 1 νþr 1 K þ ½Glc μ þ r phos redox

which we can express in the form 2 3 ½Glc 2 3 V max 6 ½ATP 7 6 7 6 7 6 7 6 K 7 7 ψ ¼ f ðC; θÞ; C ¼ 6 7: 6 ½ADP 7; θ ¼ 6 4 μ 5 6 ½NAD þ 7 4 5 ν ½NADH Consider now a simulation of the discretized, radially distributed model, which, with a given activation pattern, produces a vector of glycolysis fluxes in neuron at given time instances tk. The computed

Table 3 Oxidative phosphorylation: values of the Michaelis–Menten parameter of the spatially distributed model (True) used for data generation, and estimates of the parameters for the lumped model using data with activation in proximity of the blood vessel (Near) and far away from it (Far), with the percentage change from the corresponding values for the distributed model. μ

ν

K

Neuron True 700.1367 Near 700.1367 (0%) Far 700.1363 (0%)

0.1 0.0515 ( 48%) 0.3032 ( þ203%)

0.01 0.1041 (þ 941%) 0.0002 ( 98%)

0.8574 0.1894 ( 78%) 1.1680 ( þ36%)

Astrocyte True 157.0148 Near 157.0155 (0%) Far 157.0137 (0%)

0.1 0.0090 ( 61%) 0.3376 ( þ 237%)

0.01 0.6482 ( þ6382%) 0.0003 ( 97%)

0.8574 0.0466 ( 95%) 1.0543 ( þ 23%)

V max

Neuron

Astrocyte

ΨGlycolysis

Curve: act. unit 1 Data: act. unit 1 Curve: act. unit 5 Data: act. unit 5

0.3 0.25 0.2

0

10

20 time(min)

30

40

Averaged fluxes(mM/min)

Averaged fluxes(mM/min)

ΨGlycolysis 0.35

ð29Þ

0.24 Curve: act. unit 1 Data: act. unit 1 0.22

Curve: act. unit 5 Data: act. unit 5

0.2

0.18

0

10

20 time(min)

30

40

Fig. 9. The dashed curves show the time course of the glycolysis in neuron (left) and in astrocyte (right) using the estimated lumped parameter values. The simulated data is plotted with dots. Blue dots and blue curve correspond to activation in the unit nearest to the capillary, red ones to the activation of the unit farthest away from the capillary. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)


61

Table 4 Michaelis–Menten parameter values for the reactions in neuron and astrocyte. Most of these values are obtained starting from the corresponding values in Calvetti and Somersalo (2011) and tuning them so that the system supports a stable steady state. Reaction

V max ðmM=minÞ

K ðmMÞ

μ

ν

Neuron

Glycolysis Lactate formation Lactate reduction PDHþ TCA PCr - Cr Cr - P Cr OxPhos ATP hydrolysis Gln - Glu

19.76 8.616e04 9.479e04 2.57 1e06 1e06 700.1 25.94 1.23

4.47 2.15 23.70 0.0125 528 495 0.8574 2.2 0.003

0.09 – – 0.01 100 0.01 0.01 – –

10 0.1 10 10 – – 0.1 – –

Astrocyte

Glycolysis Lactate formation Lactate reduction PDHþ TCA PCr - Cr Cr - P Cr OxPhos ATP hydrolysis Glu - Gln

19.59 2.496e05 1.947e05 0.576 1e06 1e06 157.01 3.12 2.46

2.49 6.24 48.67 0.0124 528 495 0.8574 2.2 0.03

0.09 – – 0.01 100 0.01 0.01 – 100

10 0.1 10 10 – – 0.1 – –

Table 5 Parameter values for the passive and facilitated exchange rates. Most of these values are obtained starting from the corresponding values in Calvetti and Somersalo (2011) and tuning them so that the system supports a stable steady state. Exchange

T max ðmM=minÞ

M ðmMÞ

λ ð1=minÞ

BBB

Glucose (GLUT-3) Lactate (MCT) O2 CO2

1.33 10 – –

4.47 5 – –

– – 463.07 0.0355

Neuron

Glucose (GLUT-3) Lactate (MCT) O2 CO2 Glutamate Glutamine

6667 4000 – – 2.46 2.46

6.667 0.4 – – 98 7e 05

– – 304.3 0.9948 – –

Astrocyte

Glucose (GLUT-1) Lactate O2 CO2 Glutamate Glutamine

667 4000 – – 2.46 2.46

1667 0.4 – – 7e-5 0.07

– – 68.35 0.22 – –

Table 6 Input values for the system, specifying the volume fractions of the compartments, blood and blood flow parameters, and diffusion coefficient. Volume fractions

Flow parameters

Blood parameters

Diffusion coefficients

Neuron

0.45

Q (mL/min)

0.5

Hct

0.45

Dglc

4:4 10 6 cm2 =s

Astrocyte

0.25

F

2/3

Hbrbc

5.18

Dlac

4:4 10 6 cm2 =s

ECS

0.20

KH (mM)

26 1:4e 3

D O2

1:7 10 5 cm2 =s

Cleft

0.01

K1 (mmol)

45.9

DCO2

Blood

0.04

K2 (mmol) K3 (mmol) K4 (mmol)

23.9 243.1 1.52

rKrogh

1:7 10 5 cm2 =s 50 μm

output therefore consists of the pairs of concentration vectors and reaction fluxes in each of the N units:

ψ j ðt k Þ;

C j ðt k Þ;

1 r j rN; 1 rk rK;

the subindex referring to the jth neuron compartment. The averaged flux and concentrations in the tissue sample, based on the

simulations, are

ψ ðt k Þ ¼

N 1X ψ ðt Þ; Nj¼1 j k

C ðt k Þ ¼

N 1X C ðt Þ: Nj¼1 j k

Assume that the lumped model also uses a Michaelis–Menten model for the glycolytic flux in neuron, and let θ A R4 denote the

62


Table 7 Steady state concentration values for the tracked substrates and intermediates. Most of these values are obtained starting from the corresponding values in Calvetti and Somersalo (2011) for the single unit case and tuning them so that the system supports a stable steady state. Com.

Concentration (mmol=L) Neuron

Astrocyte

Unit

1

2

3

4

5

1

2

3

4

5

Glc Pyr Lac O2 CO2 Pi PCr Cr Glu Gln ATP ADP NAD þ NADH

1.2585 0.5528 1.1751 0.0378 88.0628 12.8301 7.1671 3.1663 14.0020 0.0010 2.2028 0.0092 0.0322 0.0008

1.2437 0.2295 1.1775 0.0169 88.0752 13.6493 6.3517 3.9817 14.0014 0.0010 2.1990 0.0130 0.0313 0.0017

1.2404 0.1668 1.1782 0.0123 88.0729 14.1870 5.8159 4.5175 14.0010 0.0010 2.1960 0.0160 0.0308 0.0022

1.2390 0.1415 1.1786 0.0104 88.0696 14.5459 5.4607 4.8727 14.0007 0.0010 2.1936 0.0184 0.0304 0.0026

1.2385 0.1324 1.1788 0.0097 88.0677 14.7088 5.2988 5.0345 14.0005 0.0100 2.1924 0.0196 0.0303 0.0027

0.8219 0.4872 1.1759 0.0378 88.0787 12.8667 7.1306 3.2027 0.0065 0.0100 2.2026 0.0094 0.0322 0.0008

0.7763 0.2206 1.1783 0.0169 88.0778 13.6600 6.3409 3.9924 0.0075 0.0100 2.1990 0.0130 0.0313 0.0017

0.7499 0.1635 1.1790 0.0123 88.0660 14.1676 5.8349 4.4984 0.0082 0.0100 2.1961 0.0159 0.0308 0.0022

0.7323 0.1400 1.1793 0.0104 88.0564 14.5028 5.5058 4.8275 0.0089 0.0100 2.1939 0.0181 0.0304 0.0026

0.7242 0.1314 1.1795 0.0097 88.0517 14.6495 5.3575 4.9758 0.0092 0.0100 2.1929 0.0191 0.0303 0.0027

Com. Unit

ECS 1

Glc Lac O2 CO2

1.2587 1.1755 0.0436 86.2856

Blood 2

3

1.2440 1.1779 0.0226 86.3066

1.2407 1.1786 0.0180 86.3112

4 1.2393 1.1789 0.0161 86.3131

5 1.2388 1.1791 0.0154 86.3138

corresponding model parameters of the lumped model. Define the residual function between the model and the averaged data: Rðθ Þ ¼

K X

2

ψ ðt k Þ f ðC ðt k Þ; θ Þ ;

k¼1

and let the optimal lumped model parameter vector θ opt satisfy n o θ opt ¼ arg min Rðθ Þ subject to ðθ ℓ 4 0Þ: We estimate the optimal lumped parameters using a standard Levenberg–Marquardt least squares method; see Dennis and Schnabel (1996). Moreover, since the parameters are non-negative, we will project the negative values into the non-negative cone during the optimization process. If the algorithm fails to find a optimized parameter set, we take it as a sign that the current model cannot capture the averaged simulation data, that is the data gained from the mean approximation of the distributed model. Below, we present selected results obtained when estimating the parameters for the lumped model from the data set where the activation is near the capillary and from the data set where the activation is farthest from the blood supply. Consider first the estimation of the Michaelis–Menten parameters of the glycolysis, as discussed in the previous section. Table 2 shows the estimated lumped parameters corresponding to glycolysis both in astrocyte and in neuron, using the synthetic data corresponding to activation near the capillary and far away from it. We also indicate how much the estimated values differ from the underlying values used to generate the data. To assess how well the estimated parameters capture the data, we show in Fig. 9 the time courses of glycolysis predicted by the lumped model with the two sets of estimated parameters, as well as the mean approximation data used in the parameter estimation process. In neuron, the estimate of V max using the simulation data from the distributed model when the activity takes place in the first unit is almost the same value as in the distributed model, while its estimate using the simulation data from the distributed model when the activity takes place in the fifth unit is 85% smaller. The estimated μ is 106% larger when the first unit is activated and 85 % smaller when the fifth

4.4988 1.1481 6.2903 25.8603

unit is activated, while ν is 65% larger when the first unit is activated and 43% smaller when the fifth unit is activated. The estimates for K are approximately 75% smaller when the first unit is activated and the decrease is 15% when the fifth unit is activated. In astrocyte, the estimates for V max are 81% and 70% smaller than in the distributed model in correspondence of activity in the first and the fifth unit, respectively, and the decrease in the estimated μ is similar. The estimate of ν is slightly smaller when the activity occurs near the blood vessel, and approximately 58% smaller when it is far away. The estimates of K are 36% smaller than in the distributed model when the first unit is activated, while it is 191% larger when the fifth unit is activated. Similar to the glycolysis, we consider the parameters related to the oxidative phosphorylation: O2 þ 2 NADH þ 5 Pi þ 5 ADP-2 NAD þ þ 5 ATP þ 2 H2 O; which is modeled with a Michaelis–Menten form similar to the glycolysis. The results of the estimation both in neuron and astrocyte are given in Table 3. We see that similar to what happens for the glycolysis, the deviations, in particular for the affinity parameters of the cofactors (μ and ν ), are quite significant. Summarizing, we conclude that without some prior information about the parameters, the estimation of lumped parameters by mere optimization is highly sensitive to the underlying unresolved fine-scale details. This finding highlights some of the problems that could result when using in lumped models of neuron–astrocyte interactions parameter values validated on patch of the cellular membrane. In summary, parameter values that are reasonable at a given scale may no longer be valid at a different one.

5. Discussion and outlook In this paper, we introduce a novel multi-domain formalism for three dimensional distributed system of cellular metabolism using reaction–diffusion equations in each of the domains, with the possibility of including additional details, for example separate mitochondria. This formalism allows us to track the biochemistry


Glucose

63

Lactate

0.25

0.4 Neuron Astrocyte

0.2

ECS to Tissue[mM/min]


0.3


0.1

0.05

0.2 0.1 0 −0.1 −0.2 −0.3

0

1

2

3

4

−0.4

5

1

2

Unit index O2



5

0


1

0.5

1

4

CO2

2

0

3 Unit index

2

3

4

5

−0.5

Neuron −1

Astrocyte

−1.5

−2

1

2

Unit index

3

4

5

Unit index

Fig. 10. The rate of glucose uptake in subunits is 0.3280, 0.3589, 0.3830, 0.4006 and 0.4091 mM/min, yielding an average uptake rate of 0.3759 mM/min; the rate of oxygen uptake in subunits is 2.1682, 2.1548, 2.1442, 2.1365 and 2.1327 mM/min, yielding an average uptake rate of 2.1473 mM/min; the rate of neuronal lactate uptake in subunits is 0.2823, 0.2424, 0.2103, 0.1864 and 0.1746 mM/min and the astrocyte lactate production in subunits is 0.2156, 0.2419, 0.2617, 0.2754 and 0.2819, yielding an average uptake rate of 0.0361 mM/min. The tissue uptake of glucose increases in the units further away from the blood while the uptake of oxygen decreases.

in mutually interacting domains without the need of detailed geometric modeling of the microscopic tissue structure. Due to the complexity of the multidimensional model and the lack of parameter values, we perform model reductions to lower the dimensionality, and making use of the available parameter values in the literature, we show simulated results in a simplified distributed model of a single Krogh cylinder. This model is used to elucidate the problems of drawing conclusions about the underlying details from a lumped model that lacks the resolution of finescale details. The results of our simulations suggest that well-mixed models that ignore diffusion may go astray in predictions that involve oxygen. More than the availability of glucose, the availability of oxygen seems to affect strongly the local behavior of the system. Modeling in life sciences requires an understanding of the multiscale nature of living systems, and it is important to connect data and parametric models in different scales so as to understand the limitations of models, and avoid erroneous extrapolations across the scales. The computed results reveal a significant sensitivity of parameters to scale change, which is an indication that simple parameter estimation methods by model fitting may be inadequate, and more sophisticated techniques using prior information need to be considered. The simulations may also help us to understand why experimental data collected under similar conditions may give rise to different conclusions, if the applied model lacks resolution. The simulations seem to reinforce the idea that a proper description of biological model parameters is not a single value, but rather a distribution of values (Somersalo et al., 2012). The numerical experiments in the current paper deal with a simplified and rather crude distributed model. Future directions include an implementation of a more general distributed model in two or in three dimensions. Possible model refinements include the investigation of the importance of direct metabolite uptake of

astrocytes through the end feet, and inclusion of diffusion in neurons and astrocytes, which may require a modeling of anomalous diffusion using fractional calculus rather than the standard diffusion model, see, e.g. Holcman and Schuss (2014) for a recent review of the theoretical aspects, and de Santis et al. (2011) and Özarslan et al. (2006) for experimental results with magnetic resonance imaging. Furthermore, the metabolic model allows refining, e.g., by the inclusion of separate mitochondria, as well as by the enrichment of the spectrum of the substrates that the cells may use for energy metabolism. The computational complexity, however, makes these extension quite challenging. Also, the introduction of more capillaries will give us more insights about how the capillary distribution and strength affect the system's behavior (Wang and Bassingthwaighte, 2001).

Acknowledgements The work of Daniela Calvetti was partially supported by a Grant from the Simons Foundation (#246665 to Daniela Calvetti). The work of Erkki Somersalo was partially supported by NSF-DMS Grant 1312424.

Appendix A This appendix contains the details about the coupling of different compartment equations in the radially distributed model. In particular, we give the parameter values defining the system of differential equations, the input parameters, as well as the time dependent activation model.

64


As explained in Section 3.3, the dynamic model consists of a system of 34N þ4 differential equations, written as dC ¼ FðC; θ; C a Þ; dt

CðtÞ ¼ C 0 ;

A.1. Activation pattern The computational model simulates an increased energy demand in the neuron–astrocyte system due to synaptic activity over a specified synaptic cleft. We assume that in connection with the neuronal activity, three phenomena will occur: (a) Glutamate is released from neuron into the synaptic cleft, (b) a significant amount of ATP is dephosphorylazed to ATP to restore the membrane potential due to an increased ion influx in neuron, and (c) the capillary blood flow increases as a result of neurovascular signalling. Following Calvetti and Somersalo (2011), we simulate a neuronal activity starting at time t ¼ t 1 and fading exponentially, using the use a time-dependent dimensionless activity function: (

0;

t o t1 ;

aðt t 1 Þexpð γ ðt t 1 ÞÞ;

t Z t1 ;

where t 1 ¼ 3 min, γ ¼ 1 min 1 , and a ¼1.5. To capture the physiological changes which are known to occur in conjunction with neuronal activation, following Calvetti and Somersalo (2011), we integrate u(t) into our model so that (i) The blood flow increases with u(t), according to the rule Q ðtÞ ¼ Q low þ uðtÞΔQ ;

ΔQ ¼ Q low ;

(ii) Elevated values of u(t) upregulate the transfer of glutamate from neuron to the synaptic cleft, to capture the efflux of glutamate by the presynaptic neuron. If φGlu;n-sc ðtÞ denotes the transport rate of glutamate from the neuron to the synaptic cleft

φGlu;n-sc ðtÞ ¼ TðtÞ

½Glun ðtÞ ; MðtÞ þ ½Glun ðtÞ

where TðtÞ ¼ T 0 ð1 þ 2uðtÞÞ;

f ð½Gluc Þ ¼

ð30Þ

where C ¼ CðtÞ A R34N þ 4 contains all the concentrations, C a A R4 are the arterial values given in Table 6, C 0 A R34N þ 4 are the initial values listed in Table 7, and θ is a vector containing the model parameters. In Tables 4 and 5, we list the model parameters related to reaction fluxes and exchange rates. The autonomous system (30) is used to find an initial steady state: the initial values in Table 7 are found as explained in Appendix B. We simulate the neuronal activity introducing a time dependent forcing, which we implement by introducing time dependency to some of the parameters, that is, θ ¼ θðtÞ. The details are given in the following section.

uðtÞ ¼

where f is a sigmoidal saturation function of the form

MðtÞ ¼ M 0 ð1 uðtÞ=2Þ;

in which T0 and M0 represent the values corresponding to low activity. (iii) A rise in the glutamate concentration in the synaptic cleft increases the ATP hydrolysis in neuron, simulating the response of postsynaptic neuron to calcium and sodium influx through glutamate sensitive ion channels. Denoting by ψ ATPase ðtÞ the ATP hydrolysis flux in neuron, we let

ψ ATPase ðtÞ ¼ ψ base þ V max f ð½Gluc Þ

½ATPn ; K þ ½ATPn

½Glu2c kGlu þ ½Glu2c

:

(iv) The increased glutamate concentration in the synaptic cleft increases glutamate uptakes by astrocyte, where it is transformed into glutamine. Since this process requires ATP, the energetic cost in the astrocyte increases. Fig. 3 summarizes all the different aspects of the activation protocol.

Appendix B. Parameters of simulation with the radially symmetric model In the following, we discuss the parameters used in the computed examples with the radially symmetric model. The parameter values correspond as much as possible to a model of the human brain, see Gjedde (2007) and the discussion in Calvetti and Somersalo (2011). Since all the subunits are assumed to be physiologically similar with the same volume, we assign the same values for the diffusion coefficients, membrane permeabilities, and the Michaelis–Menten parameters for each subunit. Diffusion coefficients: In the radially distributed model, only the ECS domain was assumed to be diffusive. Following Nicholson and Syková (2008), we express the coefficients which regulates molecular self-diffusion in the brain as D ¼ Dfree =λ, where Dfree is the free diffusion coefficients and λ is the tortuosity, which, for glucose, has been experimentally measured as λ ¼ 1:6, see Simpson et al. (2007). We use this tortuosity value also for the gases O2 and CO2 . The diffusion parameter values, listed in Table 6, are taken from Simpson et al. (2007) for glucose and lactate, and from Kasischke et al. (2011) for the gases. Initial steady state: We start by determining a reasonable initial state that corresponds to a steady state. For the sake of definiteness, consider first oxygen in ECS. A literature value for oxygen uptake in the cerebral cortex, or cerebral metabolic rate (CMRO2 ) (Gjedde, 2007), is CMRO2 ¼ 1:6 μmol=min per gram tissue: The total flux from blood to Krogh cylinder is then Φ ¼ ðV=V 0 ÞCMRO2 , where V 0 is the volume of one gram tissue. By further expressing the volume of a unit in terms of one gram tissue, V u =V 0 ¼ N 1 ðV=V 0 Þ, we conclude that in (25) we may use the value

ϕ0 ¼ CMRO2 : Assume, a priori, that the oxygen is consumed at the same rate by each N units, so that the scaled flux leaving the jth ECS compartment to neuron and astrocyte is 1 N

ϕj ¼ ϕ0 ;

1 r jr N:

Consider now the system of Eq. (25). By fixing the flux values as above, and starting from an average oxygen concentration of the lumped model (Calvetti and Somersalo, 2011), we drive the system into a steady state, i.e., we find oxygen concentrations in the subunits that are capable of sustaining the steady state. In the same way, we find the steady state values for all four species (Glc, Lac, O2 and CO2 ) that our model tracks in the ECS. To find initial concentrations for the four substrates in neuron and astrocyte, we fix the exchange flux parameters across the cell membranes equal to the values suggested for the lumped model in Calvetti and Somersalo (2011), after which we adjust the


concentrations in the cell compartments so that they are capable to maintain the steady state. The initial values of the concentrations in cells of the biochemical species that are not exchanged between ECS and cell compartment are set equal to those in the lumped model steady state. The Michaelis–Menten parameters for the reactions are set as in Calvetti and Somersalo (2011). Since the concentrations of the exchangeable substrates are adjusted so as to sustain the prescribed transfer fluxes, there is no guarantee that the system is a steady state. Therefore, we start the system (30) with the values found so far and we run it for a long time until it reaches a steady state. The resulting steady state concentration values are listed in Table 7. The oxygen concentration drops significantly with the distance of the unit from blood while the concentration of glucose, lactate and carbon dioxide are pretty stable in all units. Each ECS subunit communicates with its neuron and astrocyte unit, and the rates of which the two cells uptake or release the four substrates in each unit are shown in Fig. 10. While the exchange of O2 and CO2 does not depend on the distance of the unit from the capillary, units further away from it uptake glucose at high rate. Moreover, the rate of lactate uptake of neuron decreases with the distance from the capillary, while the rate at which lactate is released from astrocyte increases. References Aubert, A., Costalat, R., 2002. A model of the coupling between brain electrical activity, metabolism and hemodynamics: application to the interpretation of functional neuroimaging. NeuroImage 17, 1162–1181. Aubert, A., Costalat, R., Magistretti, P.J., Pellerin, L., 2005. Brain lactate kinetics: modeling evidence for neuronal lactate uptake upon activation. Proc. Natl. Acad. Sci 102, 16448–16453. Boron, W.F., Boulpaep, E.L., 2012. Medical Physiology: A Cellular and Molecular Approach. Elsevier Saunders, Philadelphia, PA. Calvetti, D., Somersalo, E., 2011. Dynamic activation model for a glutamatergic neurovascular unit. J. Theor. Biol. 274, 12–29. Calvetti, D., Heino, J., Somersalo, E., Tunyan, K., 2007. Bayesian stationary state flux balance analysis for a skeletal muscle metabolic model. Inverse Probl. Imaging 1, 247–263. Calvetti, D., Somersalo, E., 2012. Ménage à trois: the role of neurotransmitters in the energy metabolism of astrocytes, glutamatergic, and GABAergic neurons. J. Cereb. Blood Flow Metab. 32, 1472–1483. Calvetti, D., Somersalo, E., 2013. Quantitative in silico analysis of neurotransmitter pathways under steady state conditions. Front. Endocrinol. 8http://dxdoi.org/ 10.3389/fendo.2013.00137. Cloutier, M., Bolger, F.B., Lowry, J.P., Wellstead, P., 2009. An integrative dynamic model of brain energy metabolism using in vivo neurochemical measurements. J. Comput. Neurosci. 27, 391–414.

65

Colli Franzone, P., Pavarino, L.F., Taccardi, B., 2005. Simulating patterns of excitation, repolarization and action potential duration with cardiac bidomain and monodomain models. Math. Biosci. 197, 35–66. Dennis, J.E., Schnabel, R.B., 1996. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. SIAM, Philadelphia. de Santis, S., Gabrielli, A., Bozzali, M., Maraviglia, B., Macaluso, E., Capuani, S., 2011. Anisotropic anomalous diffusion assessed in the human brain by scalar invariant indices. Magn. Res. Med. 65, 1043–1052. DiNuzzo, M., Mangia, S., Maraviglia, B., Giove, F., 2010. Changes in glucose uptake rather than lactate shuttle take center stage in subserving neuroenergetics: evidence from mathematical modeling. J. Cereb. Blood Flow Metab. 12, 586–602. Gjedde, A., 2007. Coupling of brain function to metabolism: evaluation of energy requirements. In: Handbook of Neurochemistry and Molecular Neurobiology. Springer Verlag, Berlin Heidelberg. Gruetter, R., Seaquist, E.R., Ugurbil, K., 2001. A mathematical model of compartmentalized neurotransmitter metabolism in the human brain. Am. J. Physiol. Endocrinol. Metab 281, E100–E112. Henriquez, C., 1993. Simulating the electrical behavior of cardiac tissue using the bidomain model. Crit. Rev. Biomed. Eng. 21, 177. Holcman, D., Schuss, Z., 2014. Time scale of diffusion in molecular and cellular biology. J. Phys. A: Math. Theor. 47, 173001. Kasischke, K.A., Lambert, E.M., Panepento, B., Sun, A., Gelbard, H.A., Burgess, R.W., Foster, T.H., Nedergaard, M., 2011. Two-photon NADH imaging exposes boundaries of oxygen diffusion in cortical vascular supply regions. J. Cereb. Blood Flow Metab. 31, 68–81. Kreft, M., Bak, L.K., Waagepetersen, H.S., Schousboe, A., 2012. Aspects of astrocyte energy metabolism, amino acid neurotransmitter homoeostasis and metabolic compartmentation. ASN Neuro 4, http://dx.doi.org/10.1042/AN20120007. Krogh, A., 1919. The number and distribution of capillaries in muscles with calculations of the oxygen pressure head necessary for supplying the tissue. J. Physiol. 52, 409–415. Nicholson, C., Syková, E., 2008. Diffusion in brain extracellular space. Physiol. Rev. 88, 1277–1340. Occhipinti, R., Somersalo, E., Calvetti, D., 2008. Astrocytes as the glucose shunt for glutamatergic neurons at high activity: an in silico study. J. Neurophysiol. 101, 2516–2527. Occhipinti, R., Somersalo, E., Calvetti, D., 2010. Energetics of inhibition: insights with a computational model of the human GABAergic neuron-astrocyte cellular complex. J. Cereb. Blood Flow Metab. 30, 1834–1846. Özarslan, E., Basser, P.J., Shepherd, T.J., Thelwall, P.E., Vemuri, B.C., Blackband, S.J., 2006. Observation of anomalous diffusion in excised tissue by characterizing the diffusion-time dependence of the MR signal. J. Magn. Res 183, 315–323. Roth, B.J., Wikswo, J.P., 1994. Electrical stimulation of cardiac tissue: a bidomain model with active membrane properties. IEEE Trans. Biomed. Eng. 41, 232–240. Simpson, I.A., Carruthers, A., Vannucci, S.J., 2007. Supply and demand in cerebral energy metabolism: the role of nutrient transporters. J. Cereb. Blood Flow Metab. 27, 1766–1791. Somersalo, E., Cheng, Y., Calvetti, D., 2012. The metabolism of neurons and astrocyte through mathematical models. Ann. Biomed. Eng. 40, 2328–2344. Wang, C.Y., Bassingthwaighte, J.B., 2001. Capillary supply regions. Math. Biosci. 173, 103–114. Wyss, M.T., Jolivet, R., Buck, A., Magistretti, P.J., Weber, B., 2011. In vivo evidence for lactate as a neuronal energy source. J. Neurosci. 31, 7477–7485.

A computational model of liver iron metabolism.

Computational Model Predicts the Effects of Targeting Cellular Metabolism in Pancreatic Cancer.

A cellular perspective on brain energy metabolism and functional imaging.

Spatially distributed dendritic resonance selectively filters synaptic input.

A Bayesian nonparametric model for spatially distributed multivariate binary data with application to a multidrug-resistant tuberculosis (MDR-TB) study.

Spatially variable coevolution between a haemosporidian parasite and the MHC of a widely distributed passerine.

Cellular compartmentation of brain metabolism and its functional significance.

The spread of a parasitic infection in a spatially distributed host population.

Ciliopathy is differentially distributed in the brain of a Bardet-Biedl syndrome mouse model.

Pain: a distributed brain information network?

A Bayesian probit model with spatially varying coefficients for brain decoding using fMRI data.

Erratum to: Bacteriophages affect evolution of bacterial communities in spatially distributed habitats: a simulation study.

Bacteriophages affect evolution of bacterial communities in spatially distributed habitats: a simulation study.

Modeling evolution of spatially distributed bacterial communities: a simulation with the haploid evolutionary constructor.

Modeling the isotopic evolution of snowpack and snowmelt: Testing a spatially distributed parsimonious approach.

Electromechanics of a membrane with spatially distributed fixed charges: flexoelectricity and elastic parameters.

Regional brain metabolism in a murine systemic lupus erythematosus model.

A distributed pore size model for capillaries.

A Computational Model of Hepatic Energy Metabolism: Understanding Zonated Damage and Steatosis in NAFLD.

Power output and fatigue properties using spatially distributed sequential stimulation in a dynamic knee extension task.

Estimating spatially distributed turbulent heat fluxes from high-resolution thermal imagery acquired with a UAV system.

Computational modeling of sphingolipid metabolism.

Alteration of heme metabolism in a cellular model of Diamond-Blackfan anemia.

A spatially-averaged mathematical model of kidney branching morphogenesis.