Spatial control of neuronal metabolism through glucosemediated mitochondrial transport regulation
Abstract
Eukaryotic cells modulate their metabolism by organizing metabolic components in response to varying nutrient availability and energy demands. In rat axons, mitochondria respond to glucose levels by halting active transport in high glucose regions. We employ quantitative modeling to explore physical limits on spatial organization of mitochondria and localized metabolic enhancement through regulated stopping of processive motion. We delineate the role of key parameters, including cellular glucose uptake and consumption rates, that are expected to modulate mitochondrial distribution and metabolic response in spatially varying glucose conditions. Our estimates indicate that physiological brain glucose levels fall within the limited range necessary for metabolic enhancement. Hence mitochondrial localization is shown to be a plausible regulatory mechanism for neuronal metabolic flexibility in the presence of spatially heterogeneous glucose, as may occur in long processes of projection neurons. These findings provide a framework for the control of cellular bioenergetics through organelle trafficking.
https://doi.org/10.7554/eLife.40986.001eLife digest
Cells are equipped with power factories called mitochondria that turn nutrients into chemical energy to fuel processes in the cell. Hundreds of mitochondria move throughout the cell, shifting their positions in response to energy demands. This happens via molecular motors that pick the mitochondria up and carry them to new locations. Such movements enable the mitochondria to accumulate in parts of the cell with the greatest energy needs.
Mitochondria of nerve cells or neurons have a particular challenging job, as neurons can be very long and different parts within the cells can have different energy needs. It has been shown that mitochondria stop in regions where nutrients such as sugar are most concentrated. So far, it has been unclear whether this regulated stopping helps control energy balance in neurons.
Here, Agrawal et al. used a computational model of rat neurons to find out whether sugar levels are sufficient in guiding mitochondria. The results showed that the mitochondria only accumulated in highnutrient regions when the sugar concentrations were moderate – not too low and not too high. A specific range of sugar levels was necessary to make this mechanism useful for increasing the efficiency of energy production. Such concentrations match the ones observed in healthy rat brains.
When neurons are unable to meet their energy demands, they stop working and sometimes even die. This is the case in many diseases, including diabetes, dementia, and Alzheimer’s disease. Computer models allow us to explore the complex energy regulation in detail. A better understanding of how neurons regulate their energy production and demand may help us discover how they become faulty in these diseases.
https://doi.org/10.7554/eLife.40986.002Introduction
Cellular metabolism comprises an intricate system of reactions whose finetuned control is critical to cell health and function. A number of quantitative studies have focused on metabolic control through modulating reactant and enzyme concentrations and turnover rates (Grima and Schnell, 2006; Amar et al., 2008). However, these studies generally neglect the spatial organization of metabolic components within the cell. By localizing specific enzymes in regions of high metabolic demand (Laughton et al., 2007; Zecchin et al., 2015), as well as clustering together consecutively acting enzymes (O'Connell et al., 2012), cells have the potential to substantially enhance their metabolism.
Spatial organization is particularly critical in highly extended cells, such as mammalian neurons, whose axons can grow to lengths on the meter scale. Metabolic demand in neurons is spatially and temporally heterogeneous, with especially rapid ATP turnover found in the presynaptic boutons (Rangaraju et al., 2014), and ATP requirements peaking during synaptic activity and neuronal firing (Shulman et al., 2004; Ferreira et al., 2011; Weisová et al., 2009). Neurons rely primarily on glucose as the energy source for meeting these metabolic demands (Peppiatt and Attwell, 2004). Due to the long lengths of neural processes, the glucose supply can vary substantially over different regions of the cell (Ferreira et al., 2011; Weisová et al., 2009; Hall et al., 2012). In myelinated neurons, for instance, it has been speculated that glucose transport into the cell is localized primarily to narrow regions around the nodes of Ranvier (Magnani et al., 1996; Harris and Attwell, 2012; Rosenbluth, 2009), which can be spaced hundreds of microns apart (Ibrahim et al., 1995; Butt et al., 1998). Glucose transporters in neurons have also been shown to dynamically mobilize to active synapses, providing a source of intracellular glucose heterogeneity (Ashrafi et al., 2017). Furthermore, varying levels of activity in the mammalian brain may lead to varying extracellular glucose levels, resulting in spatially heterogeneous nutrient access (Hawkins et al., 1979). Individual axons have been shown to span across multiple regions of the brain (Matsuda et al., 2009), enabling them to encounter regions with different glucose concentrations.
Most ATP production in neurons occurs within mitochondria: motile organelles that range from interconnected networks to individual globular structures that extend throughout the cell. As energy powerhouses and metabolic signaling centers of the cell, mitochondria are critical for neuronal health (Nunnari and Suomalainen, 2012). Their spatial organization within the neuron plays a pivotal role in growth and cell physiology (Li et al., 2004). Defects in mitochondrial transport are involved in the pathologies of several neurological disorders such as peripheral neuropathy and CharcotMarieTooth disease (Baloh, 2008; Baloh et al., 2007).
A number of studies have shown that mitochondria are localized preferentially to regions of high metabolic demand, such as the synaptic terminals (Li et al., 2004; Chang and Reynolds, 2006). Such localization can occur via several molecular mechanisms, mediated by the MiroMilton mitochondrial motor adaptor complex that links mitochondria to the molecular motors responsible for transport (Mishra and Chan, 2016). Increased Ca^{2+} levels at active synapses lead to loading of calcium binding sites on Miro, releasing mitochondria from the microtubule and thereby halting transport (Wang and Schwarz, 2009; Macaskill et al., 2009). High glucose levels can also lead to stalling, through the glycosylation of motor adaptor protein Milton by the glucoseactivated enzyme OGlcNAc transferase (OGT) (Pekkurnaz et al., 2014). This mechanism has been shown to lead to mitochondrial accumulation at glucoserich regions in cultured neurons (Pekkurnaz et al., 2014). It is postulated to regulate mitochondrial spatial distribution, allowing efficient metabolic response to heterogeneous glucose availability.
Mitochondrial positioning relies on an interplay between heterogeneously distributed diffusive signaling molecules (such as Ca^{2+} and glucose), their consumption through metabolic and other pathways, and their effect on motor transport kinetics. While the biochemical mechanisms and physiological consequences of mitochondrial localization have been a topic of much interest in recent years (MacAskill and Kittler, 2010; Mishra and Chan, 2016), no quantitative framework for this phenomenon has yet been developed.
In this work we focus on glucosemediated regulation of mitochondrial transport, developing quantitative models to examine the consequences of this phenomenon for metabolism under spatially varying glucose conditions. Our approach relies on a reactiondiffusion formalism, which describes the behavior of species subject to both consumption and diffusion. Reactiondiffusion systems have been applied to describe the spatial organization of a broad array of cellular processes (Kondo and Miura, 2010), ranging from protein oscillations in E. coli (Howard et al., 2001), to coordination of mitotic signalling (Chang and Ferrell, 2013), to pattern formation in developing embryos (Bunow et al., 1980; Gregor et al., 2005). The response of actively moving particles to spatially heterogeneous, diffusive regulators has also been extensively investigated in the context of chemotaxis (Van Haastert and Devreotes, 2004). In contrast to most chemotactic cells, however, mitochondria have no currently known mechanism for directly sensing glucose gradients. Instead, they are expected to accumulate in response to local glucose concentration only. Our goal is to delineate the regimes in which such a crude form of chemotaxis can lead to substantial spatial organization and enhancement of metabolism.
Specifically, we model the modulation of mitochondrial density with glucose concentration in a tubular axonal region, focusing on two forms of spatial heterogeneity. In one case, we consider an axonal domain between two localized regions of glucose entry, representing the internodal region between nodes of Ranvier in myelinated neurons (Figure 1a). The second case focuses on an unmyelinated cellular region with continuous glucose permeability, embedded in an external glucose gradient (Figure 1b). In both cases, we show that mitochondrial accumulation and enhanced metabolic flux is expected to occur over a limited range of glucose concentrations, which overlaps with physiological brain glucose levels. Our simplified quantitative model allows identification of a handful of key parameters that govern the extent to which glucosemediated mitochondrial halting can modulate metabolism. We establish the region of parameter space where this mechanism has a substantial effect, and highlight its potential importance in neuronal metabolic flexibility and ability to respond to spatially varying glucose.
Results
Minimal model for mitochondrial and glucose dynamics
We begin by formulating a quantitative model to describe the spatial localization of mitochondria that halt in a glucosedependent manner, in the presence of localized sources of glucose. This situation arises in myelinated neurons, which have glucose transporters enriched at the nodes of Ranvier, leading to highly localized sources of glucose spaced hundreds of micrometers apart within the cell (Saab et al., 2013).
Neuronal glucose transporters are known to be bidirectional (Simpson et al., 2007), allowing glucose concentration within the cell to equilibrate with external glucose. For simplicity, we assume rapid transport of glucose through these transporters, so that the internal concentration of glucose at the nodes where transporters are present is assumed to be fixed. The cellular region between two glucose sources is modeled as a onedimensional interval of length L with glucose concentration fixed to a value c_{0} at the interval boundaries (Figure 1a). Glucose diffuses throughout this interval with diffusivity D, while being metabolized by hexokinase enzyme in the first step of mammalian glucose utilization (Figure 1c) (Wilson, 2003).
The concentration of glucose is thus governed by the reactiondiffusion equation,
where $k(x)$ describes the spatial distribution of the hexokinase enzyme as well as the rate of consumption. In the case of spatially uniform, linear consumption [$k(x)=k$, a constant] this equation can be solved directly, yielding a distribution of glucose that falls exponentially from each source boundary, with a decay length $\lambda =\sqrt{D/k}$ (Kholodenko, 2006).
Hexokinase 1 (HK1), the predominant form of hexokinase expressed in neurons, is known to localize preferentially to mitochondria (John et al., 2011), which in mammalian axons can form individual organelles approximately 1 µm in length (Fawcett, 1981). We carry out numerical simulations of Equation 1 where consumption is limited to locations of individual discrete mitochondria, represented by short intervals of length $\mathrm{\Delta}$. Specifically, we define the mitochondria density as $M(x)=n(x)/(\pi {r}^{2}\mathrm{\Delta})$, where $n(x)$ is the number of mitochondria overlapping position $x$, and $r$ is the axon radius. The phosphorylation of glucose by mitochondrial hexokinase is assumed to follow MichaelisMenten kinetics, described by
where ${K}_{M}$ is the saturation constant and k_{g} is the turnover rate of glucose (per unit time per mitochondrion). The turnover rate k_{g} incorporates both the catalytic rate of hexokinase and the number of hexokinase enzymes per mitochondrion. This expression reduces to the case of constant linear consumption when glucose concentration is low ($G\ll {K}_{M}$) and mitochondria are uniformly distributed throughout the region.
In general, glucose consumption depends on the location of mitochondria within the domain. Mitochondrial distribution in neurons is known to be mediated through regulation of their motordriven motility (Chang and Reynolds, 2006; Pekkurnaz et al., 2014). Individual mitochondria switch between processively moving and paused states, modulated by the interplay between kinesin and dynein motors and the adaptor proteins that link these motors to the mitochondria (Schwarz, 2013). In our model, we simulate mitochondria as stochastically switching between a processive walking state that moves in either direction with velocity $v$ and a stationary state. The rate of initiating a walk (k_{w}) is assumed to be constant, while the halting rate (${k}_{s}(x)$) can be spatially heterogeneous. For simplicity, we assume the mitochondria are equally likely to move in the positive (+) or negative () direction each time they initiate a processive walk (Figure 1b).
It has recently been demonstrated that the key motor adaptor protein (Milton) is sensitive to glucose levels, halting mitochondrial motility when it is modified through OGlcNAcylation by the OGT enzyme (Pekkurnaz et al., 2014). Our model employs a highly simplified description of mitochondrial dynamics, which assumes that all pauses are associated with such an OGlcNAcylation event. Recovery from the pause at the constant rate k_{w} corresponds to removal of the modification through the activity of the complementary enzyme OGlcNAcase (OGA). Although there is evidence indicating longterm glucose deprivation can reduce OGA expression (Zou et al., 2012), for simplicity we assume in our model that OGA activity is independent of glucose levels. In vivo axonal mitochondria have been observed to undergo shortlived sporadic pausing while continuing to move processively in their previous anterograde or retrograde direction (Russo et al., 2009; Wang and Schwarz, 2009). Such pauses are subsumed into an effective processive velocity $v$ in our model. Other sources of pausing, such as Ca^{2+}regulated motor disengagement, PINK1/Parkinmediated detachment of motors, and anchoring to the microtubules by syntaphilin (Schwarz, 2013), are not considered here in order to focus specifically on the effect of glucosedependent mitochondrial spatial organization.
Upon entry into the cell, the first ratelimiting step of glucose metabolism is its conversion into glucose6phosphate by hexokinase. Further downstream metabolic pathways split, with much of the flux going to glycolysis while a small fraction is funneled into the pentose phosphate pathway and the hexosamine biosynthetic pathway (HBP). The HBP produces UDPGlcNAc, the sugar substrate for OGlcNAcylation (Figure 1c) (Hart et al., 2011). In our model, we assume that the rate of UDPGlcNAc production equals the rate of glucose conversion by hexokinase, scaled by the fraction of G6P that is channeled into the hexosamine biosynthetic pathway. This assumption is valid if, at each point of pathway branching, the MichaelisMenten saturation constants for the two branches are similar. This in fact appears to be the case for both the branching of the pentose phosphate pathway and glycolysis from the hexosamine biosynthetic pathway which is the focus of this work (see Appendix 2). Consequently, saturation of the initial glucose conversion step will imply saturation of the entire hexosamine biosynthetic pathway. We therefore model the kinetics of Milton modification using the same MichaelisMenten form as for hexokinase activity, with the pathway flux leading to Milton modification subsumed within a rate constant for mitochondrial stopping (k_{s}).
We note that the subcellular organization of the intermediates in the conversion from glucose into OGlcNAcylated Milton is largely unknown. In our model, we make the extreme case assumption that all intermediates are localized to mitochondria, with only the initial glucose substrate capable of diffusing through the cytoplasm. We note that cytoplasmic diffusion of any of the pathway intermediates would attenuate the effect on mitochondrial localization. Our simplified model thus gives an upper limit on the extent to which mitochondria can localize at high glucose regions through the Milton modification mechanism. Following these simplified assumptions, we treat the kinetics of mitochondrial halting as dependent only on the local glucose concentration, according to the functional form
where ${K}_{M}$ is the MichaelisMenten constant of hexokinase.
We proceed to evolve the simulation forward in time, with glucose consumption localized to regions within $\pm \mathrm{\Delta}/2$ of each discrete mitochondrial position (details in Materials and methods). A snapshot of one simulation run is shown in Figure 2a, highlighting the accumulation of stationary mitochondria in the high glucose regions near the ends of the domain.
We are interested primarily in investigating the steadystate distribution of mitochondria and glucose in this system, averaged over all possible mitochondrial trajectories. We thus proceed to coarsegrain our model by treating the distribution of mitochondria as a continuous field $M(x)={W}_{+}(x)+{W}_{}(x)+S(x)$, where ${W}_{+}(x)$ is the distribution of mitochondria walking in the positive direction, ${W}_{}(x)$ is the distribution of those walking in the negative direction, and $S(x)$ is the distribution of stationary mitochondria. We can then write down the coupled differential equations governing the behavior of the mitochondrial distributions as:
The glucose distribution evolves according to Equation 1 with consumption rate $k(x)$ given by Equation 2. The boundary conditions at the ends of the domain are assumed to be reflective for the mitochondrial distributions, and to have a fixed glucose concentration ${c}_{0}$. The stationary state for this system can be calculated numerically (see Materials and methods). The formulation with a continuous mitochondrial density faithfully represents the behavior of simulations with discrete mitochondria, as illustrated in Figure 2b.
The steadystate spatial distribution of mitochondria and glucose in the continuous system depend on six parameters: ${k}_{s}/{k}_{w},{K}_{M},{c}_{0},D,L,{k}_{g}\overline{M}$ where $\overline{M}$ is the average mitochondrial density in the axon (number of mitochondria per unit volume) . Estimates of physiologically relevant values are provided in Table 1. Dimensional analysis indicates that three of these parameters can be used to define units of time, length, and glucose concentration, leaving three dimensionless groups. We choose to use the following three dimensionless parameters, each of which has an intuitive physical meaning:
Here $\widehat{\lambda}$ is the lengthscale of glucose decay relative to the domain length, ${\widehat{c}}_{0}$ is the boundary glucose concentration relative to the saturation constant ${K}_{M}$, and ${\widehat{k}}_{s}$ is the ratio of stopped to walking mitochondria at high glucose levels. We proceed to explore the steadystate distribution of mitochondria and glucose as a function of these three parameters.
Mitochondrial localization requires limited range of external glucose
In order for mitochondria to preferentially accumulate at the source of glucose via a glucosedependent stopping mechanism, three criteria must be met. First, the glucose concentration needs to be higher at the source than in the bulk of the cell, as occurs when the decay length due to consumption is much smaller than the size of the domain ($\widehat{\lambda}\ll 1$). Second, if glucose levels become too high (${\widehat{c}}_{0}\gg 1$) then both glucose consumption rates and stopping rates of the mitochondria become saturated, leading to a flattening of glucose and mitochondrial distributions (Figure 3). There is thus an upper limit on the possible external glucose concentrations that will yield mitochondrial localization at the edges of the domain. Finally, the mitochondria must spend a substantial amount of time in the stationary state, since walking mitochondria will be broadly distributed throughout the domain. Because the stopping rate is itself dependent on the glucose concentration, this criterion implies that very low concentrations will also not allow mitochondrial localization. Figure 3 shows the distribution of glucose and mitochondria at different values of the external glucose ${\widehat{c}}_{0}$, illustrating that accumulation of mitochondria at the edges requires intermediate glucose levels.

Figure 3—source data 1
 https://doi.org/10.7554/eLife.40986.007
To characterize the distribution of mitochondria along the interval, we introduce an accumulation metric $A$, defined by
where ${\sigma}^{2}$ is the variance in the mitochondrial distribution. This metric scales from $A=0$ for a uniform distribution to $A=1$ for two narrow peaks at the domain edges. Mitochondrial distributions with several different values of the accumulation metric are shown in Figure 3a. We use a cutoff of $A=0.2$ to define distributions where the mitochondria are localized at the glucose source.
We explore the dependence of the mitochondrial accumulation on the three dimensionless parameters defining the behavior of the system: the stopping rate constant ${\widehat{k}}_{s}$, the glucose decay length $\widehat{\lambda}$, and the external concentration ${\widehat{c}}_{0}$. Because only the stopped mitochondria localize near the glucose sources, increasing the fraction of mitochondria in the stopped state (increased ${\widehat{k}}_{s}$) inevitably raises the overall accumulation (Figure 4a). The fraction of mitochondria in the stopped state will depend on both ${\widehat{k}}_{s}$ and the overall levels of glucose, as dictated by ${\widehat{c}}_{0}$ (Figure 4b). Experimental measurements indicate that at high glucose concentrations, approximately 95% of mitochondria are in the stationary state (Pekkurnaz et al., 2014). We are thus interested primarily in the parameter regime of high stopping rates: ${\hat{k}}_{s}\gtrsim 10$. The limited range of concentrations that lead to mitochondrial accumulation at the edges of the domain can be seen in Figure 4a.

Figure 4—source data 1
 https://doi.org/10.7554/eLife.40986.009

Figure 4—source data 2
 https://doi.org/10.7554/eLife.40986.010

Figure 4—source data 3
 https://doi.org/10.7554/eLife.40986.011
For a high stopping rate (${\widehat{k}}_{s}=10$), we then calculate the mitochondrial accumulation as a function of the remaining two parameters: $\widehat{\lambda},{\widehat{c}}_{0}$. Here, again, we note that only intermediate glucose concentrations result in accumulation, with the range of concentrations becoming narrower as the decay length $\widehat{\lambda}$ becomes comparable to the domain size (Figure 4c). We can establish the concentration range within which substantial accumulation is expected, by setting a cutoff $A=0.2$ on the accumulation metric and calculating the resulting phase diagram (Figure 4d). Below the lower concentration cutoff, insufficient mitochondria are in the stationary state and so no localization is seen. This lower cutoff disappears in the limit of infinite ${\widehat{k}}_{s}$. At intermediate concentrations, mitochondria are localized near the domain edges. Above the upper concentration cutoff, no localization is observed due to saturation of the MichaelisMenten kinetics.
Using empirically derived approximations for the rate of glucose consumption by mitochondria and the diffusivity of glucose in cytoplasm (see Table 1), we estimate the decay length parameter as $\widehat{\lambda}\approx 0.03$. The mitochondria are then expected to localize near the glucose source only if ${\hat{c}}_{0}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}66$. Because the saturation concentration for hexokinase is quite low (${K}_{M}\approx 0.03$mM) (Wilson, 2003), we would expect mitochondrial accumulation for glucose concentrations below about $2$ mM. We note that physiological brain glucose levels have been measured at 0.7 − 1.3 mM, depending on the brain region (McNay et al., 2001), implying that glucosedependent halting of mitochondrial transport would be expected to result in localization of mitochondria at nodes of Ranvier.
Glucosedependent halting can increase metabolic flux under physiological conditions
Localizing mitochondria to the glucose entry points is expected to increase the flux of glucose entering the cell, thereby potentially enhancing the overall metabolic rate. We calculate the overall effect of transportbased regulation on the net metabolic flux within the simplified model with localized glucose entry. Figure 5 shows the effect of increasing mitochondrial stopping rates (${\widehat{k}}_{s}$) on the total rate of glucose consumption in the interval between nodes of glucose influx. At low ${\widehat{k}}_{s}$ values, mitochondria are distributed uniformly throughout the interval. At high ${\widehat{k}}_{s}$ values and at sufficiently low glucose concentrations, the mitochondria cluster in the regions of glucose entry, increasing the overall consumption rate by up to 40% at physiologically relevant glucose levels (c_{0} = 1 mM). We note that in hypoglycemic conditions, glucose levels can drop to 0.1 mM (Silver and Erecińska, 1994), further increasing the magnitude of this effect.

Figure 5—source data 1
 https://doi.org/10.7554/eLife.40986.013
In the case of limited glucose transport into the cell, intracellular glucose levels could be significantly below the concentrations outside the cell. Measurements of intracellular glucose in a variety of cultured mammalian cell types indicate internal concentrations within the range of 0.07 − 1mM, up to an order of magnitude lower than glucose concentrations in the medium (John et al., 2008). However, neuronal cells are known to express a particularly efficient glucose transporter (GLUT3) (Simpson et al., 2008), and these transporters have been shown to be highly concentrated near the nodes of Ranvier (Magnani et al., 1996; Rosenbluth, 2009). We therefore assume that glucose import into the nodes is not rate limiting for myelinated neurons in physiological conditions. Introducing a finite rate of glucose transport would effectively decrease the intracellular glucose concentration at the nodes ${c}_{0}$, increasing the enhancement in metabolic flux due to mitochondrial localization. In subsequent sections, we explore the role of limited glucose import in unmyelinated axons with spatially uniform glucose permeability.
Model for spatial organization in a glucose gradient
Extracellular brain glucose levels exhibit substantial regional variation, particularly under hypoglycemic conditions where more than tenfold differences in local glucose concentrations have been reported (Paschen et al., 1986). Because individual neurons can traverse multiple different brain regions (Matsuda et al., 2009), a single axon can be subjected to heterogeneous glucose levels along its length. This raises the possibility that glucosedependent mitochondrial localization can play a role in neuronal metabolic flexibility even in the case where glucose entry into the cell is not localized to distinct nodes. We thus extend our model to quantify the distribution of mitochondria in an axon with limited but spatially uniform glucose permeability that is subjected to a gradient of external glucose. This situation is relevant, for instance, to unmyelinated neurons in infant brains, as well as to in vitro experiments with neurons cultured in a glucose gradient (Pekkurnaz et al., 2014).
In this model, the extracellular environment provides a continuous source of glucose whose influx is limited by the permeability of the cell membrane. Intracellular glucose dynamics are then defined by the reactiondiffusion equation
where the first term corresponds to diffusive glucose spread, the second to a spatially varying metabolism of glucose, and the third to the entry of glucose into the cell. Here, $G}_{\mathrm{e}\mathrm{x}\mathrm{t}$ is the external glucose concentration, and $P(x)$ is the membrane permeability to glucose, which we assume to depend in a MichaelisMenten fashion on the difference between external and internal glucose concentration:
where $P$ is the spatially uniform permeability constant in units of length per time. This functional form incorporates two known features of glucose transporters: (1) they are bidirectional, so that the overall flux through the transporter at low glucose levels should scale linearly with the difference between external and internal glucose (Carruthers, 1990); (2) neuronal glucose transporters saturate at high glucose levels (GLUT3 ${K}_{M\phantom{\rule{thinmathspace}{0ex}}P}\approx 3$mM (Maher et al., 1996), with an even higher saturation constant for GLUT4 (Nishimura et al., 1993). When the difference in glucose levels is low, the overall flux of glucose entering the cell reduces to $P({G}_{\mathrm{e}\mathrm{x}\mathrm{t}}(x)G(x))$. Mitochondria dynamics are defined as before (Equation 4), and we again assume MichaelisMenten kinetics for glucose metabolism by hexokinase localized to mitochondria (Equation 2).
We note that the dynamics in Equation 6 are governed by three timescales: the rate of glucose transport down the length of the axon, rate of glucose consumption, and rate of glucose entry. The first of these rates becomes negligibly small in the limit $L\gg \sqrt{D(G+{K}_{M})/({k}_{g}\overline{M})}$. Because internal glucose levels can never exceed the external concentrations, in the range where ${G}_{\mathrm{e}\mathrm{x}\mathrm{t}}<10$mM, the rate of diffusive transport should become negligible for $L\gg 150\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. In the limit where intracellular glucose is much less than K_{M}, this criterion reduces to $\widehat{\lambda}\ll 1$, indicating that glucose diffuses over a very small fraction of the interval before being consumed. The interval length $L$ in this model represents an axonal length which can range over many orders of magnitude. We focus on axon lengths above several hundred microns, allowing us to neglect the diffusive transport of intracellular glucose (see Appendix 3).
The steadystate glucose profile can then be determined entirely by the local concentration of mitochondria and external glucose. For a given mitochondrial density $M(x)$ and external glucose profile ${G}_{\mathrm{e}\mathrm{x}\mathrm{t}}(x)$, the corresponding intracellular glucose concentration can be found directly by solving the quadratic steadystate version of Equation 6 without the diffusive term. However, the steadystate mitochondrial distribution cannot be solved locally, because the limited number of mitochondria within the axon couples the mitochondrial density at different positions. We thus employ an iterative approach to numerically compute the steadystate solution for both glucose and mitochondrial density under a linear external glucose gradient $G}_{\mathrm{e}\mathrm{x}\mathrm{t}}={G}_{\mathrm{m}\mathrm{i}\mathrm{n}}+({G}_{\mathrm{m}\mathrm{a}\mathrm{x}}{G}_{\mathrm{m}\mathrm{i}\mathrm{n}})\frac{x}{L$ (see Materials and methods).
For parameter combinations where intracellular glucose concentrations are above K_{M} but well below $G}_{\mathrm{e}\mathrm{x}\mathrm{t}$, the entry and consumption processes for glucose are both saturated. There is then a steep transition between two different regimes. In one regime, glucose entry exceeds consumption and internal glucose levels approach the external concentrations. In the other, consumption dominates and glucose levels drop below saturating concentrations. The key dimensionless parameter governing this transition can be defined as the ratio of entry to consumption rates:
This ratio can be modulated in the cell either by recruiting varying amounts of glucose transporters (adjusting $P$) or changing the total amount of active hexokinase (adjusting ${k}_{g}\overline{M}$).
The remaining dimensionless parameters determining the behavior of this simplified model are the external glucose concentration relative to the hexokinase saturation constant ($\hat{G}}_{\mathrm{e}\mathrm{x}\mathrm{t}}={\overline{G}}_{\mathrm{e}\mathrm{x}\mathrm{t}}/{K}_{M$), the relative magnitude of the glucose gradient, $\mathrm{\Delta}{\hat{G}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=({G}_{\mathrm{m}\mathrm{a}\mathrm{x}}{G}_{\mathrm{m}\mathrm{i}\mathrm{n}})/{\overline{G}}_{\mathrm{e}\mathrm{x}\mathrm{t}}$, the ratio of stopped to walking mitochondria $\widehat{{k}_{s}}={k}_{s}/{k}_{w}$, and the saturation constant for glucose transporters ${K}_{M\phantom{\rule{thinmathspace}{0ex}}P}/{K}_{M}\approx 96$. The last parameter is expected to remain approximately constant in neuronal cells. The average external glucose concentration and glucose gradient are expected to vary substantially depending on the glycemic environment to which the neuron is exposed. We note that $\mathrm{\Delta}{\hat{G}}_{\mathrm{e}\mathrm{x}\mathrm{t}}$ has a maximum possible value since the minimal glucose concentration cannot drop below $0$zero. We proceed to analyze the limiting case where the glucose gradient is as steep as possible for any given value of average external glucose ($\mathrm{\Delta}{\hat{G}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=2$).
Mitochondrial arrest enables metabolic enhancement under glucose gradient
We quantify the amount of mitochondrial accumulation at the high glucose side of the domain by calculating the total mitochondrial density within the distal 10% of the interval compared to a uniform distribution, in analogy to experimental measurements (Pekkurnaz et al., 2014). Substantial enrichment in the high glucose region occurs when glucose entry into the cell cannot keep up with consumption ($\gamma \ll 1$) and the intracellular glucose levels drop below the hexokinase saturation concentration K_{M}, as can be seen in the glucose and mitochondrial distributions plotted in Figure 6a–c. The interplay between external glucose levels and the entry/consumption rates is illustrated in Figure 6d. For external glucose concentrations well above K_{M} there is a sharp transition to mitochondrial enrichment at $\gamma \phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1$. At the lowest levels of intracellular glucose, accumulation is again reduced because a very small fraction of mitochondria are found in the stopped state. In the limit of high k_{s}, mitochondrial accumulation would occur for arbitrarily low values of $\gamma $ (Figure 6—figure supplement 1). We note that because glucose entry and turnover are much faster than diffusive spread for biologically relevant parameter regimes, the model results do not depend on the cell length $L$ (Appendix 3).

Figure 6—source data 1
 https://doi.org/10.7554/eLife.40986.016

Figure 6—source data 2
 https://doi.org/10.7554/eLife.40986.017

Figure 6—source data 3
 https://doi.org/10.7554/eLife.40986.018
Experimental measurements of mitochondrial enrichment in cultured neurons subjected to a gradient of $0$ to $5$mM glucose have indicated an approximately 20% enrichment in mitchondrial counts at the axonal region exposed to high glucose. We note that using published estimates of typical glucose permeability and mitochondrial glucose turnover for mammalian cells (Table 1) yields a ratio of entrance and consumption rates of $\gamma \approx 1.9$ for this experimental system. Because this ratio is above 1, we would not expect to see substantial mitochondrial enrichment. To result in the experimentally observed enrichment at high glucose, the ratio $\gamma $ would need to be reduced by approximately a factor of 2, implying the existence of additional regulatory mechanisms. Modulation of $\gamma $ could be achieved by either decreasing the number of glucose transporters in the cell (reducing $P$) or upregulating total hexokinase levels (increasing k_{g}). Neurons are believed to regulate both the density of glucose transporters and hexokinase activity in response to external glucose concentrations and varying metabolic demand (Fujii and Beutler, 1985; Robey et al., 1999; Duelli and Kuschinsky, 2001). In particular, adaptation to glycemic levels well above physiological values, as well as possibly reduced synaptic activity in a cultured environment, may result in downregulation of glucose transporters, lowering the value of $\gamma $. The discrepancy between model prediction and observed mitochondrial accumulation highlights the existence of additional regulatory pathways not included in the current model whose role could be explored in further studies that directly quantify glucose entry and consumption rates in cultured neurons.
Physiological brain glucose levels have been measured at 0.7 mM  1.3 mM (McNay et al., 2001), with hypoglycemic levels dipping as low as 0.1 mM and hyperglycemic levels rising up to $4$mM (Silver and Erecińska, 1994). Axons that stretch across different brain regions with varying glucose levels can thus be subject to a glucose gradient with $\overline{G}}_{\mathrm{e}\mathrm{x}\mathrm{t}$ on the order of 1 mM (white line on Figure 6d). We note that the physiological range overlaps substantially with the region of high mitochondrial accumulation, indicating that glucosedependent halting can modulate mitochondrial distribution under physiologically relevant glycemic levels.
By accumulating mitochondria at the cellular region subjected to higher external glucose, the metabolic flux in that region can be substantially enhanced. In (Figure 6e) we plot the enhancement in glucose consumption rates (compared to the case with uniformly distributed mitochondria) within the $10\%$ of cellular length subjected to the highest glucose concentrations. Metabolic enhancement occurs within a narrow band of the $\gamma $ parameter. The dropoff in enhancement at low values of the internal glucose concentration (low $\gamma $) is due to the coupling between glucose levels and mitochondrial localization. Specifically, mitochondrial accumulation at the region subject to high glucose concentration increases the local rate of consumption in that region, driving down local internal glucose levels. Consequently, the difference in internal glucose concentrations between the two ends of the cell is decreased when internal levels fall substantially below _{M} (Figure 6b), reducing the enhancement of metabolic flux. Although mitochondrial accumulation decreases metabolic flux in the low glucose region, the total rate of glucose consumption integrated throughout the cell is enhanced by up to approximately 14% when $\gamma \approx 1$ (Figure 6f).
It is interesting to note that the typical physiological range of external glucose levels spans the narrow band of parameter space where metabolic enhancement is expected (white lines on Figure 6e,f). These results implicate glucosedependent mitochondrial stopping as a quantitatively plausible mechanism of metabolic flexibility, increasing metabolism in regions with high nutrient availability for axonal projections that span between hypoglycemic and euglycemic regions. The magnitude of this effect can be tightly controlled by the cell through modulating overall rates of glucose entry and consumption. Thus, by coupling mitochondrial transport to local glucose levels, wholecell changes in hexokinase or glucose transporter recruitment can be harnessed to tune the cell’s response to spatially heterogeneous glucose concentrations.
Discussion
The minimal model described here provides a quantitative framework to explore the interdependence of glucose levels and mitochondrial motility and their combined effect on neuronal metabolic flux. Glucosemediated halting of mitochondrial transport is shown to be a plausible regulatory mechanism for enhancing metabolism in cases with spatially heterogeneous glucose availability in the neuron.
We have quantitatively delineated the regions in parameter space where such a mechanism can have a substantial effect on mitochondrial localization and metabolic flux. Specifically, mitochondrial positioning requires both sufficient spatial variation in intracellular glucose and sufficiently low absolute glucose levels compared to the saturation constant of the hexokinase enzyme. In the case of tightly localized glucose entry (as at the nodes of Ranvier), intracellular spatial heterogeneity requires a small value of the dimensionless length scale for glucose decay ($\widehat{\lambda}=\sqrt{D{K}_{M}/{k}_{g}\overline{M}{L}^{2}}\ll 1$). For physiologically estimated values, mitochondrial localization to the nodes is expected to occur for glucose levels below approximately 2 mM, comparable to physiological brain glucose concentrations (McNay et al., 2001; John et al., 2008). In the case where glucose can enter homogeneously throughout the cell surface (as with unmyelinated axons), heterogeneity can arise from an external glucose gradient. We show that metabolic enhancement through mitochondrial positioning occurs in a narrow range of the key parameter $\gamma =(2P{K}_{M\phantom{\rule{thinmathspace}{0ex}}P}{\overline{G}}_{\mathrm{e}\mathrm{x}\mathrm{t}})/({k}_{g}\overline{M}({K}_{M\phantom{\rule{thinmathspace}{0ex}}P}+{\overline{G}}_{\mathrm{e}\mathrm{x}\mathrm{t}}))$, which describes the ratio of glucose entry to glucose metabolism, and that this narrow range intersects with physiological estimates.
The model developed here is intentionally highly simplified, encompassing a minimal set of parameters necessary to describe glucosedependent mitochondrial localization. Other regulatory pathways that determine mitochondrial positioning are not included in this basal model. In particular, we do not consider here calciumbased transport regulation, which is known to localize mitochondria to regions of synaptic activity (Zhang et al., 2010; Wang and Schwarz, 2009; MacAskill and Kittler, 2010; Macaskill et al., 2009). Upregulating OGT signaling in cultured cells has been shown to decrease the fraction of motile mitochondria by a factor of three, while reducing endogenous OGT nearly doubles the motile fraction, indicating that a substantial number of stationary mitochondria are stopped as a result of OGT activity (Pekkurnaz et al., 2014). Our model assumes the extreme case where all stopping events are triggered in a glucosedependent manner, thereby isolating the effect of glucose heterogeneity. Stopping mechanisms dependent on neuronal firing activity could alter mitochondrial distribution in concert with glucosedependent halting, increasing the density of mitochondria at presynaptic boutons or near areas of localized calcium influx as at the nodes of Ranvier (Zhang et al., 2010). We note that mitochondria have previously been shown to accumulate at spinal nodes of Ranvier in response to neuronal firing activity (Fabricius et al., 1993; Zhang et al., 2010). The mechanism described here provides an additional driving force for mitochondrial localization near the nodes even in quiescent neurons.
Additional metabolic feedback loops, not included in our model, may result in a more complex dependence of mitochondrial stopping on glucose concentration. In particular, both the pentose phosphate pathway and glycolysis generate intermediates that feed back into UDPGlcNAc production by the hexosamine biosynthetic pathway (Kruger and von Schaewen, 2003; Shirato et al., 2011). Furthermore, several of the enzymes involved in the metabolic pathways linking glucose levels to Milton OGlcNacylation may be regulated in a glucosedependent manner. For example, the activity of the fructose6phosphate metabolizing enzyme GFAT is believed to be regulated by intermediates in the hexosamine pathway (Traxinger and Marshall, 1991) and OGlcNAc transferase (OGT) itself is directly regulated by UDPGlcNAc levels (Hart et al., 2007). Other enzymes, such as the deGlcNAcylating enzyme OGA exhibit long term regulation of expression in response to altered glucose levels (Zou et al., 2012). These regulatory mechanisms provide additional potential routes of metabolic control through mitochondrial positioning.
Several key parameters that regulate mitochondrial localization in response to glucose heterogeneity can be dynamically regulated in neurons. Specifically, the rate of glucose consumption (${k}_{g}\overline{M}$) can be tuned by modulating the concentration or activity of hexokinase within mitochondria or by altering total mitochondrial size and number. This parameter controls both the glucose decay length $\widehat{\lambda}$ in the case of localized glucose influx and the ratio of glucose entry to consumption $\gamma $ in the case of spatially distributed entry. We note that our model assumes hexokinase to be localized exclusively to mitochondria. The predominant form of hexokinase in the brain (HK1) is known to bind reversibly to the mitochondrial membrane, with exchange between a mitochondriabound and a cytoplasmic state believed to contribute to the regulation of its activity (Golestani et al., 2007). Release of hexokinase into the cytoplasm would result in more spatially uniform glucose consumption, negating the metabolic enhancement achieved through mitochondrial localization.
An additional parameter known to be under regulatory control is the rate of glucose entry into the neuron ($P$). The glucose transporters GLUT3 (Simpson et al., 2008; Duelli and Kuschinsky, 2001; Weisová et al., 2009) and GLUT4 (Ashrafi et al., 2017) have been shown to be recruited to the plasma membrane in response to neuronal firing activity. Interestingly, transporter densities are themselves spatially heterogeneous, concentrating near regions of synaptic activity (Ashrafi et al., 2017; Ashrafi and Ryan, 2017). The model described in this work quantifies the extent to which a locally increased glucose influx can enhance total metabolic flux, given the ability of mitochondria to accumulate at regions of high intracellular glucose.
A number of possible feedback pathways linking glucose distribution and mitochondrial positioning are not included in our basic model. For instance, hexokinase release from mitochondria into the cytoplasm (potentially altering k_{g}) is known to be triggered at least in part by glucose6phosphate, the first byproduct in glucose metabolism (Crane and Sols, 1954). Chronic hypoglycemia has been linked to an upregulation in GLUT3 in rat neurons (Uehara et al., 1997), which would in turn lead to an increased glucose uptake ($P$). The fraction of glucose funneled into the hexosamine biosynthetic pathway (incorporated within k_{s}) can also be modified through feedback inhibition of GFAT by the downstream metabolic product UDPGlcNAc (Li et al., 2007). Such feedback loops imply that several of our model parameters ($P$, k_{g}, k_{s}) are themselves glucosedependent and may become spatially nonuniform in response to heterogeneous glucose. Incorporating these effects into a spatially resolved model of metabolism would require quantifying the dynamics of both the feedback pathways and mitochondrial positioning, and forms a promising avenue for future study.
Control of glucose entry and consumption underlies cellular metabolic flexiblity, and defects in the associated regulatory pathways can have grave consequences for neuronal health. Misregulation of hexokinase has been highlighted as a contributor to several neurological disorders, ranging from depression (Regenold et al., 2012) to schizophrenia (Shan et al., 2014). Neuronal glucose transporter deficiency has been linked to autism spectrum disorders (Zhao et al., 2010) and Alzheimer’s disease (Liu et al., 2008). Furthermore, defects in mitochondrial transport, with the consequent depletion of mitochondria in distal axonal regions, contribute to peripheral neuropathy disorders (Baloh, 2008).
Glucosedependent mitochondrial localization provides an additional layer of control, beyond conventionally studied regulatory mechanisms, which allows the cell to respond to spatial heterogeneity in glucose concentration. Our analysis paves the way for quantitative understanding of how flexible regulation of metabolism can be achieved by controlling the spatial distribution of glucose entry and consumption.
Materials and methods
Source code (in MATLAB [The MathWorks Inc, 2015]) for all simulations and numerical calculations is available at: https://github.com/lenafabr/mitoManuscriptCodes (copy archived at https://github.com/elifesciencespublications/mitoManuscriptCodes).
Discrete mitochondria simulations
Request a detailed protocolWe simulate the internodal space of the axon, between localized nodes of glucose entry, as a onedimensional domain for a reaction diffusion system with motile reaction sinks. The glucose concentration field is discretized over 100 equidistant points along the domain. Its dynamics are governed by the reaction diffusion equation (Equation 1), evolved forward over timesteps of $\delta t$ using the forward Euler method. Because forward Euler methods have stringent conditions for stability and convergence, we use a timestep that is much smaller than both the glucose decay timescale and the timescale associated with diffusion over our spatially discretized grid (see below).
The number of mitochondria in the domain is calculated according to $N=\overline{M}L\pi {r}^{2}\approx 38$, where the mitochondrial density $\overline{M}$, internodal distance $L$, and axonal radius $r$ are estimated from published data (Table 1; Appendix 1). The mitochondria are treated as discrete intervals of length $\mathrm{\Delta}$ = 1 μm, with the position of each mitochondrial center updated at each timestep. Over each time step, every motile mitochondrion moves a distance of $\pm v\delta t$, (with transport velocity $v$ = 1 μm/s) and switches to a stationary state with probability $1\mathrm{exp}({k}_{s}\delta t)$, where ${k}_{s}(x)$ is a function of the center position of that mitochondrion (Equation 3). Mitochondria that reach within a distance of $\mathrm{\Delta}/2$ from the ends of the domain are reflected, reversing their velocity while remaining motile. Analogously, every stationary mitochondrion switches to a motile state on each timestep with probability $1\mathrm{exp}({k}_{w}\delta t)$. Processive walks are initiated with equal probability in either direction.
At any given time, the spatial density of mitochondria is calculated from the location of mitochondrial centers at positions ${x}_{1},\mathrm{\dots}{x}_{N}$, according to $M(x)=n(x)/(\pi {r}^{2}\mathrm{\Delta})$, where
$n(x)=\sum _{i=1}^{N}\left[\theta (x{x}_{i}+\mathrm{\Delta}/2)\theta (x{x}_{i}\mathrm{\Delta}/2)\right]$,
is the number of mitochondria overlapping spatial position $x$ and $\theta $ is the Heaviside step function.
We integrate the simulation forward in timesteps of $\delta t=0.2\frac{\mathrm{\Delta}{x}^{2}}{D}$, where $\mathrm{\Delta}x$ is the spatial discretization. This timescale is much smaller than the relevant decay time for glucose consumption $\left[{\tau}_{g}={\left(\frac{{k}_{g}\overline{M}}{{K}_{M}}\right)}^{1}\right]$. Using these small timesteps allows for stability and robust convergence with the forward Euler method. The simulation proceeds for 10^{7} steps. Simulations are repeated 100 times to obtain the histogram shown in Figure 2. Convergence to steadystate is established by comparing to calculations with the continuum model described in the subsequent sections.
Mitochondrial distribution for spatially varying stopping rate
Request a detailed protocolFor an arbitrary spatial distribution of stopping rates ${k}_{s}(x)$ the corresponding steadystate mitochondrial distribution can be calculated directly by solving the equations for mitochondrial transport (Equation 4):
Because our model assumes symmetry between anterograde and retrograde mitochondrial transport, as well as equal glucose concentrations at either boundary of the domain, we take ${W}_{}={W}_{+}$, implying that the population of walking mitochondria must be spatially constant. Consequently, the population of stopped mitochondria is proportional to the stopping rate ($S=C{k}_{s}(x)/{k}_{w}$). The constant $C$ can be calculated from the normalization condition,
The overall steadystate distribution of mitochondria is then given by,
Because the stopping rate is an explicit function of glucose concentrations $\left[{k}_{s}(x)=\frac{{k}_{s}G(x)}{{K}_{M}+G(x)}\right]$, this approach allows us to find the steadystate mitochondrial distribution for any fixed distribution of glucose.
Numerical solution for steadystate distributions with localized glucose entry
Request a detailed protocolWe solve for steadystate glucose and mitochondrial distributions using a numerical method that evolves the glucose concentration forward in time while explicitly setting the mitochondrial concentration to its steadystate value at each step. The glucose distribution is initialized according to the steadystate solution for uniform consumption (Equation 13). Mitochondrial density $M(x)$ is calculated from the glucose distribution according to Equation 11 and Equation 3.
The glucose distribution $G(x)$, in turn, evolves according to the mitochondrial distribution as given by Equation 1 and Equation 2. The glucose profile is integrated forward with a timestep $\delta t={10}^{5}{L}^{2}/D$. The distributions are assumed to be converged once the root mean squared rate of glucose change drops below the minimal cutoff: ${10}^{6}{k}_{g}\overline{M}$. Results of the continuous mitochondrial distribution model are shown to match the discrete mitochondria simulations (Figure 2b). All subsequent analysis is done in the continuum limit.
Analytical solution for low glucose limit
Request a detailed protocolWe validate our numerical calculations by comparing to the analytically tractable solution in the limit of low glucose and nearly uniform mitochondrial distribution. In the limit of spatially uniform, linear consumption, the steadystate reactiondiffusion equation for glucose can be expressed as
where $k={k}_{g}\overline{M}/{K}_{M}$ is the constant consumption rate.
Assuming fixed glucose concentrations (c_{0}) at the boundaries of the domain, the steadystate glucose distribution is then given by
with $\lambda =\sqrt{\frac{D}{k}}$ defining the glucose decay lengthscale. This quantity is a measure of how far glucose diffusively penetrates into the domain before being consumed by hexokinase. It is scaled by the size of the domain to give the dimensionless decay length scale $\widehat{\lambda}=\sqrt{\frac{D{K}_{M}}{{k}_{g}\overline{M}{L}^{2}}}$ used as a key parameter in our model with localized glucose entry:
Steadystate distribution with uniform permeability in the slow diffusion limit
Request a detailed protocolFor the model with spatially uniform glucose permeability, we solve directly for the steady state distributions of glucose and mitochondria in the limit of slow diffusivity. When diffusion along the domain is slow compared to the timescales of glucose consumption and glucose import, the steadystate equation for glucose concentration is given by a simplified form of Equation 6:
Substituting $k(x)=\frac{{k}_{g}M(x)G(x)}{G(x)+{K}_{M}}$ and $P(x)=\frac{(2/r)P\phantom{\rule{thinmathspace}{0ex}}{K}_{M\phantom{\rule{thinmathspace}{0ex}}P}}{{K}_{M\phantom{\rule{thinmathspace}{0ex}}P}+{G}_{ext}(x)G(x)}$, we get a quadratic equation in $G(x)$;
For a given mitochondrial profile, this quadratic equation is solved to find $G(x)=G(M(x))$. The mitochondrial distribution, $M(x)$ is then updated according to Equation 11 and Equation 3. We thus arrive at an iterative solution for $G(x)$ and $M(x)$.
Appendix 1
Estimating physiological parameter values
In this appendix we describe our approach to estimating the parameter values summarized in Table 1 from published experimental data.
Glucose diffusivity (D)
Glucose is a small molecule of comparable molecular weight to ATP, which has a diffusion coefficient of 140 µm^{2}/s (Mironov, 2007; Vendelin et al., 2000)
Glucose consumption rate per mitochondrion (k_{g})
The oxidative capacity of muscle mitochondria has been measured at 5.8 mL of O${}_{2}$ per min per mL mitochondria (Harris and Attwell, 2012; Schwerzmann et al., 1989). We assume $6$ glucose molecules are consumed per molecule of oxygen, and a volume of 0.3 µm^{3} for globular mitochondria (Posakony et al., 1977). The corresponding glucose turnover rate of a mitochondrion is then calculated as 1.3 × 10^{5} glucose per second per mitochondrion.
Axon radius (r)
The thickness of mammalian brain axons varies widely from 0.1 µm to 10 µm (Perge et al., 2012). Statistical measurements in the human brain show that most axon diameters fall below 1 µm, with a longtailed distribution of substantially thicker axons (Liewald et al., 2014). We take as our estimate a median diameter of 0.8 µm, which is consistent with measurements in human brain regions (Liewald et al., 2014), in the rat corpus collosum (Barazany et al., 2009), guinea pig retinal neurons (Perge et al., 2009), and in a number of other mammalian tracts (Perge et al., 2012).
Internodal distance (L)
Typical internodal lengths vary widely from 200 μm to 1500 μm (Jacobs, 1988; Court et al., 2004). We use a value of L = 250 µm as measured in the axons of rat anterior medullary velum (Ibrahim et al., 1995).
Mitochondrial density ($\overline{M}$)
Measurements of mitochondrial concentration in human spinal muscular nerves give a linear density of about 15 mitochondria per 100 µm of axon (Xu et al., 2016). Similar densities are observed in Figures 1 and 2 of (Pekkurnaz et al., 2014). Assuming an axonal radius of r ≈ 0.4 µm gives a corresponding density of 0.3 µm^{−3}. EM measurements in rat brain neurons indicate that mitochondria occupy approximately $8\%$ of the neuronal cytoplasmic volume (Pysh and Khan, 1972). Assuming a mitochondrial volume of 0.3 µm^{−3} (Posakony et al., 1977) would give the same density estimate of 0.3 mitochondria per µm^{3}.
Hexokinase MichaelisMenten constant (K_{M})
The MichaelisMenten constant for glucose phosphorylation by the neuronal isoform of hexokinase (HKI) has been measured as K_{M} = 0.03mM (Wilson, 2003).
Ratio of stopped to moving mitochondria (k_{s}/k_{w})
In Pekkurnaz et al., 2014, mammalian neurons grown under high (30 mM) glucose conditions were found to have mitochondria that spent approximately 5% of their time in motion. This fraction should correspond to ${k}_{w}/({k}_{s}+{k}_{w})\approx 0.05$ under our simplified model for mitochondrial motility.
Membrane permeability to glucose ($P,{K}_{MP}$)
The neuronal glucose transporter GLUT3 in rat cerebellar granule neurons has been measured to have a turnover rate of ${k}_{\mathrm{g}\mathrm{l}\mathrm{u}\mathrm{t}3}=853$s${}^{1}$ and a MichaelisMenten constant of $K}_{M,\text{glut3}}=3\phantom{\rule{thinmathspace}{0ex}}\mathrm{m}\mathrm{M$ (Maher et al., 1996). In the same study, the density of GLUT3 channels was measured as 18 pmol/mg cell membrane. We assume the cell membrane has a density of order 1 g/cm^{3} and forms a sheet of thickness 4 nm. This allows us to calculate the area density of GLUT3 channels in cerebellar neurons as approximately a = 43 transporters/µm^{2}.
In the case where the difference in external and internal glucose concentration ($\mathrm{\Delta}G$) is below $K}_{M,\mathrm{g}\mathrm{l}\mathrm{u}\mathrm{t}3$, we can approximate the net flux into the cell as,
allowing an estimate of the permeability $P=\frac{{k}_{\mathrm{g}\mathrm{l}\mathrm{u}\mathrm{t}3}a}{{K}_{M,\text{glut3}}}\approx 0.02\phantom{\rule{thinmathspace}{0ex}}\mu \text{m/s}$
Appendix 2
Effective MichaelisMenten kinetics for glycosylation of Milton
We assume individual steps in glucose metabolism follow classic MichaelisMenten kinetics, with the rate of product formation given by $dP/dt={v}_{i}S/({K}_{Mi}+S)$. We further assume that all pathways considered here are operating in steadystate, with a stationary concentration of all intermediates. When several MichaelisMenten reactions are connected in series (eg: A $\to $B $\to $C), steady state requires that the dependence of final product formation C on the initial reactant A is given by,
where ${v}_{AB},{K}_{AB}$ are the maximum rate and saturation constant for the initial A $\to $B reaction.
If two pathways branch from a single intermediate, as occurs when the hexosamine biosynthetic pathway splits off first from the pentose phosphate pathway and then from glycolysis, then we have an additional reaction B $\to $D that alters the rate of C formation. We make the key assumption that the saturation constants in the first step of both branching pathways are comparable (${K}_{BC}\approx {K}_{BD}$). Steady state then requires
Thus, if the saturation concentrations of the splitting reactions are similar, then the formation of the final product occurs at a rate proportional to the initial substrate consumption, with the same saturation constant. As illustrated in the pathway schematic (Appendix 2—figure 1B), the branching of both the pentose phosphate pathway and glycolysis from the pathway leading to UDPGlcNAc formation involves similar values of the MichaelisMenten constant. We therefore assume that the rate of Milton glycosylation by OGT ($dC/dt$) is proportional to the rate of initial glucose consumption by hexokinase ($dA/dt$). This assumption justifies our use of the same K_{M} for both glucose consumption and mitochondrial stopping. The fraction of metabolic flux funneled into Milton glycosylation is subsumed into the effective rate constant k_{s}.
Appendix 3
Effect of domain length $L$ in uniform permeability model
Dimensional analysis of Equation 6 indicates that the diffusive term for glucose dynamics will be negligible compared to the consumption and entry terms in the limit $L\gg \sqrt{D(G+{K}_{M})/({k}_{g}\overline{M})}$. We assume that external glucose concentrations are well below $10$mM, indicating that the diffusive term is irrelevant for $L\gg 150\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$. If diffusion is neglected, the only length units in the model are found within external glucose and mitochondrial concentrations, both of which are fixed parameters independent of axonal length. We therefore expect in this limit that the model results will not depend on the interval length $L$.
To verify the accuracy of this limit, we plot glucose and mitochondrial distributions for the full model (including diffusion) for interval lengths of $L=100\phantom{\rule{thinmathspace}{0ex}}\mu \mathrm{m}$, 1 mm, and 1 cm (Appendix 3—figure 1a,b). All other parameters are from our physiological estimates in Table 1. We note that for lengths well above $100\mu $m, the distributions are independent of $L$ and are nearly identical to those expected for the model with diffusion excluded. We also plot mitochondrial accumulation and metabolic enhancement in the distal 10% of the interval obtained from solutions of the full model with diffusive transport, which match well to the plots in the main text (Figure 6) that neglect diffusive transport. We can thus assume that glucose entry and turnover are much faster than diffusive spread for biologically relevant parameter regimes.

Appendix 3—figure 1—source data 1
 https://doi.org/10.7554/eLife.40986.025

Appendix 3—figure 1—source data 2
 https://doi.org/10.7554/eLife.40986.026

Appendix 3—figure 1—source data 3
 https://doi.org/10.7554/eLife.40986.027

Appendix 3—figure 1—source data 4
 https://doi.org/10.7554/eLife.40986.028

Appendix 3—figure 1—source data 5
 https://doi.org/10.7554/eLife.40986.029
Data availability
MATLAB code for implementing the models described in this study has been made available on Github: https://github.com/lenafabr/mitoManuscriptCodes (copy archived at https://github.com/elifesciencespublications/mitoManuscriptCodes). Source data files for Figures 3, 4, 5, 6 and the appendix figure are provided in the manuscript and supporting files.
References

Glucose metabolism in nerve terminalsCurrent Opinion in Neurobiology 45:156–161.https://doi.org/10.1016/j.conb.2017.03.007

Mitochondrial dynamics and peripheral neuropathyThe Neuroscientist 14:12–18.https://doi.org/10.1177/1073858407307354

Pattern formation by reactiondiffusion instabilities: application to morphogenesis in DrosophilaJournal of Theoretical Biology 84:629–649.https://doi.org/10.1016/S00225193(80)800245

Facilitated diffusion of glucosePhysiological Reviews 70:1135–1176.https://doi.org/10.1152/physrev.1990.70.4.1135

Mitochondrial trafficking and morphology in healthy and injured neuronsProgress in Neurobiology 80:241–268.https://doi.org/10.1016/j.pneurobio.2006.09.003

The noncompetitive inhibition of brain hexokinase by glucose6phosphate and related compoundsThe Journal of Biological Chemistry 210:597–606.

BookThe Cell: Its Organelles and Inclusions: An Atlas of Fine StructureSaunders (W.B.) Co Ltd.

Activitydependent regulation of surface glucose transporter3Journal of Neuroscience 31:1991–1999.https://doi.org/10.1523/JNEUROSCI.185009.2011

A systematic investigation of the rate laws valid in intracellular environmentsBiophysical Chemistry 124:1–10.https://doi.org/10.1016/j.bpc.2006.04.019

The energetics of CNS white matterJournal of Neuroscience 32:356–371.https://doi.org/10.1523/JNEUROSCI.343011.2012

Dynamic compartmentalization of bacteria: accurate division in E. coliPhysical Review Letters 87:278102.https://doi.org/10.1103/PhysRevLett.87.278102

Relationship between myelin sheath diameter and internodal length in axons of the anterior medullary velum of the adult ratJournal of the Neurological Sciences 133:119–127.https://doi.org/10.1016/0022510X(95)00174Z

Dynamic modulation of intracellular glucose imaged in single cells using a FRETbased glucose nanosensorPflügers Archiv  European Journal of Physiology 456:307–322.https://doi.org/10.1007/s004240070395z

The kinetics of phosphoglucoisomeraseThe Journal of Biological Chemistry 235:2178–2184.

Cellsignalling dynamics in time and spaceNature Reviews Molecular Cell Biology 7:165–176.https://doi.org/10.1038/nrm1838

The oxidative pentose phosphate pathway: structure and organisationCurrent Opinion in Plant Biology 6:236–246.https://doi.org/10.1016/S13695266(03)000396

Control of mitochondrial transport and localization in neuronsTrends in Cell Biology 20:102–112.https://doi.org/10.1016/j.tcb.2009.11.002

Substrate specificity and kinetic parameters of GLUT3 in rat cerebellar granule neuronsBiochemical Journal 315:827–831.https://doi.org/10.1042/bj3150827

Fluctuations in brain glucose concentration during behavioral testing: dissociations between brain areas and between brain and bloodNeurobiology of Learning and Memory 75:325–337.https://doi.org/10.1006/nlme.2000.3976

ADP regulates movements of mitochondria in neuronsBiophysical Journal 92:2944–2952.https://doi.org/10.1529/biophysj.106.092981

Metabolic regulation of mitochondrial dynamicsThe Journal of Cell Biology 212:379–387.https://doi.org/10.1083/jcb.201511036

Kinetics of GLUT1 and GLUT4 glucose transporters expressed in Xenopus oocytesThe Journal of Biological Chemistry 268:8514–8520.

Dynamic reorganization of metabolic enzymes into intracellular bodiesAnnual Review of Cell and Developmental Biology 28:89–111.https://doi.org/10.1146/annurevcellbio101011155841

Regional differences in brain glucose content in graded hypoglycemiaNeurochemical Pathology 5:131–142.

How the optic nerve allocates space, energy capacity, and informationJournal of Neuroscience 29:7917–7928.https://doi.org/10.1523/JNEUROSCI.520008.2009

Why do axons differ in caliber?Journal of Neuroscience 32:626–638.https://doi.org/10.1523/JNEUROSCI.425411.2012

Mitochondrial growth and division during the cell cycle in HeLa cellsThe Journal of Cell Biology 74:468–491.https://doi.org/10.1083/jcb.74.2.468

Regulation of mesangial cell hexokinase activity by PKC and the classic MAPK pathwayAmerican Journal of PhysiologyRenal Physiology 277:F742–F749.https://doi.org/10.1152/ajprenal.1999.277.5.F742

Multiple functions of the paranodal junction of myelinated nerve fibersJournal of Neuroscience Research 87:3250–3258.https://doi.org/10.1002/jnr.22013

Drosophila Miro is required for both anterograde and retrograde axonal mitochondrial transportJournal of Neuroscience 29:5443–5455.https://doi.org/10.1523/JNEUROSCI.541708.2009

The role of myelin and oligodendrocytes in axonal energy metabolismCurrent Opinion in Neurobiology 23:1065–1072.https://doi.org/10.1016/j.conb.2013.09.008

Mitochondrial trafficking in neuronsCold Spring Harbor Perspectives in Biology 5:a011304.https://doi.org/10.1101/cshperspect.a011304

Hypoxic regulation of glycosylation via the Nacetylglucosamine cycleJournal of Clinical Biochemistry and Nutrition 48:20–25.https://doi.org/10.3164/jcbn.11015FR

Energetic basis of brain activity: implications for neuroimagingTrends in Neurosciences 27:489–495.https://doi.org/10.1016/j.tins.2004.06.005

Supply and demand in cerebral energy metabolism: the role of nutrient transportersJournal of Cerebral Blood Flow & Metabolism 27:1766–1791.https://doi.org/10.1038/sj.jcbfm.9600521

The facilitative glucose transporter GLUT3: 20 years of distinctionAmerican Journal of PhysiologyEndocrinology and Metabolism 295:E242–E253.https://doi.org/10.1152/ajpendo.90388.2008

SoftwareMATLAB, version 8.6.0The MathWorks Inc, Natick, Massachusetts.

Coordinated regulation of glutamine:fructose6phosphate amidotransferase activity by insulin, glucose, and glutamine. Role of hexosamine biosynthesis in enzyme regulationThe Journal of Biological Chemistry 266:10148–10154.

Chronic insulin hypoglycemia induces GLUT3 protein in rat brain neuronsAmerican Journal of PhysiologyEndocrinology and Metabolism 272:E716–E719.https://doi.org/10.1152/ajpendo.1997.272.4.E716

Regulation of energy metabolism in Trypanosoma (Schizotrypanum) cruzi epimastigotes. I. Hexokinase and phosphofructokinaseMolecular and Biochemical Parasitology 11:225–239.https://doi.org/10.1016/01666851(84)900689

Chemotaxis: signalling the way forwardNature Reviews Molecular Cell Biology 5:626–634.https://doi.org/10.1038/nrm1435

Regulation of mitochondrial respiration in heart cells analyzed by reactiondiffusion model of energy transferAmerican Journal of PhysiologyCell Physiology 278:C747–C764.https://doi.org/10.1152/ajpcell.2000.278.4.C747

Isozymes of mammalian hexokinase: structure, subcellular localization and metabolic functionJournal of Experimental Biology 206:2049–2057.https://doi.org/10.1242/jeb.00241

Metabolic pathway compartmentalization: an underappreciated opportunity?Current Opinion in Biotechnology 34:73–81.https://doi.org/10.1016/j.copbio.2014.11.022

Glucose Deprivationinduced Increase in Protein O GlcNAcylation in cardiomyocytes is CalciumdependentJournal of Biological Chemistry 287:34419–34431.https://doi.org/10.1074/jbc.M112.393207
Decision letter

Raymond E GoldsteinReviewing Editor; University of Cambridge, United Kingdom

Naama BarkaiSenior Editor; Weizmann Institute of Science, Israel
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Spatial control of neuronal metabolism through glucosemediated mitochondrial transport regulation" for consideration by eLife. Your article has been reviewed by two peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Naama Barkai as the Senior Editor. The following individual involved in review of your submission has agreed to reveal his identity: Timothy A Ryan (Reviewer #1).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
This manuscript describes a general framework for examining the quantitative details of how the feedforward control of mitochondrial motility by a downstream sensing of glucose metabolism (through OGlcNAc modification) might lead to both changes in the spatial distribution of mitochondria and the potential enhancement of overall metabolic flux that would ensue. The referees agreed that the manuscript was interesting, wellwritten, and timely.
Essential revisions:
1) The treatments seem to treat the OGT modification as irreversible as no terms describing the reverse reaction that would be driven by an O0GlcNAcase.
2) The authors also only consider two possible fates of Glucose 6P (glycolysis and the hexosamine pathway). The pentose phosphate pathway is another possible shunt that should probably at the very least be discussed or its omission justified.
3) A forward Euler scheme was used to solve the reactiondiffusion Equation 1. Such schemes are known to be generally unstable, so the authors need to justify their use and/or redo the numerics with a proper scheme (e.g. CrankNicholson) to confirm the results.
https://doi.org/10.7554/eLife.40986.031Author response
Essential revisions:
1) The treatments seem to treat the OGT modification as irreversible as no terms describing the reverse reaction that would be driven by an O0GlcNAcase.
We thank the reviewer for this comment and agree that this was not clearly explained in the manuscript. In our simplified quantitative model, we assume that all pauses are associated with an OGlcNAcylation event. Recovery from the pause at the constant rate k_{w} represents the removal of the modification of Milton through the activity of the complementary enzyme OGlcNAcase. We have added a paragraph to our manuscript more thoroughly explaining the rationale behind our simplified model for mitochondrial dynamics, quoted below:
“It has recently been demonstrated that the key motor adaptor protein (Milton) is sensitive to glucose levels, halting mitochondrial motility when it is modified through OGlcNAcylation by the OGT enzyme [Pekkurnaz et al., 2014]. […] Other sources of pausing, such as Ca^{2+}regulated motor disengagement, PINK1/Parkinmediated detachment of motors, and anchoring to the microtubules by syntaphilin [Schwarz, 2013], are not considered here in order to focus specifically on the effect of glucosedependent mitochondrial spatial organization.”
Furthermore, we include an additional paragraph in the Discussion pointing out simplifications in this model that ignore potential glucosedependent regulation of other steps in the metabolic pathway connecting glucose and mitochondrial stopping:
“Additional metabolic feedback loops, not included in our model, may result in a more complex dependence of mitochondrial stopping on glucose concentration. […] These regulatory mechanisms provide additional potential routes of metabolic control through mitochondrial positioning.”
2) The authors also only consider two possible fates of Glucose 6P (glycolysis and the hexosamine pathway). The pentose phosphate pathway is another possible shunt that should probably at the very least be discussed or its omission justified.
We thank the reviewer for pointing out this omission. We have now included the metabolic branch for the pentose phosphate pathway in our schematic in Figure 1C. We note that pathway branches where the saturation constants are similar to the corresponding step in the HBP, will not affect the overall kinetics in our model, except by changing the fraction of glucose that flows down the HBP pathway. This fraction is subsumed in the effective pausing rate k_{s}. We have modified the relevant paragraph in the model development section to address this issue:
“Upon entry into the cell, the first ratelimiting step of glucose metabolism is its conversion into glucose6phosphate by hexokinase. […] We therefore model the kinetics of Milton modification using the same MichaelisMenten form as for hexokinase activity, with the pathway flux leading to Milton modification subsumed within a rate constant for mitochondrial stopping (k_{s}).”
We also include a detailed justification of using the same saturation constant K_{M} for both glucose consumption and OGlcNAcylation of Milton in Appendix 2. A schematic with appropriate citations has been added to Appendix 2, listing the saturation constants of both the pentose phosphate and glycolysis branching pathways.
3) A forward Euler scheme was used to solve the reactiondiffusion Equation 1. Such schemes are known to be generally unstable, so the authors need to justify their use and/or redo the numerics with a proper scheme (e.g. CrankNicholson) to confirm the results.
We thank the reviewer for the comments, and agree that forward Euler methods have poor convergence properties. However, such methods can be used for simple systems, provided a sufficiently small timestep is used to assure convergence (smaller than the fastest timescale in the system, including the scale set by the spatial discretization). We know convergence is achieved in our model because the magnitude of the time derivative goes to zero as the simulation approaches steady state. To clarify this issue, we have added additional description of our choice of timestep and convergence criteria in the Materials and methods section, under “Discrete mitochondria simulations”:
“Because forward Euler methods have stringent conditions for stability and convergence, we use a timestep that is much smaller than both the glucose decay timescale and the timescale associated with diffusion over our spatially discretized grid (see below).”
“We integrate the simulation forward in timesteps of δt=0.2∆x2D, where ∆x is the spatial discretization. This timescale is much smaller than the relevant decay time for glucose consumption τg=kgMKM1. Using these small timesteps allows for stability and robust convergence with the forward Euler method.”
In the section titled “Numerical solution for steadystate distributions with localized glucose entry” we note the convergence criterion for the continuum model: “The distributions are assumed to be converged once the root mean squared rate of glucose change drops below the minimal cutoff: 10^{−6}k_{g}M.”
https://doi.org/10.7554/eLife.40986.032Article and author information
Author details
Funding
National Institutes of Health (R35GM128823)
 Gulcin Pekkurnaz
University of California, San Diego (Chancellor's Research Excellence Scholarship)
 Anamika Agrawal
Alfred P. Sloan Foundation (FG201810394)
 Elena F Koslover
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Saurabh Mogre for fruitful discussions, Manho Tang for numerical modeling advice, and David Kleinfeld for helpful comments on the manuscript.
Senior Editor
 Naama Barkai, Weizmann Institute of Science, Israel
Reviewing Editor
 Raymond E Goldstein, University of Cambridge, United Kingdom
Publication history
 Received: August 10, 2018
 Accepted: December 17, 2018
 Accepted Manuscript published: December 18, 2018 (version 1)
 Version of Record published: January 7, 2019 (version 2)
Copyright
© 2018, Agrawal et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 3,697
 Page views

 515
 Downloads

 13
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Cell Biology
 Physics of Living Systems
Motile cilia are hairlike cell extensions that beat periodically to generate fluid flow along various epithelial tissues within the body. In dense multiciliated carpets, cilia were shown to exhibit a remarkable coordination of their beat in the form of traveling metachronal waves, a phenomenon which supposedly enhances fluid transport. Yet, how cilia coordinate their regular beat in multiciliated epithelia to move fluids remains insufficiently understood, particularly due to lack of rigorous quantification. We combine experiments, novel analysis tools, and theory to address this knowledge gap. To investigate collective dynamics of cilia, we studied zebrafish multiciliated epithelia in the nose and the brain. We focused mainly on the zebrafish nose, due to its conserved properties with other ciliated tissues and its superior accessibility for noninvasive imaging. We revealed that cilia are synchronized only locally and that the size of local synchronization domains increases with the viscosity of the surrounding medium. Even though synchronization is local only, we observed global patterns of traveling metachronal waves across the zebrafish multiciliated epithelium. Intriguingly, these global wave direction patterns are conserved across individual fish, but different for left and right nose, unveiling a chiral asymmetry of metachronal coordination. To understand the implications of synchronization for fluid pumping, we used a computational model of a regular array of cilia. We found that local metachronal synchronization prevents steric collisions, cilia colliding with each other, and improves fluid pumping in dense cilia carpets, but hardly affects the direction of fluid flow. In conclusion, we show that local synchronization together with tissuescale cilia alignment coincide and generate metachronal wave patterns in multiciliated epithelia, which enhance their physiological function of fluid pumping.

 Immunology and Inflammation
 Physics of Living Systems
One of the feats of adaptive immunity is its ability to recognize foreign pathogens while sparing the self. During maturation in the thymus, T cells are selected through the binding properties of their antigenspecific Tcell receptor (TCR), through the elimination of both weakly (positive selection) and strongly (negative selection) selfreactive receptors. However, the impact of thymic selection on the TCR repertoire is poorly understood. Here, we use transgenic Nur77mice expressing a Tcell activation reporter to study the repertoires of thymic T cells at various stages of their development, including cells that do not pass selection. We combine highthroughput repertoire sequencing with statistical inference techniques to characterize the selection of the TCR in these distinct subsets. We find small but significant differences in the TCR repertoire parameters between the maturation stages, which recapitulate known differentiation pathways leading to the CD4+ and CD8+ subtypes. These differences can be simulated by simple models of selection acting linearly on the sequence features. We find no evidence of specific sequences or sequence motifs or features that are suppressed by negative selection. These results favour a collective or statistical model for Tcell self nonself discrimination, where negative selection biases the repertoire away from self recognition, rather than ensuring lack of selfreactivity at the singlecell level.