Computational Modeling of In Vitro Swelling of Mitochondria: A Biophysical Approach

Swelling of mitochondria plays an important role in the pathogenesis of human diseases by stimulating mitochondria-mediated cell death through apoptosis, necrosis, and autophagy. Changes in the permeability of the inner mitochondrial membrane (IMM) of ions and other substances induce an increase in the colloid osmotic pressure, leading to matrix swelling. Modeling of mitochondrial swelling is important for simulation and prediction of in vivo events in the cell during oxidative and energy stress. In the present study, we developed a computational model that describes the mechanism of mitochondrial swelling based on osmosis, the rigidity of the IMM, and dynamics of ionic/neutral species. The model describes a new biophysical approach to swelling dynamics, where osmotic pressure created in the matrix is compensated for by the rigidity of the IMM, i.e., osmotic pressure induces membrane deformation, which compensates for the osmotic pressure effect. Thus, the effect is linear and reversible at small membrane deformations, allowing the membrane to restore its normal form. On the other hand, the membrane rigidity drops to zero at large deformations, and the swelling becomes irreversible. As a result, an increased number of dysfunctional mitochondria can activate mitophagy and initiate cell death. Numerical modeling analysis produced results that reasonably describe the experimental data reported earlier.


Introduction
Mitochondria are widely accepted as the powerhouse and provide up to 90% of the ATP necessary for the cell. They are also important for the regulation of ionic homeostasis, cell growth, redox signaling, and cell death [1][2][3]. The electron transport chain (ETC) generates an electrochemical gradient across the inner mitochondrial membrane (IMM). This gradient, induced by the mitochondrial membrane potential (∆Ψ m ) and the H + gradient, is known as the proton-motive force that drives ATP synthesis. Mitochondria contain two membranes, the outer mitochondrial membrane (OMM) and the IMM, with an intermembrane space (IMS) between them. The OMM is a semi-permeable membrane allowing the passage of species with a molecular weight up to 6 kDa. The IMM is impermeable under normal aerobic conditions and transports selected ions and solutes through specific channels and exchangers.
An increase in mitochondrial Ca 2+ and reactive oxygen species (ROS) in response to oxidative and energetic stresses stimulates mitochondria-mediated cell death. Also, proteolysis due to protein oxidation can alter the structural and functional integrity of the IMM. The main mechanism of Calcium is a key inducer of both PTP-dependent and PTP-independent swelling of mitochondria ( Figure 1). Cytoplasmic Ca 2+ is rapidly taken up by mitochondria due to the proton-motive force generated by the ETC [23]. The mitochondrial calcium uniporter (MCU) is the main channel that transports Ca 2+ to the matrix when PTPs are closed. In fact, the MCU is a low-affinity and high-capacity channel that mediates Ca 2+ influx in a ∆Ψ m -dependent manner [24]. Mitochondrial Ca 2+ efflux occurs in three ways: (i) the mitochondrial Na + /Ca 2+ exchanger, the main contributor, most relevant in muscle mitochondria [25]; (ii) the H + /Ca 2+ antiporter, an Na-independent pathway [20,26,27]; and (iii) opening of the PTPs. The PTP opening occurs due to Ca 2+ overload, inducing swelling of mitochondria caused by the influx of water and ions through the open pores [25,28].
In addition to Ca 2+ , H + also plays an important regulatory role in PTP induction and, thus, mitochondrial swelling. Likely, H + regulates the sensitivity of the PTP towards Ca 2+ ; low pH inhibits pore opening, whereas pH above 7.0 stimulates permeability transition [26,29]. High H + levels have been shown to inhibit the PTP induction through cyclophilin D [30] and F-ATP synthase [31]. Therefore, pharmacological agents that reduce pH and create temporary acidosis are able to prevent PTP opening and exert beneficial effects in pathological conditions associated with mitochondrial Ca 2+ overload (e.g., ischemia-reperfusion) [32,33].
Modeling of mitochondrial swelling is important for the simulation and prediction of in vivo events under physiological and pathological conditions. The lack of knowledge on the molecular identity of the PTP complex obscures the understanding of the physical and chemical mechanisms of pore opening and mitochondrial swelling. Previous studies proposed several kinetic approaches to developing an optimized model of mitochondrial swelling. However, they do not provide an adequate mechanism underlying the transition from a reversible to an irreversible swelling state and should be improved by including additional factors. In particular, previous models did not take into consideration ∆Ψ m , the major parameter that regulates mitochondrial PTP induction, swelling, and cell death [34][35][36]. Simultaneous analysis of swelling and the dynamics of ∆Ψ m loss are important for understanding the relationship between changes in the matrix volume, ion flux, and mitochondrial bioenergetics. Also, previous models used the kinetic limit for the analysis of the transport of ions and solutes. As a result, diffusion was assumed to be faster than the respective kinetic rates. Therefore, an improved model should include diffusion-limited transport of Ca 2+ from the cytosol to the matrix to generalize the modeling approaches.
We have earlier reported [37] a simple model to interpret the experimental data on Ca 2+ -induced mitochondrial swelling. In the current study, we further develop modeling analysis including the transport dynamics of different ions and species across the IMM, and their effects on matrix metabolism, with the objective to describe mitochondrial swelling based on detailed physical and chemical characteristics. Our model describes the mechanism of transition from reversible to irreversible swelling in mitochondria. We developed a FORTRAN code to carry out numerical simulations, with varying values of the model parameters. Our modeling approach is based on the experimental data reported earlier by different authors. It includes the mechanism of matrix swelling, where the IMM rigidity is described by the rigidity tensor, its components dependent on the IMM deformation scale. Simultaneously, the osmotic pressure in the matrix is compensated for by the IMM deformation. Based on earlier suggestions [38], our model couples the PTP opening dynamics with the matrix pH. We conclude that the model is suitable for the analysis of the transition of mitochondria from the reversible to irreversible swelling. However, the applicability range of the model is quite limited, as it only takes into account Ca 2+ , K + , and H + to describe the irreversible swelling. In future, we plan to include a description of irreversible swelling in the more complex models with all of the relevant species, which currently disregard the possibility of irreversible swelling.

Model Description
There are two different approaches to the modeling of mitochondrial swelling dynamics [39]. The first approach introduces phenomenological kinetic parameters, which are fitted for better reproduction of the experimental kinetic results. The second approach develops a biophysical model based on the ionic and hydrodynamic fluxes across the IMM, and changes in ∆Ψ m and mitochondrial respiration rates [19,40].
The role of Ca 2+ in mitochondrial swelling has been assessed using a simple kinetic approach in rat liver mitochondria [41]. The kinetics of the permeability transition were analyzed within a model based on the assumption that the transition follows first-order kinetics and that the solute diffusion rate depends on the pore conformation. In addition, several theoretical studies concerning mitochondrial swelling dynamics were reported [39,[42][43][44][45][46]. Another study reported a model of mitochondrial ion transport and its application to the analysis of different modes of Ca 2+ uptake by mitochondria [47]. However, this model failed to reproduce the swelling and the role of PTP opening in swelling. Likewise, several systematic studies reported modeling of Ca 2+ transport across the IMM, and also PTP induction [42][43][44][45]. Their approach is based on a kinetic step that may be presented as Ca 2+ out ↔ Ca 2+ in , describing Ca 2+ ion exchange between the outside and inside of the IMM.

A Basic Model of Mitochondrial Dynamics
A simple scheme of the proposed model of Ca 2+ -induced mitochondrial swelling is shown in Figure 2. This model is an extended version of our earlier approach [37]. The Ca 2+ concentration affects ∆Ψ m = Ψ m, out -∆Ψ m, in , the proton concentration in the matrix C H, in , the PTP opening dynamics, described by the P PTP, op parameter, matrix concentration of the respiration activator (A), matrix osmotic pressure P os and matrix swelling. Mitochondrial swelling is induced by ionic/neutral species transport in/out of the matrix, which creates osmotic pressure inside the matrix. The osmotic pressure is compensated for by the IMM deformation induced by swelling. The present modeling approach aims to reproduce the IMM swelling dynamics, which is included in the transport model for the ionic and neutral species going in or out of the matrix, along with PTP opening dynamics. We used the Goldman equation to describe ionic species transport in and out of the matrix [47], which may be presented as follows: where J i is the flux of i-th ionic species, J i in and J i out are the i-th ionic species fluxes coming in and going out of the matrix, k B is the Boltzmann constant, T is the absolute temperature, z i is the relative ionic charge, |e| is the absolute electron charge, C in i and C out i are the concentrations of the i-th ionic species inside and outside the matrix, respectively, and p i is the reduced permeability of the IMM to the i-th ionic species. The transport of neutral species may also be described by Equation (1), substituting ∆Ψ m → 0, i.e., where p i (0) is the diffusion coefficient of the respective species in the uniport channels. We assumed p i (0) = 0 for the uniport of the ionic species. Our model uses both Equations (1) and (2). We used published values of the model parameters in our numerical analysis, with the respective data and references listed in Table 1.

Ca 2+ Transport across the IMM
Calcium influx through the IMM may be described by Equation (1), taking into account its background cytosolic concentration C 0,cb = 0.5 µM [48]. The effective reduced permeability of Ca 2+ uniport across the IMM is listed in Table 1. Our model includes the transport of Ca 2+ only. The IMM is impermeable to Ca 2+ , K + , Na + , H + and other ions; specific channels and exchangers control their flux and concentrations in the mitochondrial matrix. Specifically, mitochondrial K + balance is controlled by ATP-dependent (mitochondrial K ATP channel) and Ca 2+ -dependent channels responsible for influx, and by K + /H + exchanger responsible for removal of the excess matrix K + . Sodium balance is controlled by Na + /Ca 2+ (influx) and Na + /H + (efflux) exchangers [20]. As we already stated, the present model does not include Na + influx/efflux mechanisms. Therefore, the respective processes were omitted. However, we included the K + /H + exchange mechanism because K + transport at low [Ca 2+ ] (low-conductance PTP) induces oscillations of the mitochondrial volume [53,54]. This effect resulting from K + transport across the IMM is appropriately described by modified Equation (1), which may be presented as follows: where the permeability for the K + and H + ion transport through the membrane depends linearly on H + and K + concentrations:  Table 1.

PTP Transport of Ionic and Neutral Species
As we already mentioned, PTP opening also controls transport of different ions and neutral species to and from IMM. Factors controlling PTP opening include changes in [Ca 2+ ] and ∆Ψ m , and matrix pH dynamics [48]. Hence, for simplicity, we assume that Ca 2+ affects the PTP opening dynamics only indirectly, by inducing changes in ∆Ψ m or matrix pH. Presently, we follow the modeling approach proposed earlier [38], adding the model component required to describe the mitochondrial swelling dynamics. According to this approach, we assume that ∆Ψ m has no direct influence on the PTP opening dynamics in the autoinduced release of mitochondrial Ca 2+ . It is known that PTPs are completely open at pH ≥ 7.3-7.5, and completely closed at pH ≤ 7.0 [48]. Based on these facts, we define the probability of PTP opening as follows: This expression correctly produces P PTP = 0 at pH < 7.0, and P PTP = 1 at pH > 7.5. Figure 3 shows the probability dependence on pH at different values of α and n. On the other hand, the PTP opening rate and the PTP closing rate may be presented as follows: where k i is the respective rate constants (min −1 ), n m PTP is the number of PTP on the IMM of a single mitochondrion, n PTP,op is the number of open PTP, and n PTP,cl is the number of closed PTP. Using these notations, the equation for the PTP opening dynamics may be presented as follows: The latter relationship describes the dynamics of PTP opening. However, it does not include explicitly any effects of the matrix pH; these effects should then be included into the k i parameters. This makes the analysis quite complex. To simplify the calculations, we use Equation (3) instead of Equation (5). The matrix pH is now a time-dependent function, which we introduce below. Ion fluxes through the PTP may also be described by the modified Goldman Equation (1), and presented as follows: In case of H + transport, both C in H+ and C out H+ are variable parameters. The permeability coefficient for H + coming through PTP is unknown. We assume that its value is close to that for Ca 2+ , with both values listed in Table 1, where p H+,PTP = p Ca2+ , PTP /2 due to the difference in ionic charges [55]. We also assume that the value of K + permeability for PTP transport is the same as that for H + . Thus, we included H + , K + , and Ca 2+ transport through PTP into the current model. In our calculations, we used the conservation relationships: for K + and Ca 2+ , and used a constant external hydrogen ion concentration outside of the matrix C out H+ ( Table 2). We also used the same initial [H + ] in the matrix in all of the simulation runs and recorded its evolution in time. In the case of depolarized IMM, we may represent the ion transport through the PTP by the relationship similar to Equation (2): where p i,PTP (0) is the effective diffusion coefficient of the respective species through PTP, unknown for live cells. However, we neglect the thermal diffusion of ions through PTP in a depolarized IMM, Equation (7), as the ion release, in this case, is dominated by the hydrodynamic flux induced by mitochondrial rigidity.

Weak Acid Dissociation
Since the probability of PTP opening is directly dependent on the matrix [H + ] (C H+ ), we include the dissociation equilibrium for a weak acid AH. We use the latter as a simplified model for the second dissociation of phosphoric acid, falling into the pH interval of interest with pK a = 7.20. In this case, we obtain: where K AH is the equilibrium constant of the acid and C AH its initial concentration. Assuming that there are other sources of H + producing its background concentration C bg H+ , Equation (8) becomes: Next, we substitute concentrations by activities in Equations (1)- (8): where the activity coefficient γ I is determined using the Debye-Hückel theory [56]: where I is the ionic force, given by: We calculated activity coefficients using Equation (11). The mechanism involving electrodiffusion of anions and H + through PTP may be extended to the transport of other weak organic acids to/from the matrix through PTP. The respective kinetic model for a weak acid AH is described by the following reactions: The dissociation constants of weak acids are known quite well [56]. We assume that the dissociation equilibriums, described by Reactions (1) and (4), are much faster than diffusion, described by Reactions (3) and (4). We use the same values of the equilibrium constants inside and outside the matrix. Taking into account the relations (6), and (9) to (12), we described AH transport through PTP. Our model runs use the permeability value for the PTP transport of A − reported earlier [38] (see Table 1).

Effect of Respiration (ETC Activity) on H + Generation
Taking into account that ∆Ψ m directly affects the respiration and depends on the ETC activity [46,47], we used W A , the generation rate of the respiration activator A, in our model in the following form [55]: where k A [min −1 mM −1 ] is the effective respiration rate constant, δ(∆Ψ m ) is a variable parameter and a A,0 is the activator activity at zero time. The respiration rate is constant at ∆Ψ m = 200 mV [38].
This was achieved in the model by requiring that W A = 0.99k A a A,0 at ∆Ψ m = 200 mV, which results in δ(∆Ψ m ) = 43.42 mV, the value used in all model calculations. The respiration activator flux may be presented by the modified Equation (7): with the p A,PTP value listed in Table 1.
Mitochondrial respiration induces the formation of H + in the IMS due to its ejection from the matrix by ETC complexes I, III, and IV. Our model does not include the detailed mechanism of H + transport into the IMS. Therefore, it controls [H + ] by its generation via an irreversible process: with the S − product not participating in any other reactions. We determined the H + transport using Equation (6) and now, discuss the implementation of the effects of H + on ∆Ψ m . We calculated the rate of ETC-induced H + generation: where k A,H+ is the H + formation rate constant in the respiration cycle.

Model Implementation
The dynamics of ∆Ψ m is determined by the matrix ion fluxes, calculated as explained above. Thus, the system of equations describing the evolution of the model in time may be written as: where F is the Faraday constant, C is the IMM electric capacity, a i is the activity of the respective ion, S m the mitochondrial surface area, and B = 3 × 10 5 [38]. Here we assume that the evolution of P PTP , Equation (3) and that of W A , Equation (13), occurs much faster than the processes presented in Equations (16)- (20). We calculated the IMM electric capacity using the relationship C = S m h , where h = 10 nm is the IMM thickness, used as a constant value in our analysis.

Mitochondrial Swelling Dynamics
We complement our model with the mitochondrial swelling induced by osmosis. Osmotic pressure results from the concentration differences of various species inside and outside the matrix, written as follows [51,57,58]: where N is Avogadro's number. The higher osmotic pressure inside the matrix induces mitochondrial swelling, while the IMM deformation compensates for osmotic pressure up to a certain limit, maintaining equilibrium: where ∆P IMM is created by the IMM deformation. In the analysis of mitochondrial deformation, we assume that osmotic water transport in/out of the matrix is much faster than the rates of processes described by Equations (16)- (20). Our model calculations also assume ellipsoid shape for the mitochondrial matrix, with the IMM rigidity described by the second-order tensor: g xx g xy g xz g yx g yy g yz g zx g zy g zz with the matrix shape described by a prolate ellipsoid where the rigidity tensor may be diagonalized in the reference system shown in Figure 2: Using the rigidity tensor, we express the pressure created by the IMM deformation: where and S m (t) = 2π(a + ∆x) 2 The matrix volume changes and the resulting changes in the concentrations/activities are given by: where ∆t is the time step used in the numerical integration of Equations (16)- (20). To take into account the irreversibility of the mitochondrial swelling at large IMM deformations, the tensor components vanish at high deformations: g xx = g yy = g 9 = g 00 1 − β 0 ∆r n 1 1+β 0 ∆r n 1 r = x or y g zz = g zz,0 1 − β z ∆z n 1 1+β z z n 1 , (29) where g 00 , g zz,0 , β 0 , β z , and n 1 are adjustable parameters.
We also need to add to Equations (16)- (20) the terms describing the effluxes of the modeled species through PTP due to excess matrix pressure. According to Poiseuille's law, the volume of homogeneous fluid escaping through a cylindrical tube is given by [44]: where r is the effective PTP radius, η is the viscosity of the matrix fluid, and l is the IMM thickness. We define the effective PTP radius as: where r 0 ≤ 1.5 nm [26]. Thus, the hydrodynamic efflux of i-th species is given by: Thus, our model calculations used Equations (16)- (20) with added hydrodynamic flux terms according to Equation (30). Assuming spherically symmetric mitochondria (g xx = g yy = g zz ), the equations describing mitochondrial swelling dynamics become much simpler. The volume change rate will be given by: where η is the coefficient defining the rate of water transfer inside the IMM, and κ the coefficient describing the water transfer outside the IMM due to rigidity, vanishing at high deformations: Here, κ 0 is the linear rigidity coefficient, and γ and m the parameters characterizing membrane properties. At small volume changes, Equation (32) becomes: Its solution is still quite complex, because ∆P os is a time-dependent function, including ion transport and dilution dynamics created by volume dynamics. Using a linear approximation to Equation (33), the latter equation may be solved using Equation (22), with the result given by: Assuming κ −1 is much shorter than the characteristic times of other processes included in the model, Equation (34) takes the form: The earlier approach presented to mitochondrial swelling dynamics [46] is based on Equation (35). This approach, however, has limited utility, as swelling involves dilution of various components; additionally, Equation (35) is only applicable to spherically symmetric systems with homogeneous swelling. This equation also disregards the irreversibility of large deformations, introduced here by Equation (32).
All of the above considerations apply to the analysis of mitochondrial volume dynamics as a function of different system parameters. Our main objective was to gain an understanding of irreversible mitochondrial swelling, testing the new modeling approach. We are considering a set of equivalent mitochondria, with their number density given by n mit (cm −3 ), an additional model parameter.

Numerical Experiments
Since mitochondrial dimensions vary significantly between biological species and cell types, our model analysis specifically targets mitochondria of the adult rat cardiomyocytes. We implemented the numerical model in a homemade FORTRAN code, with the input parameter values listed in Table 2.
Our numerical experiments used constant values for the α, n, C AH,0 , δ(∆Ψ m ), C A,0 , k A , and k A,H+ parameters, selected ( Table 2) according to previous reports [38,43,46,51,55,[57][58][59][60]. We estimated the values of g 00 and g zz,0 for spherical mitochondria based on their average mass and expected oscillation frequency [61]. We used the experimental criterion of irreversible swelling observed upon 50% increase in volume; this translates into ∆r = 0.144a and ∆z = 0.144c, and the rigidity tensor components vanishing at a = 0.5 µm and c = 1 µm [61]. Thus, and using Equation (28), we selected β 0 , β z , and n 1 values as listed in Table 2. We also introduced an outside [H + ] = C out H+ , kept constant during the model run. We also kept constant C in H+,0 in , the matrix [H + ], in each of the numerical experiments. The values of these two parameters were selected using previous reports [61], with the respective data listed in Table 2. We selected K + and Ca 2+ concentrations, which we varied independently, to cover their entire natural range. We solved the system of equations defining our model using a homemade FORTRAN code, and the parameters of Table 2, producing the following output parameters: with PTP opening probability P PTP and ζ = V m (t)/V m (t = 0) − 1.
All of the calculations used ∆Ψ m (t = 0) = 200 mV, C in Ca2+ (t = 0) = 0.5 µM, and C in K+ (t = 0) = 0. In all of the numerical experiments, we kept the ratio of the cell volume to the initial mitochondrial volume constant. Usually an adult rat cardiomyocyte contains~5000 mitochondria, with the cell to mitochondrial volume ratio of ca. 2.50-2.86 [62]. Thus, the free cell volume per each mitochondrion is (2.50 − 2.86)V m (t = 0). Our numerical analysis additionally used the values of various constants as reported earlier [52,63,64] (Table 1).

Results and Discussion
According to Table 2, we varied only two variable parameters: the initial K + out and Ca 2+ out concentrations. We made a total of 100 independent numerical experiments with the two concentrations of K + and Ca 2+ that were varied independently. However, we limited the analysis to only two concentrations of K + out (0.01 and 10 mM) and 10 different Ca 2+ out concentrations (Table 2), producing a total of the 20 most interesting numerical experiments. We calculated the dynamics of the H + in , K + in , and Ca     [49,50], and described theoretically [38] in previous studies. Experimental measurements of [K + ] oscillations and a theoretical analysis of oscillations were reported earlier, with oscillation periods of 20 s and 93 s [38,46].
These results are quite similar to those produced by our model, in particular for pH oscillations. Notably, our model differs from previous studies since it includes the swelling of mitochondria. Figure 4b,d,f plot matrix [Ca 2+ ], [K + ] and matrix pH for different initial cytosolic [Ca 2+ ]. An interesting result is observed at t < 600 s: the matrix [Ca 2+ ] never exceeds 100 µM, while the initial cytosolic [Ca 2+ ] is 500 µM. The model generates correct distribution between the "in" and "out" Ca 2+ at low cytosolic [Ca 2+ ] < 200 µM, in agreement with earlier reports [24,43]. However, strong deviations from the expected equilibrium distribution between "in" and "out" Ca 2+ are produced at high cytosolic Ca 2+ concentrations, between 300 and 500 µM. These deviations may be attributed to the PTP efflux of Ca 2+ from the matrix, generated by the excess pressure within the matrix due to IMM deformation. This effect has not been reported before. All of the data plotted in Figure 4b,d,f also demonstrate that irreversible processes arise at initial cytosolic [Ca 2+ ] in the 200-300 µM range. This result agrees with acceptable accuracy with the earlier reported experimental data [37,41].  Figure 5a plots [38,41]. Figure 5b plots the dynamics of ∆Ψ m , with weak oscillations at the initial cytosolic [Ca 2+ ] of 1.0 µM. This result agrees with earlier model data [38], reporting weak oscillations of ∆Ψ m in similar conditions. However, previous studies demonstrated different results reporting deep modulation of ∆Ψ m , from 0 to 175 mV [46]. We will address these discrepancies in a future publication by the inclusion of the swelling dynamics. It should be noted that our swelling model is quite different, whereas the frequency of oscillation reported in our study is in agreement with the earlier reported model [46]. In conclusion, ∆Ψ m remains reversible for the initial cytosolic [Ca 2+ ] in 25-200 µM range; however, it becomes irreversible at cytosolic [Ca 2+ ] exceeding 300 µm, and completely depolarized at 500 µM. These results agree with the data reported earlier [26,43,65]. Figure 5c plots the dynamics of the matrix volume as V 0 − 1 for clarity. Note that the volume oscillations appear already at the lowest cytosolic [Ca 2+ ]. The same results were obtained earlier both in experimental measurements [54] and theoretical modeling analysis [46]. Once again, we found an acceptable agreement between the oscillation periods obtained here and reported earlier [46]. Apparently, reversible volume changes arise at initial cytosolic [Ca 2+ ] < 300 µM, followed by irreversible swelling at high (>300 µM) [Ca 2+ ]. The maximum volume increase of 50% of the initial mitochondrial volume was observed, as expected on the basis of the model parameterization. This volume increase corresponds to vanishing rigidity constant, Equation (28) and the data of Table 2, and, therefore, to the irreversible deformation and mitochondrial death.
We report significant differences between the results obtained at  Figure 6a shows that all of the matrix concentrations reveal weak oscillation dynamics, with periods similar to those obtained for the initial cytosolic [K + ] = 0.1 µM. However, the modulation depth at [K + ] = 10 µM is smaller than that obtained at the initial cytosolic [K + ] = 0.1 µM. The same is apparent for the ∆Ψ m and mitochondrial volume dynamics. Thus, the changes in the initial cytosolic [K + ] produce only minor changes in the dynamics of the key system parameters.
Thus, our numerical model produces results that are different from those obtained in previous studies using a more sophisticated model [46]. The strongest discrepancies in the modulation depth of different model parameters appear at low [Ca 2+ ]. For instance, we report ∆Ψ m modulation depth of 5-7% of the maximum value, while the previous study [46] reports 100% modulation depth, although still in the reversible mode. Apparently, even complete IMM depolarization does not affect the system dynamics in the respective model. This issue will be addressed in a future study, with the presently proposed swelling mechanism included in a more complex model. The swelling mechanism proposed in the previous report only considers the osmotic pressure, disregarding the interaction between the IMM rigidity and IMM deformation.

Conclusions
In the present study, we implemented a simplified model of mitochondrial swelling that includes the contribution of osmosis and the IMM rigidity. This is a novel biophysical approach to mitochondrial swelling dynamics, which we will incorporate into more complex and sophisticated models. The numerical experiments in the framework of the developed model produced reasonable results, which describe the earlier reported experimental data with an acceptable accuracy.

Limitations of the Study
The modeling approach explored in the present study is limited by the simplified treatment of the biophysical and chemical processes, including transport of only three ionic species and a simplified respiration mechanism. However, we intentionally focused on the development of a more complex mechanism of mitochondrial swelling, and its usage in the analysis of the more sophisticated modeling approaches reported previously.