Cable energy function of cortical axons

Accurate estimation of action potential (AP)-related metabolic cost is essential for understanding energetic constraints on brain connections and signaling processes. Most previous energy estimates of the AP were obtained using the Na+-counting method, which seriously limits accurate assessment of metabolic cost of ionic currents that underlie AP conduction along the axon. Here, we first derive a full cable energy function for cortical axons based on classic Hodgkin-Huxley (HH) neuronal equations and then apply the cable energy function to precisely estimate the energy consumption of AP conduction along axons with different geometric shapes. Our analytical approach predicts an inhomogeneous distribution of metabolic cost along an axon with either uniformly or nonuniformly distributed ion channels. The results show that the Na+-counting method severely underestimates energy cost in the cable model by 20–70%. AP propagation along axons that differ in length may require over 15% more energy per unit of axon area than that required by a point model. However, actual energy cost can vary greatly depending on axonal branching complexity, ion channel density distributions, and AP conduction states. We also infer that the metabolic rate (i.e. energy consumption rate) of cortical axonal branches as a function of spatial volume exhibits a 3/4 power law relationship.

this equation describes the AP-related energy consumption rate per unit membrane area (cm −2 s −1 ) at any axonal distance and any time. The individual terms on the right-hand side of the equation represent the contributions of the sodium, potassium, leak, and axial currents, respectively. Calculations based on this function distinguish between the contributions of each item toward total energy consumption.

Figure 1. Effect of Axonal Length on AP-related Energy Consumption and Efficiency. (A) Cable model of a
Hodgkin-Huxley-type cortical axon, where axial current, i a , flows through axial resistance, r a , within a uniform cylinder. The membrane currents consist of i C , i K , i Na , and i L , through c m , g K , g Na and g L , respectively. V K , V Na and V L , Nernst potentials. APs (60 Hz) were initiated by a stimulating current (19.1 μ A/cm 2 , 1 s) at one end of the uniform axon (X = 0 μ m) at 37 °C.(B) Traces of AP (black), -i Na (red, inverted Na + current), i K (blue) and i a (green) in the compartments of a uniform axon as in (A) (1.5 μ m in diameter, 1000 μ m in axonal length, 50 μ m in compartmental length (Δ x), unmyelinated) at X = 0, 200 and 500 μ m. (C) Traces of the total energy consuming power (black, P tot = P Na + P K + P a + P L ) and its main components, i.e., the energy consuming power of the sodium conductance (red, P Na = i Na × (V − V Na )), potassium conductance (blue, P K = i K × (V − V K )), and axial conductance (green, P a = π ∂ ∂ i a a V x 1 2 ), in the compartments at X = 0, 200, and 500 μ m. P L is the leak conductance (negligible). See also Figure S2. (D) Distribution of energy cost per unit membrane area per AP along axons of different lengths (same diameter, APs firing at 60 Hz, at 37 °C). Note that the longer the axon, the higher the energy cost at the same distance. Inset, increase of the energy cost at the AP initiation site with axonal length. (E) Distribution of the excess Na + entry ratio (γ ), i.e., the ratio of the total Na + flux during an AP to the theoretically minimal Na + load needed to generate the upstroke of an AP, along axons of different length. Note that the longer the axon, the higher the γ value at the same distance. Inset, increase of the value of γ at the AP initiation site with axonal length. (Theoretical minimum Na + load: Na + flux integrated from dv/dt = 20 to dv/dt = 0.) (F) Comparison of the energy calculation based on Na + counting (E Na c ) to the one based on the energy function (E). See also Figure S2.
Scientific RepoRts | 6:29686 | DOI: 10.1038/srep29686 Effect of Axonal Length on AP-related Energy Consumption and Efficiency. Next, we applied the above energy function to estimate the distribution of energy cost along a single-cable axon (Fig. 1A). In the HH-type cortical axon model that is based on experimental data 18,30,31 , AP propagation (60 Hz, at 37 °C) was simulated by injecting a steady current (19.1 μ A/cm 2 , 1 s) at one end of the axon ( Fig. 1A; simulations in this paper were run in MATLAB). From recordings of APs, transmembrane ion currents, and axial currents at different axonal distances (Fig. 1B), we calculated the energy consumption associated with each ionic current (P ion = i ion × (V − V ion ), Fig. 1C, Figure S2) as well as the total energy cost ∫ = ∑ (E P dt) ion as a function of axon length (see Fig. 1D). Although the Na + current (i Na ) amplitude during an AP increased with distance, the width of i Na decreased (Fig. 1B, red), resulting in a reduction in total Na + influx and the energy consumed by the sodium conductance (integration of the metabolic consuming power (P Na = i Na × (V − V Na )) of the sodium conductance over the duration of an AP, Fig. 1C, red). In contrast, the K + current (i K ) almost leveled off along the axon (Fig. 1B, blue), making the metabolic consuming power of the potassium conductance (P K = i K × (V − V K )) stable (Fig. 1C, blue). Additionally, the amplitude of the axial current (i a ) rose sharply along the axon, resulting in a considerable increase in the amount of energy consumed by the axial conductance (integration of = π ).
To study the effect of axon length on energy consumption, we carried out simulations of unbranched cable axons ranging from 0 (point neuron model) to 1,500 μ m in length (with the same diameter, AP firing at 60 Hz, at 37 °C). As shown in Fig. 1D, the energy cost along an axon of any length is generally highest at the AP initiation point, decreases sharply in the initial 200 μ m, increases slowly and finally saturates at a value below that of the initial section. More crucially, the longer the axon, the higher the energy cost per unit membrane area per AP at the same distance (Fig. 1D). Even more importantly, in the cable model, the cost at the AP initiation site (E int ) can be up to 15% higher than the cost in the point model (Fig. 1D, inset), indicating that more energy is needed to promote AP conduction across longer axons than across shorter ones. By analyzing the components that contribute to the energy consumption, we found that the higher cost in the longer axons was due to the higher cost of the axial and sodium conductances (see Fig. 1B,C), which is reasonable because AP propagation in longer axons demands a larger Na + influx and a larger axial current.
We also examined how axonal length affects energy efficiency, as measured by the excess Na + entry ratio (γ ). A γ value greater than 1 indicates excess Na + influx, and larger γ values, in turn, reflect a less efficient use of energy. As illustrated in Fig. 1E, the AP initiation location of the axon had the lowest efficiency, with values that ranged from 1.5 to 2; the longer the axon, the lower the efficiency at the same distance. In particular, the decrease in efficiency at the initiation site was up to 33% greater than the value of 1.5 obtained using the point model (note a value of 1.3 is obtained using the Na + counting method), see Fig. 1E inset, indicating that long axons will pay the price of reduced efficiency.
To compare the energy-function approach with the traditional Na + -counting method, we applied both approaches to the cable model and the point model. Figure 1F shows that the cost estimated using Na + counting was significantly lower in both a cable model with a length of 450 μ m (27% lower) and in a point model of the same length (42% lower), suggesting that the Na + -counting method severely underestimates metabolic consumption. The discrepancy of the two results is calculated as− , where E Na c and E are the results of the sodium counting method and our approach, respectively. Our results show that factors including spiking rate, length of cable, temperature and gmax have limited effects on the discrepancy (< 10%), indicating that the discrepancy is quite stable. See Discussion for detailed results.

Effect of AP Conduction State on AP-related Energy Consumption. Neurons exhibit various axonal
branching patterns. Different AP conduction states (i.e., successful, blocked, and reflected conduction) have been observed at branch points in different types of neurons 32 . Factors that affect AP propagation through branch points include the geometry of the branches connected at the branch point [33][34][35] , the AP spike rate 32,36,37 and the temperature 38 . To compare the energy costs of different AP conduction states, we examined a model consisting of a parent axon connected with a pair of identical branches ( Fig. 2A). The parent axon is stimulated at the proximal end so that it generates APs at a controlled spike rate (5 or 60 Hz) at 37 °C. We determined the smallest child diameter at which AP propagation failure occurred for at least one AP in the train. Consider the geometrical ratio (GR) 39,40 , which is defined as where d D and d M are the diameters of the child branches and the parent axon, respectively. Figure 2B illustrates the influence of the GR on APs (top) and energy consuming power (bottom) at the initiation site (a) and at 25 μ m after the branch point (b) for a firing rate of 5 Hz. Successful propagation into child branches was observed when the GR was ≤ 9; however, when the GR reached 10, conduction failed at the branch point, indicating that the critical GR (GR c ), or the largest GR for successful conduction, is between 9 and 10. Similar effects of GR on propagation were observed with an AP spike rate of 60 Hz but with a lower GR of 7 (Fig. 2C). Figure 2D,E show the energy cost along the axon for two spike rates. For AP conduction through the branch point, less energy is used when the GR is small, while a larger GR generally requires more energy during the transition. For GR values greater than the critical value GR c , AP propagation failed in the child branch, and very little energy is consumed by the child branches. Actually, the rapid change in energy cost around the branch point at GR ≠ 1 (blue, red, and green) was mainly caused by AP reflection around the branch point. This phenomenon was not observed at GR = 1 because at that time, the APs were propagating in a branching system equivalent to a homogeneous parent axon in accordance with Rall's 3/2 power law 40 . Our simulations predict that the GR c decreases nonlinearly with increasing AP firing rate (Fig. 2F).
Scientific RepoRts | 6:29686 | DOI: 10.1038/srep29686 Distribution of AP-related Energy Consumption in an Axonal Branching Tree. We next applied the energy-function method to a more complexly branched axon whose structure was that of a binary tree with GR = 1 (Fig. 3A). We found that the energy cost per unit length per AP decreased monotonically with distance from the injection site (Fig. 3A) and BL (Fig. 3B). Nonetheless, the energy cost per unit membrane area decreased sharply in the initial 400 μ m before rising and then leveling off (Fig. 3C), resembling the distribution along a single-cable axon (Fig. 1D) and indicating that this distribution pattern of energy cost per unit membrane area against distance might apply to all types of axonal morphology.
Then, to examine how the energy consumption of an entire axonal branching system varies with geometry (i.e., axonal volume, membrane surface and branching complexity), we changed the geometry by increasing the total BL of the system (BL sys ) from one to four while keeping the diameter of the highest level of child branch fixed at 0.2 μ m (geometry versus BL sys , see Figure S1A-S1C). Interestingly, the log-log plot (Fig. 3D) shows that the dependence of the system's metabolic rate on its volume followed an allometric scaling law with a scaling exponent 0.75 for GR = 1 here (the exponent changes from 0.3 to 2 when GR changes from 0.5 to 2). Assuming that the density of the axon (1.05 g/mL) is uniform, we found that the allometric scaling of the metabolic rate linked to AP versus the mass of the axonal tree ( Figure S1E) is the same as the empirically observed scaling of the metabolic rate versus the mass of the entire organism 23,24 . The excellent match between our modeling and the experimental data validates the theory that uses a fractal structure model to explain the quarter-law scaling that is common in biological systems 22 . Additionally, this allometric scaling arose only when we changed the axonal tree's volume by varying the total BL of the tree and keeping the diameter of the highest level of child branch unchanged (the results generated when the volume was changed in other ways are described in the DISCUSSION), suggesting that there may be a lower limit for the diameter of real axonal branches in any type of neuron, consistent with previous anatomical findings and computational predictions 41 ; a comprehensive survey of anatomical data showed that the lower limit for AP-conducting axons is 0.08-0.2 μ m in diameter, and stochastic simulations demonstrated that due to channel noise, the limiting diameter for mammalian pyramidal cell axons is 0.15 μ m.
Moreover, the scaling exponent for the relationship between energy consumption and membrane surface area is close to 1 rather than 3/4 ( Fig. 3E). This finding is in agreement with results from computational models of single-compartment spiking neurons 7 . This relationship can be attributed to the uniform spread of ion conductances over the membrane surface of the system. Additionally, our results showed that during AP propagation, the metabolic cost of the entire axonal branching system increased exponentially with BL sys ( Figure S1D).

Effects of Ion Channel Density and Distribution on AP-related Energy Consumption. This
section investigates how AP-related energy consumption is influenced by the density and distribution of ion channels. First, to study the impact of ion channel density, we changed the Na + and K + channel density of the and unreliable (right), respectively. In (B), AP spike rate of 5 Hz; GR ≤ 9, reliable state; GR ≥ 10, unreliable state. In (C), AP spike rate of 60 Hz; GR ≤ 6, reliable state; GR ≥ 7, unreliable state; red asterisk indicates propagation failure. (D,E) Variation in energy cost per unit membrane area per AP along the axon (branch point is at X = 0) at four different GR values when the AP spike rate is 5 Hz (D) and 60 Hz (E). GR c represents the critical GR, i.e., the largest GR for successful conduction. Note that in the child branch (X > 0), the energy cost per unit membrane area per AP was much lower in unreliable propagation states. (F) GR c decreases exponentially with AP frequency.
Scientific RepoRts | 6:29686 | DOI: 10.1038/srep29686 single-cable axon (Fig. 1A) by manipulating the maximal Na + and K + channel conductance (g Na max from 50 to 650 mS/cm 2 , g k max from 3 to 100 mS/cm 2 , uniformly distributed over the axon) and then stimulated APs (60 Hz, 37 °C) that propagate along the axon. As illustrated in Fig. 4A-C, at the initiation site of the axon, when g k max was fixed, larger g Na max led to increased Na + entry during the rising phase of the AP and increased K + exit during the falling phase of the AP, increasing energy costs. A similar situation was observed when g Na max was fixed and g k max was changed (Fig. 4A,D,E). Because the combination of smaller g Na max and smaller g k max lowered energy costs and increased energy efficiency simultaneously (Fig. 4F,G), we propose that within the range of ion-channel densities sufficient for AP generation and propagation, there is an optimal combination of ion channel densities that will result in both minimal energy consumption and maximal efficiency. Actually, the effect of ion channel density found in the cable model is consistent with previous results obtained through the Na + -counting method for point neuronal models, including biophysical modeling and dynamic clamping of neocortical fast-spiking interneurons 8 , as well as single-compartment models of the squid giant axon and rat interneurons 7,42 . Next, to investigate the influence of ion channel distribution, we carried out simulations in the unbranched axon with a nonuniform Na + channel expression pattern (Fig. 5A, top) and then compared the results with those of the uniform axon. In the nonuniform axon, the initial 50 μ m, where g Na max is three times larger than the value for the rest of the axon, showed a dramatic decrease in Na + entry (Fig. 5B, red, note the amplitude) and K + exit during an AP (blue). This finding accounts for the sharp decrease in two of the main components of energy expenditure, E Na (62% drop) and E K (33% drop) (see Fig. 5A, bottom; Fig. 5C), and consequently, the significant decrease in total energy expenditure (36% drop, see Fig. 5C,D, black). Noticeably, both the energy cost (Fig. 5D, black) and the excess Na + entry ratio (Fig. 5F, black) reached minimum values right at the end of the high Na + density area. Importantly, comparing the energy costs along the nonuniform axon and the uniform one, based on either the energy-function (Fig. 5D) or Na + -counting (Fig. 5E) method reveals that the nonuniform Na + channel pattern resulted in an energy savings of approximately 20% for the total axon. Similarly, comparing the excess Na + entry ratios (Fig. 5F) shows that the nonuniform Na + channel pattern enhanced energy efficiency by approximately 10%, , (E), g K max from 3 to 100 mS/cm 2 , g Na max kept the same. (F,G) Color map of the energy cost per unit membrane area per AP at the initiation site (F) and the excess Na+ entry ratio at the initiation site (G) as a function of g Na max and g K max . White: low cost or low ratio. Note that the minimum energy cost (F) and the maximum energy efficiency (G) were achieved when both the g Na max and g K max were smallest.
with a saturated Na + entry-ratio value of approximately 1.5 compared with the uniform pattern (in gray line). The low energy usage and high energy efficiency of the nonuniform pattern concluded here by the simulations could be a novel prediction for further experimental investigation. It actually suggests a nontrivial rule that actual neuronal axons may have nonuniform distributed ionic channel distributions in an energy efficient strategy.

Discussion
Based on the energy function used to estimate the metabolic cost of APs in HH point neuron models [19][20][21] , we established the energy-function method for the cable model of axons. Through this approach, we provided the first demonstration that AP propagation requires 15% more energy at the initiation site (Fig. 1D, inset) than the amount predicted by the point model. In addition, the allometric scaling relationship between the total energetic rate of the axonal branching tree (P tot ) versus the axonal volume (V), ∝ . P V tot 0 75 , suggested that the metabolic rate of a single neuron is very likely to scale to the 3/4 power of its mass, just as the metabolic rate of an entire organism scales with its body mass in animals 23 , plants and microbes 24 . Additionally, the allometric scaling relationship suggests an invariant minimum diameter for axonal branches in any type of neuron; this suggestion is supported by anatomical findings and stochastic simulations that take channel noise into account 41 . Moreover, by enabling the energy costs to be precisely distributed over the complicated branching pattern (Fig. 3A), this method can be applied to branched dendrites as well.
Advantages of the Cable Energy Function. The cable energy function derived here provides an accurate method for calculating the energy cost of electric signals conducted in any type of neuron with an axonal or dendritic arbor. This theoretical framework allows us to estimate the energy used by cortical axons with many types of ion channels more accurately than is possible with the Na + -counting method, especially in the presence of some subtypes of Ca 2+ and K + channels. By taking into account all of the energy-consuming components in the cable equation, including the axial and transmembrane leak currents, the energy function provides a precise evaluation at any axonal distance and any time (see METHODS). The distribution of energy usage among all of the energy-consuming components in the cortical neuron shows that potassium and sodium conductances consumed the most energy; together, these conductances accounted for 84.4% of the total energy usage in the cable model ( Figure S2) and 96% in the point model. This finding is comparable to the result of the point model of the squid giant neuron 19 .
In contrast, the Na + -counting method tends to underestimate energy consumption (Fig. 1F). In fact, because it is based on the number of ATP molecules required by Na + /K + -ATPase pumps to restore the Na + and K + concentration gradients after an AP, this method merely calculates the energy cost related to the transmembrane and nonuniform (black) Na + channel distribution along the cortical axon. Bottom, variation of the main components of the total energy cost per unit membrane area per AP, i.e., the energy costs of the K + (blue), Na + (red) and axial (green) conductances, along the nonuniform axon when APs propagated at 60 Hz at 37 °C. (B) Traces of AP (black), -i Na (red, inverted Na + current), i K (blue) and i a (green) at X = 0, 50 and 500 μ m. (C) Traces of the total energy consuming power (black) and its main components, i.e., the energy consuming power of the sodium (red), potassium (blue), and axial (green) conductances at X = 0, 50 and 500 μ m. (D,E) Comparison of the energy cost per unit membrane area per AP along axons with nonuniform (black) and uniform (grey) Na + channel expression patterns. The energy cost was obtained by the energy-function (D) and Na + -counting (E) methods. (F) Comparison of the excess Na + entry ratio along axons with nonuniform (black) and uniform (grey) Na + channel expression patterns. Note that in (D-F), the values of the nonuniform axon were lower than those of the uniform axon, except in the initial 50 μ m.
Scientific RepoRts | 6:29686 | DOI: 10.1038/srep29686 flowing of Na + and K + ions. However, even when the energy usage contributed by sodium and potassium channels only was considered for the energy-function method, the energy usage estimate derived from the Na + -counting method was 20% lower for the cortical neuron cable model and 32% lower for the point model. More importantly, we have also examined how the discrepancy between the two methods is influenced by the parameters, i.e. spiking rate, length of cable, temperature and g max . The discrepancy is calculated as , where E Na c and E are the results of the sodium counting method and our approach, respectively. When one factor was examined, other factors were kept the same. Our results show that these factors have limited effects on the discrepancy (< 10%), indicating that the discrepancy is quite stable. The detailed results are as follows. First, when the spiking rate increased from 5 Hz to 125 Hz, the discrepancy of the point model and the cable model increased 7% (39~46%) and 1% (37~38%), respectively. Second, when the length of the cable was ≤ 3000 um, the discrepancy was within the range of 27~38%. Third, when the temperature was in the range of 6 °C~42 °C, the discrepancy was within the range of 32~37%. Fourth, larger g K max and g Na max leads to larger discrepancy, and within the range of ≤ g 100 mS cm / K max 2 and ≤ g 650 mS cm / Na max 2 , the discrepancy is within the range of 30~40%.
For a deeper understanding of how the two methods calculate energy differently, see Methods. In addition, the energy-function method predicts some interesting phenomena. One example is the nonuniform Na + channel expression pattern, which has been observed experimentally and was demonstrated to be vital for axonal signaling functions in previous modeling studies 9,43 . The natural beauty of the computational modeling is that it could provide novel prediction to instruct new experiment investigation and confirmation. The model conclusion of high energy efficiency with nonuniform ionic channel distribution could be a nontrivial point for new experimental examination. Other predictions, such as the prediction that longer axonal length leads to higher energy consumption (Fig. 1D) and the prediction that there is an optimal combination of ion channel densities that will minimize metabolic cost (Fig. 4F,G), need to be experimentally validated. Furthermore, this analytical approach can be used to make predictions regarding a diverse range of neurons because it can be applied to any HH-type neuron model, such as one where Ca + channels play an important role in AP generation, simply by deriving the energy function for that model.

Allometric Scaling of Energy Versus Volume. Biological systems have evolved branching networks for
information flow and energy transport. Experimental 23,[44][45][46] and theoretical studies 22,47,48 have shown that the metabolic rate of a given biological system generally scales with the 3/4 power of system volume or mass, and the potential mechanism predicted a space-filling, fractal-like branching pattern with a minimized energy distribution 22 . In this paper, we also examined the relationship between the metabolic rate and volume of cortical axons. To do this, the diameter of the final branch level was set to an invariant value of 0.2 μ m, as reported by experiments. This setting ensured a realistic value for the diameters of axonal branches in pyramidal cells 41 ; these values varied from 0.2 to 1.27 μ m (Fig. 3A, inset, note that GR = 1 for the diameter of branches at adjacent BLs). The model analysis predicted that allometric scaling with a scaling exponent of 0.75 arose only when the tree's volume was changed by varying the branch level while keeping the diameter of the final branch level constant. When we kept the diameter of the lowest level of parent branch constant in trees of different branch level ( Figure S1F) or when we changed the volume of the single uniform axon by changing the length of each branch ( Figure S1G), the scaling relation could not be obtained. Interestingly, our setting of the model is consistent with the main assumption of the general theory used to explain the widely used quarter-power scaling, i.e., keeping the final branch of the fractal network a size-invariant unit 22 .
The power law relation between the energy consumption of the whole axonal architecture and the level of branching ( Figure S1D) was also observed in computational models of the bifurcating axonal arbor of dopamine neurons of the substantia nigra pars compacta 49 . The high energy demand of the massive axonal arbors of these neurons, which was confirmed by the computational work, was proposed as a factor in the selective vulnerability of these dopamine neurons to Parkinson's disease 49 . Implications of the Analytical Approach. With extended applications from axons to more complex neuronal systems, this analytical approach could play a significant role in the estimation of energy expenditure, from subcellular to whole-organism levels. For instance, by applying this approach to the morphology-based HH neuron model that was built on experimental data, we can obtain a more realistic subcellular distribution of energy use, compared to the one that was first obtained through calculations of Na + influx 1 . Similarly, the energy function could be used to replace previous estimates of the elevation of AP-related energy, which is highly dependent on the excess Na + entry ratio, within the neural signaling energy budgets that have been constructed for grey matter 2,4 and white matter 3 .
Furthermore, this analytical method of energy calculation can be used to investigate factors that influence energy consumption and reveal trade-offs between energetic constraints and signaling performance. In this study, we examined the impact of axonal geometry, AP propagation state, and ion channel density. Other factors that strongly affect AP conduction speed and spiking frequency [7][8][9][10][11]16 , such as information rate 6,7,50 and ion channel kinetics, are worth re-investigating using the analytical approach to confirm previous conclusions.
Moreover, this analytical approach provides a valuable way of understanding energetic constraints on neuronal morphology. For instance, the repression of axonal branching, which is a crucial regulatory mechanism that can counteract positive cues to evoke branching 51 , might also be an energy-saving strategy because metabolic costs linearly increase with increasing surface area of the axon branching system (Fig. 3E). In addition, how the growth of an axonal process is coordinated with the retraction of another is unknown on a molecular level 51 , and from an energy consumption standpoint, the repression of axonal branching might help to achieve physiological function within the limits of available energy resources. Furthermore, branch points, which are thought to represent a powerful process that filters communication with postsynaptic neurons 32 , might also have an energy-conserving function. Our results, which demonstrate a reduction in metabolic cost when conduction fails compared to when it succeeds (Fig. 2D,E), support the idea first proposed by Hallermann et al. that there is a tradeoff between energy minimization and the reliability of AP propagation 1 . In their study, successful propagation of APs into collaterals of the main axon depended on the inactivation kinetics of Na + channels, while in ours, it depended on the geometrical relationship between the diameters of the branches before and after the branch point. Additionally, a previous study 52 showed that the location of somata relative to the dendritic and axonal trees affects signal attenuation and hence metabolic costs. This result can be reinvestigated through direct energy calculations based on the energy function.
We believe this analytical approach will provide new insight into the metabolic costs of brain signaling, which is considered to be energetically expensive compared with other metabolic processes, although it is remarkably efficient compared with computers. On the one hand, the brain accounts for 20% of the body's total resting oxygen consumption [53][54][55] despite accounting for only 2% of body mass. This limited energy supply imposes serious metabolic constraints on brain function and evolution, including constraints on the density and size of neurons 7,56,57 and on functional connectivity 58 , network activity 59 , and coding strategies [60][61][62][63] . On the other hand, the human brain's power consumption is approximately 20 W 64 , which is merely 0.00001% of the power of today's most advanced supercomputers 65 . Finally, accurate estimation of the metabolic cost will improve the interpretation of functional magnetic resonance imaging data [12][13][14] .

Methods
Cable Model of Hodgkin-Huxley-type Cortical Axon. To describe AP propagation along an axon, the cable equation that describes the flow of ion currents along the axon needs to be derived 66,67 . Figure 1A gives the equivalent circuit of the cable model. During AP propagation, the membrane potential, V(x, t), changes along the x-axis, , R a = 0.15 kΩ · cm is the axial intracellular resistivity, a is the axon radius, and i a (μ A) is the axial current. According to Kirchhoff 's current law, i a varies along the x-axis, where i m is the membrane-crossing current and i stim is the AP-stimulating current. The membrane-crossing current i m consists of the capacitance and the ion currents,  where φ regulates the temperature dependence, with Q 10 = 2.3. Combining Equations 1-3, we obtain, a  a  m  stim   2   2 i.e., the cable equation, is the axial conductance and I stim = i stim /2π a (μ A/cm 2 ).
Energy Consumption of AP Propagation in the Cable Model. Inspired from the procedure applied to the classical HH model to compute the electrochemical energy involved in the dynamics of point model of single neuronal circuit 19 we plan to derive energy function of neuronal axon in the case of multi-compartment neuronal model as the method used in previous reference 19 (2 ) The role of the capacitance is to store energy from the battery supplier, not to consume energy, the total energy of the capacitors of the axon cable is The total consumed energy H g all by all dissipated items could be calculated by summation of battery energy, capacitance energy and the stimulus energy based on Eq. 6. To substitute Eqs. 7-9 into Eq. 6, we have ( ( , ) ) ( ( , ) ) ( ( , ) ) ( , ) ( , ) In summary, our approach is based on calculation of energy dissipation through conductances of the cable (including the axial conductance). In methods, we demonstrate the mathematical consistency of this with energy conservation and calculation of energy supplied by the ionic reversal potentials (and stimulus). During AP firing, the summation of external stimulus energy, the ion battery energy and the capacitor storing energy should be equal to the dissipated energy by conductors. Within the circuit, ion batteries supply energy, capacitors store energy while only conductors consume energy for action potential generation and initiation. Actually, computing the ATP consumption of the Na+ /K+ -ATPase pumps is a crude way to estimate the energy discharged from the batteries during AP firing. The explanation is as follows. The amount of energy stored in the ion batteries depends on the ion concentration gradients across the membrane. During AP firing, the membrane-crossing of the ions decreases the concentration gradients, which decreases the energy stored in the ion batteries. Conversely, after AP firing, the Na + -K + -ATPase recharges the ion batteries by restoring the ion concentration gradients to the equilibrium level. The ATP used in this restoring process is supposed to be the energy required to recharge the ion batteries after AP firing, which is numerically equal to the energy discharged from the batteries during AP firing.
Based on the explanation above, there are two main reasons for the inaccuracy of the sodium counting approach. First, it is not an accurate way to estimate the energy dissipated by the conductors. Second, even as a way to compute the energy discharged by the ion batteries, it only gives a coarse estimation by calculating the ATP molecules by the ratio of the number of pumped out sodium ions to needed ATP molecules, i.e. 3:1. According to our simulation, this coarse estimation is significantly lower than the direct calculation of the energy discharged from the batteries. Furthermore, whether the ATP calculating approach should be based on Na + -counting or K + -counting is controversial. Some believed that calculating the ATP molecules by the ratio of the number of pumped in potassium ions to needed ATP molecules, i.e. 2:1, is more accurate 5 . Actually, those two methods give different results.
From the biological perspective, we want to emphasize that, besides the sodium currents which have been widely recognized as making important contribution to the metabolic cost of AP firing, the potassium currents also contribute significantly to the cost. The potassium currents not only cause energy dissipation when flowing out of the membrane to restore the membrane potential, but also may be a part of the axial currents that result in the energy dissipation of propagating AP along the axon.