Mitochondrial morphology governs ATP production rate

Energy conversion is at the core of cellular function. Finding efficient strategies to transduce energy has driven the evolution of cells and organisms. In eukaryotic cells, respiration occurs in specialized organelles known as mitochondria. There, metabolic substrates are converted to adenosine triphosphate (ATP), the energy currency of the cell (Alberts et al., 2015). Mitochondrial architecture provides an efficient framework for cellular respiration. It is composed of two distinct membranes: the outer membrane (OM) and the inner membrane (IM). The IM is used to generate an electrochemical proton gradient through a series of electron carriers known as the electron-transport chain (ETC; Fig. 1). The energy stored in this gradient drives the synthesis of ATP from adenosine diphosphate (ADP) and inorganic phosphate (Pi) through a multiprotein complex, ATP synthase (Yoshida et al., 2001). The outer membrane plays a central role in controlling the flow of molecules from and to the interior of mitochondria. Although it forms a barrier to large molecules, it contains a large number of proteins, voltage-dependent anion channels (VDAC; Colombini, 2004), that allow the passage of small metabolites and ions. In contrast, the IM is impermeable to molecules; carrier proteins and ion channels are required for molecules to cross this barrier. One of the most relevant carriers is the adenine nucleotide translocator (ANT), which transports ATP and ADP to and from the matrix, the internal region encapsulated by the IM (Fig. 1). Overall, this is a membrane-based process that relies on the lipid bilayer of mitochondria to generate the energy required for ATP synthesis.

The importance of ATP production by mitochondria means that reliable mathematical formulations are necessary to build a quantitative framework. To that effect, there have been many mathematical models developed to describe the processes involved in ATP production (Pietrobon and Caplan, 1985; Magnus and Keizer, 1997; Cortassa et al., 2003; Nguyen et al., 2007; Bertram et al., 2006; Beard, 2005; Heiske et al., 2017; Garcia et al., 2019). Pietrobon and Caplan (1985) developed a proton pump model that can be interpreted as a respiration-driven proton pump or as a model of ATP synthase (Magnus and Keizer, 1997). Later on, Magnus and Keizer (1997) assembled different components and built a steady-state model of ATP production in mitochondria, which several other models are based on (Cortassa et al., 2003; Bertram et al., 2006; Nguyen et al., 2007). However, they did not explicitly show how the thermodynamic constraints were imposed on the system¹. Some of these models also have been built to be consistent with thermodynamic laws (Beard, 2005; Heiske et al., 2017), and in most cases include more components as the activity of the TCA cycle (Cortassa et al., 2003; Nguyen et al., 2007) enzymes or the activity of other components of the ETC (Beard, 2005; Heiske et al., 2017). Almost all previously developed models are implemented with ordinary differential equations (ODEs), disregarding the spatial organization of mitochondria. Here and in our previous work (Garcia et al., 2019), our efforts were focused on developing a minimal reaction-based mitochondrion model so that we can implement it to perform spatial simulations in MCell (Stiles et al., 1996; Kerr et al., 2008; Husar et al., 2022 Preprint) or non-spatial simulations. In this work, we ensure the model is thermodynamically consistent, and for that we imposed detailed balance constraints to the reaction rate constants.

With the model developed here, we answered the following questions: under physiological conditions, what limits the amount of ATP that reaches the cytosol? Is it the rate of ATP production or is it the rate of ATP export from the mitochondria to the cytosol? To answer these questions, we developed a thermodynamically consistent model of ATP production in mitochondria—i.e., a physically plausible model. Our criteria for thermodynamic consistency was to calculate thermodynamic forces (Hill, 1989) for each cycle in the model and to equate them to the ratio of the forward-to-backward reaction rate constants. In this manner, we verified that detailed balance is satisfied for the conditions considered. Furthermore, we calculated the dependence of the reaction rate constants on the membrane potential, pH, ATP, ADP, and Pi concentrations. We conducted simulations of this model using ODEs; although simpler, this approach is informed with detailed structural calculations of surface areas and volumes from a realistic 3-D reconstruction (Garcia et al., 2019). The outcomes reveal that, under physiological conditions, ATP production in the matrix is the limiting factor. Lastly, we used the model to perform spatial simulations using nine realistic 3-D reconstructions (Mendelsohn et al., 2021). Our results link the ATP production rate to structural features of the mitochondrial geometries.

The model presented here is based on previous work (Pietrobon and Caplan, 1985; Magnus and Keizer, 1997; Metelkin et al., 2006; Garcia et al., 2019). Our efforts were focused on assembling a thermodynamically consistent model (Hill, 1989) and explicitly highlighting the dependence of the reaction rate constants on the membrane potential, pH, ATP, ADP, and Pi concentration.

The ATP synthase model is grounded on the work of Pietrobon and Caplan (1985). The kinetic model consists of a membrane protein—represented as E—that transports protons (H⁺) from the intermembrane space (IMS) to the matrix. The protein can be in six states, represented in Fig. 2 A. Each state is associated with a number from 1 to 6 (Fig. 2 B). A clockwise cycle starting in state 6 (E⁻³) represents the binding of three protons from the IMS (transition 6 → 5 in Fig. 2 B), transport of the protons to the matrix (state 5 → 4), binding of ADP and Pi (state 4 → 3), and subsequent synthesis of ATP (state 3 → 2), followed by unbinding of the protons in the matrix (state 2 → 1). We considered a constant proton concentration inside the IMS and the matrix, and the ATP and ADP concentration as variables. The list of the reactions is given in Table 1, and the model parameters are given in Tables 2. The number of ATP synthases used in the model was determined employing the density estimated at 3,070 µm⁻² (Acehan et al., 2011). Given the surface area of high curvature of the IM in the reconstructions (Mendelsohn et al., 2021), we estimated the number of ATP synthases as 267 for mitochondrial reconstruction #1 (Mito 1).

We next explain how some kinetic parameters depend on the membrane potential, the proton concentration, and the phosphate concentration.

The ANT model is based on the work of Metelkin et al. (2006). It is composed of 11 states and 20 bidirectional transitions. The kinetic diagram is shown in Fig. 2 F. The free protein (L) can bind molecules from the matrix (on the right) or IMS side (on the left), forming a triple molecular state (YLX), where Y and X represent ATP and ADP, respectively. State YLX can transition to state XLY, exporting X to the IMS and importing Y, following unbinding of the molecules. The productive transition that exports ATP is DLT → TLD, and the reversal reaction imports ATP to the matrix and exports ADP to the IMS. Two additional states (DLD’ and TLT’) were added to track futile translocation (exchange of ADP for ADP and ATP for ATP). Starting from previously fitted flux parameters (Metelkin et al., 2006) for ANTs extracted from heart mitochondria (Krämer and Klingenberg, 1982), we first estimated the parameters of the dynamic implementation of the model. With this set of parameters, we qualitatively reproduced the independent data from published work (Krämer and Klingenberg, 1982; Duyckaerts et al., 1980). The list of the reactions is given in Table 3. The rate constants of the model are presented in Tables 4 and 5, and the qualitative reproduction of the experimental data is given in Fig. S1 and Fig. S2. Additional information about the experiments and comparisons with other models can be found in the supplemental text at the end of the PDF.

This equation summarizes the imposing restrictions on the rate constants.

As before, the rate constants of the bimolecular reactions have been normalized to have proper units to compute the number of particles, i.e., $k=k′NA*V$⁠, where k′ represents the constants given in Table 4, N_A is Avogadro’s number, and V represents the volume.

As for the ATP synthase model, the total number of ANTs is a conserved quantity $(Ltot=L+DL+DLD+DLD′+TL+LT+TLT+TLT′+DLT+TLD)$⁠. L_tot is the total number of ANTs considered in our simulations; for the mitochondrial reconstruction #1, this quantity is 16,471 (this parameter is #ANTs₁ in Table 4).

The main mechanism for metabolites to cross the OM is through VDACs. To consider the ATP export from the IMS into the cytosol, we included VDACs. We implemented a rather simple model assuming VDAC proteins interact with ATP molecules and translocate them to the cytosol by the reaction $VDAC+ATPIMS⇌VDAC+ATPcyto$⁠. The list of parameter values employed is given in Table 6.

The total concentration of ATP and ADP can be distributed in different compounds or states like ATP⁴⁻, ADP³⁻, ATPMg²⁻, etc. ATP and ADP can react with different cations, or be bound or ionized. The final distributions can be estimated by coefficients representing the fraction of unbound ATP in the matrix or the other compartments, and these proportions are determined mainly by the level of Mg²⁺. For our model, mitochondrial ADP³⁻ and ATP⁴⁻ concentrations were estimated analogously to a previous work (Magnus and Keizer, 1997). These estimations are based on experimental data (Corkey et al., 1986). The approximate distributions are [ADP]_m,free = 0.8 [ADP]_m, [ATP]_m,free = [ATP]_m, [ATP⁴⁻] = 0.05 [ATP]_free, and [ADP⁻³] = 0.45 [ADP]_free.

The initial concentrations of ATP and ADP in the IMS and cytosol were set to 6.5 and 0.1 mM, respectively (Mörikofer-Zwez and Walter, 1989; Cunningham et al., 1986). The initial concentrations of ATP and ADP in the matrix were set to 13 and 2 mM (15 mM in total), respectively. The total concentration of adenylates in the matrix (ATP + ADP) of liver mitochondria is in the range of 2–27 mM (Rulfs and Aprille, 1982; Nosek et al., 1990; Austin and Aprille, 1984), in the same range for heart mitochondria (Livingston et al., 1996), and also ∼15 mM for pancreatic cells (Magnus and Keizer, 1997). Moreover, the ratio of ATP to ADP in the matrix has been estimated to be about 10 (Nicholls, 2013). These results are consistent with our assumption of an initial concentration of 13 mM of ATP and 2 mM of ADP in the matrix. In some simulations, these concentrations were kept constant. The concentrations were converted to the number of molecules using the volumes of the different compartments (see Table 7).

The system of ODEs has a total of 19 equations and 41 parameters, including chemical rate constants, volumes, and number of proteins. The volumes of the different compartments are listed in Table 7. We used reconstruction #1 for all the non-spatial simulations.

Eqs. 29, 37, and 40 were integrated using a Radau stiff integrator in PyDSTool (version 1.19.2; Clewley et al., 2007).

To perform spatiotemporal simulations in realistic 3-D reconstructions, we implemented the set of reactions in MCell version 4.0 (Kerr et al., 2008; Stiles et al., 1996; Husar et al., 2022,Preprint), an agent-based reaction–diffusion simulator. The model is built using the set of reactions described in Tables 1 and 3. Previous research showed that ATP synthases accumulate at the rim of the cristae membrane in lamellar cristae and along the length of tubular cristae (Strauss et al., 2008; Davies et al., 2012). To emulate this organization in our reconstructions, we selected the areas of high curvature in the cristae membrane (Mendelsohn et al., 2021) and we randomly distributed them in these areas. The areas of high curvature were identified as vertices of the cristae membrane with the first principal curvature higher than 70 µm⁻¹. The location of ANTs in mitochondria has not been definitively resolved, and experimental evidence shows they may form complexes with ATP synthases (Ko et al., 2003) located in the cristae membrane (CM; Wittig and Schägger, 2009; Vogel et al., 2006). Therefore, we randomly distributed ANTs in the areas of high curvature as well. VDACs were homogeneously distributed on the OM. ATP and ADP molecules were homogeneously distributed in the respective compartments. The diffusion coefficient considered for the free forms of ATP and ADP is 1.5 × 10⁻⁷ cm² s⁻¹.

In Figs. S1 and S2, we compared experimental results from Krämer and Klingenberg (1982) and Duyckaerts et al. (1980) to the results obtained with the ANT translocator model. In Fig. S3, we show ATP fluxes through ANTs generated with the model developed by Bohnensack (1982). In Fig. S4, we present spatial simulation performed with the same conditions considered in the main text, but kept the number of ATP synthases constant in all reconstructions. In Fig. S5, we compare the result with the previously developed version of the mitochondrial model (Garcia et al., 2019) with the current version. Supplemental text appears at the end of the PDF and provides additional information about the experiments and comparisons to other models.

Before studying the dynamics of the assembled model, we first analyzed the behavior of its components separately and compared them with previously developed models. In these first simulations, we assumed steady-state conditions and kept the ADP and ATP concentrations constant in the different compartments. We then conducted dynamic simulations with the ODE representation to understand which is the limiting step in ATP production. Finally, we performed spatial simulations with nine realistic 3-D geometries and related the ATP production rate to the morphology of the organelles.

First, we analyzed the dependence of the ATP flux on the ADP concentration, membrane potential, and pH for the ATP synthase model. The driving forces of cycle c (Fig. 2 C) are the membrane potential, the pH, and the ratio of the ADP and Pi concentration to the ATP concentration relative to the equilibrium constant (Eq. 26). For these simulations, the ΔpH was set to −0.4 (Tischler et al., 1977) and the total concentration of adenine nucleotides in the matrix was set to 15 mM (Austin and Aprille, 1984; Rulfs and Aprille, 1982). In Fig. 3 A, we show the ATP net flux through ATP synthases for different concentrations of ADP in the matrix. A positive flux is defined in the counterclockwise direction, favoring ATP hydrolysis. We found that at small concentrations of ADP, the flux is positive, meaning ATP is hydrolyzed. Increasing the concentration of ADP in the matrix, the ATP flux becomes negative, indicating that ATP synthesis is favored. The flux saturates at concentrations of ADP above 1 mM. Moreover, the driving forces of the cycle also depend on the membrane potential (Eq. 26). For a membrane potential above 160 mV, the ATP flux is negative for large concentrations of ADP, but for membrane potential below those values, ATP is hydrolyzed. Similar results were obtained for the ATP flux derived by Magnus and Keizer (1997) and later modified by Cortassa et al. (2003) (for comparison see Fig. S6 A in the supplementary material in Cortassa et al. [2003]; the difference in the sign of the fluxes is due to different sign conventions).

Next, we studied the dependence of the ATP flux on the pH difference between the matrix and the IMS (Fig. 3 B). Here, we fixed the pH in the matrix at 7.6 and varied the pH in the IMS. When the pH difference between the matrix and IMS is large, a smaller membrane potential is needed to observe a negative ATP flux, i.e., to synthesize ATP. This is because the proton motive force—the membrane potential minus the difference in pH (Eq. 21)—drives the protons to the matrix.

The ATP flux also saturates at higher membrane potentials, and in particular is almost saturated at 180 mV—the potential considered for the rest of this work. Similar flux dependencies were observed in the past (for comparison see Fig. S6 B in the supplementary material in Cortassa et al. [2003]), although their maximal flux is about 10 mM s⁻¹ while ours is about 0.9 mM s⁻¹. A similar dependence on the membrane potential was also observed by Nguyen et al. (2007) (Fig. 2 D in their paper), and their maximum flux is 1.125 mM s⁻¹, very close to the range in our work.

The driving force of the ANT translocation is the membrane potential and the ratio of ATP_m/ADP_m in the matrix and ATP_i/ADP_i in the IMS (Eq. 33). Therefore, we studied the dependence of the net ADP_i − ATP_m flux through ANTs at steady concentrations of ATP and ADP for different ratios. First, we calculated the net ADP_i − ATP_m flux through ANTs for different ratios ATP_m/ADP_m in the matrix (Fig. 4 A). Keeping the cytosolic ratio fixed at 65, we varied the ATP_m/ADP_m ratio. The driving force of the cycle (Eq. 33) increases for larger ATP_m/ADP_m ratios, increasing the net ADP_i − ATP_m flux (Fig. 4 A). The net ADP_i − ATP_m flux saturates for ratios ATP_m/ADP_m closer to the ratio of ATP_i/ADP_i (Eq. 33).

Next, we varied the ATP_i/ADP_i ratio in the IMS and kept the matrix ratio constant at 6.5. Increasing the ATP_i/ADP_i ratio decreases the driving force (Eq. 33) at a constant ATP_m/ADP_m ratio. Therefore, the net ADP_i − ATP_m flux decreases for larger ratios of ATP_i/ADP_i (Fig. 4 B). Moreover, the maximum flux is set by the ATP_m/ADP_m ratio. Similar flux dependencies were found by Magnus and Keizer (1997) with the model developed by Bohnensack (1982), see Fig. S3.

From this analysis, we identified that for 1 mM ADP_m, 180 mV membrane potential, and −0.4 ΔpH, the ATP flux is almost at its maximum absolute value (Fig. 3, A and B). We used these values for the rest of this work. Likewise, the initial concentration ratios considered in the rest of this work are ATP_m/ADP_m = 6.5 and ATP_i/ADP_i = 65, from which the ADP_i − ATP_m flux through ANTs is almost at its maximal value.

The amount of ATP in the cytosol depends on the production rate and its export through ANTs and VDACs. Since experimental results suggest that the ANT translocation rate constants are in the order of 100s per second (Gropp et al., 1999), we explored the impact of fast export through ANTs on the final ATP concentration in the cytosol. We held the concentration of ADP in the IMS constant and changed the parameters of the ANT model. Specifically, we changed the value of the rate constants k_p, k_cp, k_t, and k_d multiplying them by a factor of r, such that detailed balance is still satisfied. We measured the number of ADP and ATP molecules in the matrix (Fig. 5, A and B) and the number of molecules of ATP in the IMS and cytosol (Fig. 5, C and D). Faster translocation allows more ADP to accumulate in the matrix (Fig. 5 A) and therefore less ATP (Fig. 5 B), while ATP molecules reach the IMS and the cytosol faster (Fig. 5, C and D). We further calculated the percentage of the ATP increase relative to the slowest import (Fig. 5 E). The fastest ANTs (r = 10) exports 10% more ATP molecules than the slowest ones. Moreover, the amount of ATP starts to saturate for the fastest ANTs (r = 10). This means even ANTs that are faster cannot overcome the ATP production rate. Thus, the rate-limiting step is the ATP production rate, not the export.

In addition to the translocation rate of the ANTs, the ratio of ANTs to the ATP synthases can influence the export of ATP molecules into the cytosol. Therefore, we next asked how this ratio affects the availability of ATP in the cytosol for the fastest ANTs. We asked which step limits the number of ATP molecules in the cytosol under these conditions—is it the synthesis of ATP or its translocation through ANTs? As in the previous simulations, the number of ADP molecules in the IMS was constant, assuming unlimited resources. We measured the number of ADP and ATP molecules in the matrix, IMS, and cytosol over time (Fig. 6, A–D).

The amount of ATP in the cytosol after 1 s starts to saturate when the number of ANTs is 20 times the number of ATP synthases (Fig. 6 E). Thus, our model predicts that having more than 20 times ANTs to the number of ATP synthases will not have any benefit. If this ratio is bigger than 20, the limiting step is the production of ATP.

Similarly, VDACs can limit the amount of ATP in the cytosol. Therefore, we analyzed how the export of ATP through VDACs limits their amount in the cytosol. Considering similar conditions as before the constant concentration of ADP in the IMS, we changed the number of VDACs in the OM and measured the number of ADP and ATP in the different compartments (Fig. 7, A–D), along with the percentage increment of ATP in the cytosol after 1 s (Fig. 7 F). Having less than ∼500 VDACs in the OM significantly limits the amount of ATP available in the cytosol. Otherwise, the production of ATP is the limiting step.

Given the density of VDACs in the OM (De Pinto et al., 1987) and the surface area of the OM measured (Mendelsohn et al., 2021; see Materials and methods), the number of VDACs in reconstruction #1 is 6,340. This number is well above the limiting threshold of 500. Therefore, in the physiological conditions we considered, the production of ATP is the limiting step.

To illustrate the utility of our model, we conducted spatial simulations with nine realistic 3-D reconstructions (Mendelsohn et al., 2021) using the agent-based reaction–diffusion simulator MCell, and we related the ATP production rate to morphological features of mitochondria. In Fig. 8 A, we show the CM of the nine reconstructions with the total number of ATP synthases considered for the simulations (in red). In Fig. 8 B, we show the results of the ODE system and averaged traces of the spatial simulations (over 10 different realizations) of the number of ADP and ATP in the different compartments for the elongated mitochondrial reconstruction #2 for the same conditions considered as before.

Using tools from differential geometry, we recently analyzed the structural features of these reconstructions, and based on this analysis we group them into two categories: globular organelles and elongated ones (Mendelsohn et al., 2021; Fig. 8 A). Organelles in the globular group have larger surface area formed by vertices with high curvature. Large curved-surface areas correspond to more number of ATP synthases in these reconstructions.

We performed spatial simulations for all reconstructions and measured the ATP production rate in the cytosol for each reconstruction. For that, we performed a linear regression to single traces of the number of ATP molecules in the cytosol over time. Averaged rates are shown in Fig. 8 C. As expected from our analysis, the ATP production rate in the cytosol is larger for globular organelles, which have larger areas of high curvature. For comparison, we also performed simulation with the same number of ATP synthases in all reconstructions. The results can be found in Fig. S4.

Given the indispensable nature of mitochondria in cellular energy production, having accurate models of mitochondrial dynamics is critical for understanding their impact on cellular function.

In this work, we focused on developing a thermodynamically consistent model of ATP production in mitochondria to ensure the model is physically plausible. In doing so, we derived the relationships for the dependence of kinetic parameters on the membrane potential, proton, and substrate concentration.

We explored different parameter regimes to understand in which conditions ATP production or its export are the limiting steps in making ATP available in the cytosol. Under the conditions considered, we found two regimes: at a low ratio of ANTs to ATP synthases, the export of ATP is limited by the ANTs. At ratios higher than 20:1, the ATP production becomes the limiting step. Similarly, the number of VDACs does not limit the export of ATP, unless there are very few VDACs present (in the order of 100s, for the ratio of ANTs/ATP synthases considered).

Finally, we performed spatial simulations on nine mitochondrial reconstructions derived from electron microscopy data. According to our parameter exploration, the number of ANTs and VDACs in these nine reconstructions are such that the synthesis of ATP is the limiting step on the ATP production rate in the cytosol and not its export. We further measured the rate for the same conditions and relate this quantity to the structural features of the organelles. We had mapped features of the OM to features of the IM recently in Mendelsohn et al. (2021) with tools from differential geometry, and showed that globular organelles—characterized by smaller aspect ratio and smaller average first principal curvature of the OM—have larger curved CM surface areas. These organelles also have other features of high energy capacity like a large crista shape factor and a high number of crista junctions per OM surface area. Others have shown that mitochondria at high-energy-demanding locations in outer hair cells of young mice (Perkins et al., 2020), as well as mitochondria from high-spiking interneurons (Cserép et al., 2018), show a large crista shape factor and a high number of crista junctions per OM surface area. These features are also characteristic of our globular set, which suggests a high level of metabolic activity for such mitochondria. Here, using our mathematical model, we found that globular organelles generate more ATP in the cytosol than the elongated morphologies (Fig. 8 C).

The maximum rate measured with our spatial model is roughly 30 ATPs/ms (Fig. 8 C), and this corresponds to 170 µM s⁻¹ (considering the volume of the cytosol). Recent measurements in myocytes (Wescott et al., 2019) estimated the ATP production rate in the range of 200–600 µM s⁻¹ depending on the calcium concentration. The organelles considered in this manuscript are axonal and presynaptic reconstructions and are small (in the range of 0.037 to 0.012 µm³) compared with organelles in myocytes (in the order of 0.5 µm³; Laguens, 1971). This could explain our smaller ATP production rate in the cytosol.

We acknowledge that our modeling approach has a number of limitations. First, the membrane potential in our simulations is constant, and we emulate a constant state of ATP production. Second, the ATP synthase flux considered in our model saturates sigmoidally at large membrane potentials (Fig. 3 B), as it has been shown for Propionigenium modestum bacteria and Escherichia coli (Kaim and Dimroth, 1999), and is more sensitive to changes in the membrane potential than to ΔpH. However, recent experiments on heart and skeletal mitochondria (Wescott et al., 2019) show no saturation of the ATP synthase flux when increasing the membrane potential. Our model does not reproduce this behavior. Its voltage dependence is a consequence of the detailed dynamic of the ATP synthase model (Pietrobon and Caplan, 1985), which is drawn from physical assumptions about the rate constants, the membrane potential, and the pH concentration. Taking different assumptions or expanding the number of states could perhaps give a different flux membrane potential dependence. Third, the Mg²⁺-bound forms of ADP and ATP are the main substrates of ATP synthase in mitochondria. Since our goal was to generate a spatial model of a mitochondrion and Mg²⁺ is found in large concentrations in the matrix, we considered the free form of ADP as the main substrate for ATP production. This is an approximation due to the limitations of using an agent-based model. However, this approximation has been taken in many mitochondrial models in the past (Magnus and Keizer, 1997; Bertram et al., 2006; Cortassa et al., 2003), and is accepted in the community (Gnaiger, 2020).

Different approaches can be taken to model the spatial organization of mitochondria: from continuum approximations (Leung et al., 2021; Afzal et al., 2021) to discrete descriptions as we considered in our spatial simulations (with an agent-based reaction-diffusion simulator). Well-mixed continuum approximations based on concentrations like those in Leung et al. (2021) and Afzal et al. (2021) break down at low particle numbers. A discrete particle-based approach, as we considered in our spatial simulations with an agent-based reaction-diffusion simulator, is more accurate in this scenario (Bartol et al., 2015). Continuum approaches could be more challenging to integrate into complex geometries. Since we want to explore the impact of the complex mitochondrial morphology on ATP production and include the specialized location of ATP synthases at the rim of the cristae, we opted for a discrete description in our spatial simulations.

These different approaches can lead to similar results. For instance, in the past, we observed gradients of ATP with our discrete approach (Garcia et al., 2019), and so did Afzal et al. (2021), with a continuum description. Our earlier version of the model (Garcia et al., 2019) is based on the same published data (Pietrobon and Caplan, 1985; Metelkin et al., 2006; Magnus and Keizer, 1997). The major improvement of the current version is that the parameters were constrained by thermodynamic considerations (Hill, 1989). The results, however, are consistent (in Fig. S5 we compare for reconstruction #1 both versions of the model).

Finally, our results have implications for our understanding of disease states. For instance, near-total loss of ATP synthase dimerization caused by a neurodegenerative mutation generated profound disturbances of mitochondrial crista ultrastructure in neurons derived from a Leigh syndrome patient (Siegmund et al., 2018). The main architectural feature different from control cells is an increment in the angle of curvature at the tips of the cristae. This feature relates to our curvature measurement of the CM, which directly relates to the total number of ATP synthases. Altered mitochondrial structures have also been shown in aged animals (Perkins et al., 2020; Glavis-Bloom et al., 2023). Further investigation using the methods developed here will allow a deeper understanding of potentially impaired mitochondrial mechanisms in disease states.

All the material associated with this work is available in the repository: https://github.com/RangamaniLabUCSD/spatial_mito_model.

David A. Eisner served as editor.

This work was supported by the Air Force Office of Scientific Research (AFOSR) Multidisciplinary University Research Initiative (MURI) grant FA9550-18-1-0051 (to P. Rangamani) and Office of Naval Research N00014-20-1-2469 (to P. Rangamani).

Author contributions: Conceptualization: G.C. Garcia and P. Rangamani; data curation: G.C. Garcia; formal analysis: G.C. Garcia and P. Rangamani; funding acquisition: T.J. Sejnowski and P. Rangamani; investigation: G.C. Garcia and P. Rangamani; methodology: G.C. Garcia, T.M. Bartol, and P. Rangamani; project administration: G.C. Garcia and P. Rangamani; resources: G.C. Garcia and P. Rangamani; software: G.C. Garcia and T.M. Bartol; supervision: G.C. Garcia, T.M. Bartol, T.J. Sejnowski, and P. Rangamani; validation: G.C. Garcia, K. Gupta, and T.M. Bartol; visualization: G.C. Garcia; writing - original draft: G.C. Garcia and P. Rangamani; writing - review & editing: G.C. Garcia, T.M. Bartol, T.J. Sejnowski, and P. Rangamani.

All the material associated with this work is available in the repository: https://github.com/RangamaniLabUCSD/spatial_mito_model.