Mitochondrial morphology provides a mechanism for energy buffering at synapses

Mitochondria as the main energy suppliers of eukaryotic cells are highly dynamic organelles that fuse, divide and are transported along the cytoskeleton to ensure cellular energy homeostasis. While these processes are well established, substantial evidence indicates that the internal structure is also highly variable in dependence on metabolic conditions. However, a quantitative mechanistic understanding of how mitochondrial morphology affects energetic states is still elusive. To address this question, we here present an agent-based multiscale model that integrates three-dimensional morphologies from electron microscopy tomography with the molecular dynamics of the main ATP producing components. We apply our modeling approach to mitochondria at the synapse which is the largest energy consumer within the brain. Interestingly, comparing the spatiotemporal simulations with a corresponding space-independent approach, we find minor spatial effects when the system relaxes toward equilibrium but a qualitative difference in fluctuating environments. These results suggest that internal mitochondrial morphology is not only optimized for ATP production but also provides a mechanism for energy buffering and may represent a mechanism for cellular robustness.


S1.2 Local concentrations and gradients
For the scenario in which the ADP concentration is clamped at the OM, we estimated local ADP and ATP concentrations in the intracristal space (ICS), the inner boundary (IBM) and outer membrane space for the 3 different spatial arrangements (Figure 1). For this purpose, we counted with MCell the numbers of hits in an open region from which the concentrations can be estimated by assigning a small volume to this surface.
As explained in the main text, an initial drop in the ADP concentration is observed within the ICS when ANTs (adenine nucleotide translocators) are located in the CM as shown by the blue line in Figure 1C. The opposite holds for the ATP concentration which is persistently higher in the ICS due to the constant transport of ATP from the matrix to the cristae space (blue line in Figure 1G). We computed local concentration gradients between the inner boundary and outer membrane as well as between the cristae and the outer membrane ( Figure 2). To calculate these gradients, we measured the distance between the outer and inner membrane (≤ 0.02 µm 1 ) and between the middle of the mitochondrion to the outer membrane (∼ 0.15 µm), respectively, by the Blender add-on Measureit.
From the resulting data, we identified a strong gradient formed between the IBM and OM if ANTs are located in the IBM (grey line in Figure 2C) driving ATP export from the mitochondrion. In contrast, ATP import is observed in case ANTs are localized in the CM by an on average positive gradient (blue line in Figure 2C). In this case with ANTs localized in the CM, the gradient between the CM and the OM points to the exterior but is one order of magnitude smaller than the gradient between the IBM and OM. These observations underlay the results obtained in all different configurations.

S1.3 Additional variables
Here, we include all variables of the model that are not shown in the main text for both, the averaged traces of MCell simulations and the corresponding solution of the ODE system. Figure  We performed analogous simulations for the 3 different scenarios described in the main text but with a diffusion coefficient one order of magnitude smaller D = 1.5 · 10 −8 cm 2 s −1 . For this diffusion coefficient, substantial differences are observed in almost all cases (Figures 5,6 and 7). Remarkably, significant less ATP molecules (only ∼ 47) reached the cytosol after 10 ms when ANTs were located in the CM (blue line in Figure 6E). When ANTs were located in the IBM, ∼ 354 ATP molecules reached the exterior during the same period of time (grey line in Figure 6E). This diffusion limitation based effect was also found for faster diffusion as used in the main text but reduced diffusion is amplifying the differences between configurations.
Based on these simulations, we could estimate the rate at which ATP becomes available within the cytosol what is severely delayed under diffusion limitation conditions (Figure 7). For ANTs placed in the IBM, the ATP rate is 38 molecules per millisecond, and for ANTs in the CM the rate reduces to 11 molecules per millisecond whereas the rate for the spatially independent ordinary differential equation (ODE) yields 62 molecules per millisecond. If diffusion is reduced, the location of ANTs have a tremendous impact on the rate at which ATP can reach the cytosol and thus on the rate of energy supply at the synapse.

S1.5 Anomalous Diffusion
The extreme mitochondrial morphology has evoked a scientific discussion on diffusion properties within the organelle. Combining mathematical models and experimental data, the diffusive properties of proteins in the mitochondrial matrix have been measured (Partikian et al., 1998;Dieteren et al., 2011) and simulations of diffusion in obstructed objects suggested anomalous diffusion within organelles (Ölveczky and Verkman, 1998). To test this hypothesis, we used here our physiologically realistic reconstruction and simulated particle diffusion within the interior of the mitochondrion with and without the consideration of the internal structure (cristae membrane). From the resulting particle trajectories, we computed the mean-square displacement (MSD) of the molecules in both configurations. For normal diffusion, the MSD obeys Fick's law and equals 6 D t where D is the diffusion coefficient (set to D = 15 µm 2 s −1 ) and t represents time. In particular, this relation predicts a linear dependence of the MSD on time t.
We compared this theoretical prediction with the MSD measured in simulations with and without internal structure (Figure 8). For the comparison, we calculated the linear regressions of the measured MSD in both configurations and determined the slopes. The slope of the theoretical prediction for Fick's law in an open space is 15 µm 2 s −1 whereas we measured 44.7 µm 2 s −1 and 41.4 µm 2 s −1 for simulations without and with considering the internal structure, respectively. While the deviation on the large time scale is caused by the closed space of the mitochondrion, the deviation on the shorter time scale is influenced by the internal structure. Therefore, the deviation of the slope from the theoretical prediction is larger and rather instantaneous for the simulations considering the internal structure whereas the simulations only considering the closed space of the mitochondrion initially exhibits a rather similar slope and only deviates from the prediction on the longer time scale.

S2 Material and Methods
Here, we give the details on the modeling methods. For this purpose, we first introduce the molecular Markov chain models (MCMs) used to describe the ANT, the ATP synthase and the voltage dependent anion channel (VDAC) including explanations of model parameterization by experimental data. Subsequently, we describe the derivation of the corresponding ODE systems, which were used to validate the MCMs and to compare the spatiotemporal simulations with the space independent approach. Finally, we detail the estimation of the ATP production in the synapse and the time scale analysis.

ADP/ATP translocator (ANT)
The ANT model is based on published work (Metelkin et al., 2006) and it is composed of eleven kinetics states ( Figure 9) that represent the protein complex L in its different binding configurations. ATP and ADP correspond to T and D, respectively, and the index i refers to the matrix (interior) side, and o to IMS and ICS that represents together the outside space. Hence, L represents the free protein and YLX describes the molecular state with one X molecule bound from the matrix side and one Y molecule bound from the outside space. The kinetic rate constants are inferred from independent experimental data (in the following Section) and given in Table 2.
The dissociation constants of the ANT exhibit a complex dependence on the membrane potential (Klingenberg, 2008). Metelkin et al. (2006) developed a model introducing the dependence into the affinity of the adenylates and the rate constants. Our model is based on this previous work. Nevertheless more parameters needed to be estimated, and therefore others had to be slightly modified. The resulting parameters slightly deviate from those used by Metelkin et al. (2006) and are given in Table 2. The dissociation constant K To is changed from ∼ 400 µM to 500 µM in our simulations based on the parameter inference. Similarly, K Do changed from ∼ 51 µM to 25 µM . We also estimated the dissociation constants of ADP and ATP from the matrix side K Ti leading to a value of 6.25 mM and K Di of 10 mM, respectively.
Interestingly, these values for the matrix side are in the mM range, whereas the rates for the external side are in the µM range indicating a lower affinity of substrates on the matrix side. Similar values have been observed for the Michaelis-Menten constant Km for other mitochondrial carriers (Klingenberg, 2008).
From these dissociation constants, the backward and forward rate constants were estimated. The forward rate constants were set to be smaller than the diffusion limited rate (Milo and Phillips, 2015), and the backward rates were set to satisfy the dissociation constant. With these modifications, we were able to qualitatively reproduce the data from Kraemer and Klingenberg (1982); Duyckaerts et al. (1980) (data not shown). Although, the data was qualitatively well reproduced, the resulting turnover rate of 0.1 s −1 is two orders of magnitude smaller than the measured rate (Chinopoulos et al., 2009). This discrepancy might be caused by the limitations of the experimental data in Kraemer and Klingenberg (1982) that was generated with liposomes, and it is known that the lipid composition can affect the velocity of the translocator. To compensate for these limitations and to better calibrate our model, we included a third experiment (Chinopoulos et al., 2009) in which the transient appearance of ATP in the medium was measured for isolated synaptic mitochondria. From this data we were able to adjust the rate of ATP translocation (k p ) to match the measured turnover rate of 23 s −1 (Chinopoulos et al., 2009) for synaptic mitochondria. This led to a rate constant of 800 s −1 for the ATP exporting reaction DLT k p

TLD.
Similar values have also been measured in transient experiments with liposomes (Gropp et al., 1999). All the final parameters of the ANT model used for the presented simulations are summarized in Table 2.

ATP synthase
The ATP synthase model is composed of six states ( Figure 10A). Binding events from the matrix to the synthase E are indicated on the left and binding from the outside on the right, respectively. A clockwise cycle of the synthase starting from E −n corresponds to the binding of n protons from the outside and transport to the matrix followed by the binding of ADP and P i for the successive synthesis of ATP and unbinding of the n protons in the matrix side. The ATP generating step is indicated by the dashed arrow.
The states are assigned numbers ( Figure 10B), and the first order or pseudo first order rate constant of the specific transition, from state i to state j is k ij where we follow the same notation as in the original paper (Pietrobon and Caplan, 1985). State −3 E corresponds to state number 1, follow by state number 2 H 3 E, etc, following a counterclockwise rotation. All parameters values are given in Table 3.

ATP consuming reactions
To set the parameter values of the ATP consuming reactions, we used the calculations by Attwell and Laughlin (2001) for the energy demand of glutamatergic signaling. From these considerations, they estimated a demand of 11,000 ATP molecules per released vesicles, 12,000 ATP molecules for related Ca 2+ removal and ∼ 821 ATP molecules for endocytosis and exocytosis of vesicles. From the 11,000 ATP molecules needed for glutamate recycling, we only considered 1,333 (for the packing of glutamate in vesicles) since this is actually processed at the presynaptic terminal whereas the other processes are located within astrocytes. Thus, ∼ 14,000 ATP molecules are needed per vesicle released. For a firing frequency of 200 Hz and a release probability of 0.25, the total number of ATP molecules needed at the synapse equals 7 · 10 5 ATPs/s (0.25 × 200Hz × 1.4 · 10 4 ).
To consider ATP consumption within the synapse, we placed channel molecules at the synaptic membrane which consume ATP by the reaction ATP + channels k cha channels. For mimicking the action potential dependent ATP consumption, the reaction rate exhibits the form of two square pulses with a minimal and maximal value k basal and k max , respectively. The values of k max and the number of channels n cha were set to match an ATP consumption rate of 7 · 10 5 ATPs/s. The basal level of ATP consumption reflects additional housekeeping processes and is given in Table 4. For the synaptic simulations, we assumed a spontaneous firing rate of 5 Hz and vesicle release probability of 0.25 leading to a number of ∼ 1.8 · 10 4 ATP molecules per second needed at the synapse.

VDAC
In our model, we included VDACs to let ATP molecules exit the mitochondria into the cytosol. Since we were not aiming at capturing specific features of the channels, we implemented a rather simple Markov chain model that assumes that VDAC proteins interact with ATP and translocate it into the cytosol by the chemical reaction (VDAC + ATP mito VDAC + ATP cyto ). Based on experimental estimations of the VDAC density (De Pinto et al., 1987), we determined the number of channels n vdac as 1.1 · 10 4 for our specific mitochondrion model (Table 5). The rate constant of the reaction was set such that ATP export was not substantially delayed by the interaction with VDACs. Hence, in simulations with ATP molecules diffusing from a spherical region in the interior of the OM to the cytosol, the VDAC dependent export reduced the cytosolic ATP concentration to 75 % of the VDAC independent export within 10 milliseconds.

S2.2 Space independent ODE approach
For each of the model components, we derived an ODE system assuming mass action kinetics (Keener and Sneyd, 1998). The equations were integrated with PyDSTool (Clewley et al., 2007) using a Radau integrator.

ANT ODE model
The ODE system for the ANT model is given by Equations 1. The variables L, TL, LT, etc. represent the number of molecules in each state. One ANT can be found in one of the 11 states described in Figure 9.
For the 11 variables, we derived 10 differential equations and used the conservation of total number of proteins L tot = L + DL + DLD + DLD + TL + LT + TLT + TLT + DLT + TLD from which we deduced the number of molecules in the state DLT. The governing equations read: Here, D i and T i describe the number of ADP and ATP molecules in the matrix and D o and T o in the IMS and ICS building the outside space, respectively. T cyto describes the number of ATP molecules in the cytosol. The rate constants are labeled according to the reaction they drive, e.g. for reaction T o + L → TL, the corresponding rate constant is k f To and for the inverse reaction k b To . The reactions that translocate metabolites are: DLT k p TLD, exporting ATP to the outside and ADP to the inside, TLD k cp DLT bringing ATP inside and ADP outside, TLT k t TLT' transporting ATP from inside to outside and vice versa (this reaction matters for example when ATP within the media is labeled and we used this flux to determine the parameter k t from experimental data (Kraemer and Klingenberg, 1982)), and finally DLD k d DLD' that analogously exchanges to ADP from the inside and outside.
All parameters of the ANT model are summarized in Table 2.

ATP synthase ODE model
The ATP synthase model has six states represented in Figure 10A. In our simulations, n is set to 3 (i.e. 3 protons are needed to produce 1 molecule of ATP). As for the ANT, the total number of proteins is a conserved quantity (E tot ) leading to the conservation expression E −3 +EH 3 +H 3 E+H 3 ES+H 3 E * + −3 E = E tot . Thus, we derived a system of 5 ODEs (Equations 2) to describe the ATP synthase dynamic by and one more variable is deduced from the conservation expression. All parameter values are given in Table 3.

Metabolite ODE model
We also derived equations for the number of ATP and ADP molecules in the different compartments (Equations 3) given by For T cyto the rate depends on the VDAC model parameters (Table 5) as well as on ATP consuming reactions. Therefore, when ATP consuming reactions were considered, the following term was added to the equation −n cha k cha (t) · T o with k cha (t) being the rate of the reaction. This function has the form of two square pulses with t on being the activation time and t off the inactivation time for each pulse detailed in Table 4.
The rate constant of the bimolecular reactions have been normalized to have proper units to compute the number of particles, i.e. k = k NaVol , where N a is the Avogadro's number, and Vol represents the volume. Therefore, Vol is either the matrix volume or the volume of the outside space (ICS together with IMS).

S2.3 Estimation ATP production in synapses
Based on the turnover rate of ANTs and the number of translocators in the synaptic mitochondrion, we estimated the number of ATP molecules that can be exported to the cytosol per second ( Table 6). The product of the turnover rate times the number of ANTs gives the maximum number of ATP molecules that can reach the cytosol per second if the ANTs are not the limiting step. The total amount of ANT proteins in synaptic mitochondria has been estimated to 1.37 nmol/mg protein (Chinopoulos et al., 2009) assuming 1 nmol/mg protein ≈ 1.25 mM (Magnus and Keizer, 1997) leading to a concentration of 1.71 mM. The inner membrane has a volume of 0.021 µm 3 what correspond to ∼ 2 · 10 4 ANTs within the mitochondrion. This number of ANTs working at a turnover rate of 23 s −1 exports 4.6 · 10 5 ATPs/s. We further estimated the number of ATP produced by non-synaptic brain mitochondria (Table 6). This value critically depends on the volume of the mitochondria and, since no complete reconstructions are available in the literature, this is only a rough estimation based on an assumed cylindrical shaped mitochondria with a radius of 0.3 µm and length 1 µm (Perkins et al., 2001).

The number of ATP molecules produced by a synaptic mitochondrion coincides with the estimations of
Attwell and Laughlin (Attwell and Laughlin, 2001) explained above, in section ATP consuming reactions.
Moreover, a similar value of 6.02 · 10 5 ATP molecules per second has been also estimated for hair cell mitochondria (Babu et al., 2017). In Table 1 in the main text, we compared the values obtained with our simulations with estimations made by us and others.

S2.4 Time scale analysis
We calculate the spatial time scales for a diffusion coefficient of D = 1.5 · 10 −7 cm 2 s −1 by considering the characteristic length scales as the diameter of cristae junctions ∼ 25.5 nm 2 (L 1 ) and the distance from the center of the mitochondrion to the IBM ∼ 0.13 µm (L 2 ). The time required to transverse these distances is τ and is calculated as τ = L 2 /6D leading to A number of reactions occur in the ANT model. We calculated the corresponding time scales of some of the reactions to compare with the temporal scales of diffusion. The time scales of the forward reactions are denoted by τ f and for the backward reactions by τ b (Andrews and Arkin, 2006). For these calculations, we used the initial concentrations. As before L represents the free protein, D i denotes ADP molecules from the matrix side and D o ADP from the inner boundary membrane and cristae side (called outside) yielding: A similar analyzes has been done for the reactions of the ATP synthase model resulting in Different time scales are found for the reactions of the ANT and ATP synthase models that spread over a large range but most of them are between 1 · 10 −5 s and 1 · 10 −2 s ( Figure 11). We graphically represent the time scales of the reactions (blue dots) and the temporal scales of diffusion (black lines).