In silico simulation of reversible and irreversible swelling of mitochondria: the role of membrane rigidity.

Mitochondria have been widely accepted as the main source of ATP in the cell. The inner mitochondrial membrane (IMM) is important for the maintenance of ATP production and other functions of mitochondria. The electron transport chain (ETC) generates an electrochemical gradient of protons known as the proton-motive force across the IMM and thus produces the mitochondrial membrane potential that is critical to ATP synthesis. One of the main factors regulating the structural and functional integrity of the IMM is the changes in the matrix volume. Mild (reversible) swelling regulates mitochondrial metabolism and function; however, excessive (irreversible) swelling causes mitochondrial dysfunction and cell death. The central mechanism of mitochondrial swelling includes the opening of non-selective channels known as permeability transition pores (PTPs) in the IMM by high mitochondrial Ca2+ and reactive oxygen species (ROS). The mechanisms of reversible and irreversible mitochondrial swelling and transition between these two states are still unknown. The present study elucidates an upgraded biophysical model of reversible and irreversible mitochondrial swelling dynamics. The model provides a description of the PTP regulation dynamics using an additional differential equation. The rigidity tensor was used in numerical simulations of the mitochondrial parameter dynamics with different initial conditions defined by Ca2+ concentration in the sarco/endoplasmic reticulum. We were able to estimate the values of the IMM rigidity tensor components by fitting the model to the previously reported experimental data. Overall, the model provides a better description of the reversible and irreversible mitochondrial swelling dynamics.


Introduction
Mitochondria are double-membrane organelles, consisting of functionally and structurally distinct outer mitochondrial membrane (OMM) and inner mitochondrial membrane (IMM). Mitochondria contain two compartments: a small intermembrane space (IMS) between the two membranes, and the interior site (matrix) of the organelle enclosed by the IMM. The mitochondrial matrix is responsible for the mitochondrial metabolic pathways including, among others, tricarboxylic acid cycle, heme synthesis, urea cycle, and fatty acid oxidation. Changes in the matrix volume of mitochondria play a crucial role in mitochondrial life and death (Beavis et al., 1985;Bernardi, 1999;Halestrap et al., 1986;Lehninger, 1959). Mild mitochondrial swelling (MS) under physiologic conditions regulates metabolism and function of mitochondria, whereas excessive swelling causes mitochondrial dysfunction. Swelling of mitochondria is induced by external (cytoplasmic) and internal (matrix) factors that affect ion gradients across the IMM and increase osmosis in the matrix (Bernardi, 1999;Szabo and Zoratti, 2014). Among ions, Ca 2+ plays the central role in the swelling of mitochondria where strong Ca 2+ exchange system across the IMM maintains Ca 2+ homeostasis in the matrix under physiological conditions. Ca 2+ ions enter the matrix through several IMM channels including the mitochondrial calcium uniporter (MCU), rapid mode channel (RaM), and mitochondrial ryanodine receptor 1 (mRyR) (Giorgi et al., 2018;O'Rourke, 2007). The MCU is a highly selective channel that mediates Ca 2+ influx in a ΔΨ m -dependent manner (Kirichok et al., 2004). The main mechanisms responsible for Ca 2+ release (efflux) are the mitochondrial Na + /Ca 2+ exchanger, the Ca 2+ /H + exchanger, and the permeability transition pore (PTP) (Griffiths, 2009;Hoppe, 2010;Szabo and Zoratti, 2014). Also, Ca 2+ levels in the mitochondrial matrix may be regulated by physical and functional interactions between mitochondria and sarco/endoplasmic reticulum, a large reservoir of Ca 2+ (Lopez-Crisosto et al., 2017).
Increased Ca 2+ in the matrix accompanied by ATP depletion, high phosphate ions (P i ) and accumulation of ROS causes an abrupt increase in the permeability of the IMM to solutes and ions through the nonselective channels, namely PTPs. The PTP-induced swelling of mitochondria is driven by the colloid-osmotic pressure, which increases in the presence of non-diffusible matrix proteins and a high osmotic gradient across the IMM. In addition to Ca 2+ , the PTP opening depends on pH, IMM potential (ΔΨ m ), and the redox state of the cell; low pH inhibits whereas low ΔΨ m and high ROS stimulate the PTP opening (Bernardi and Di Lisa, 2015;Halestrap et al., 1998;Javadov et al., 2009). Ca 2+ competes with a variety of negative modulators such as ADP, cyclosporin A (cyclophilin D inhibitor), divalent cations (e.g., Mg 2+ ), and protons (Bernardi et al., 1992;Javadov et al., 2017). Mitochondrial PTP opening can occur in low-conductance (physiologic) and high-conductance (pathologic) modes [Reviewed in (Brenner and Moulin, 2012;Kwong and Molkentin, 2015)]. The PTP flickering in the low-conductance mode increases permeability to solutes ≤300 Da, mostly ions, and induces negligible matrix swelling . The volume-induced increases in the IMM permeability over the physiological range stimulate ETC activity, ATP production, fatty acid oxidation, and other metabolic processes (Halestrap, 1994;Scalettar et al., 1991). Also, low-conductance PTPs initiate mitochondrial depolarization spikes and generate Ca 2+ waves from one mitochondrion to another caused by Ca 2+ -induced Ca 2+ release (Ichas et al., 1997). In contrast, high-conductance PTPs allow unrestricted bi-directional movements of water and solutes ≤1500 Da across the IMM accompanied by increased ROS production, ΔΨ m loss, uncoupled oxidative phosphorylation, and ATP hydrolysis (Bernardi and Di Lisa, 2015;Halestrap et al., 1998;Javadov et al., 2009). Excessive matrix swelling causes rupture of the OMM and release of pro-apoptotic proteins (e.g., cytochrome c) into the cytoplasm (Petronilli et al., 2001), although the relative contribution of apoptosis and necrosis to cellular death depends on the ATP level in the cell.
Despite many studies, the mechanisms underlying the MS dynamics remain poorly known. One of the main challenges in understanding the mechanisms of MS is the lack of knowledge of the molecular identity of the PTP complex. Although voltage-dependent anion channel (porin), adenine nucleotide translocase, and phosphate carrier were initially proposed as the PTP core components, subsequent genetic studies demonstrated that pore opening still occurs in the absence of these proteins, suggesting that they are not involved in the PTP structure (reviewed in (Bernardi and Di Lisa, 2015;Halestrap and Richardson, 2015;Javadov et al., 2009)). Recent studies implicate F 0 F 1 -ATP synthase as a structural component of the PTP (Carraro et al., 2019), although other studies challenge this conclusion (Carroll et al., 2019). Cyclophilin D has been shown as a major regulator of PTP and pharmacological inhibition of cyclophilin D reduces pore opening and MS (Halestrap and Richardson, 2015;Javadov et al., 2017). Lack of knowledge on the molecular structure of the PTP complicates elucidation of the mechanisms of PTP-induced MS in vitro and in vivo.
Simulations of MS coupled with the bioenergetic status of mitochondria by different modeling approaches can provide a powerful alternative to exhaustive experiments . These models are particularly important for understanding the mechanisms of the transition of MS from a reversible to an irreversible state. Several in silico approaches have been used to develop biophysical and kinetic models for the MS dynamic analysis (Baranov et al., 2008;Bazil et al., 2010b;Massari, 1996;Pokhilko et al., 2006;Selivanov et al., 1998). However, most of these models describe MS as a linear process, which is unable to reproduce irreversible MS or the transition from reversible to irreversible MS. Notably, there are fundamental differences between biophysical and kinetic models: biophysical models describe continuous mitochondrial dynamics depending on the initial conditions, whereas kinetic models describe the transition dynamics between well-defined discrete mitochondrial states. Therefore, additional modeling tools are required for joining these two theoretical approaches together.
In the present study, we attempted to upgrade and further develop our recent biophysical model for the description of the MS dynamics Makarov et al., 2018). The dynamics of the PTP opening probability in the current study is described by a firstorder differential equation. This equation describes how the dynamics of the PTP opening probability depends on the matrix concentrations of Ca 2+ and H + . It illustrates the time evolution of the PTP opening status and thus provides a more correct approach in comparison with the earlier discussed modeling approaches. One of the main factors regulating the MS is the IMM rigidity, which may be described by the rigidity tensor. The IMM strain by osmotic pressure creates IMM mechanical stress, which compensates the osmotic pressure in equilibrium conditions. The values of the rigidity tensor components were obtained fitting the upgraded model to the earlier reported experimental data.
Using new values of the tensor components, we recalculated the time dependences of different system parameters, comparing the results to those reported earlier. The tensor components obtained in the present study significantly differ from our earlier estimates. The currently developed and verified tools may be further applied to the analysis of the irreversible MS dynamics described by more complex biophysical models.

Model description
We have recently developed a biophysical model that describes MS based on i) ΔΨ m dynamics, ii) ion transport through the PTP, iii) dynamics of the PTP opening probability and iv) nonlinear dynamics of MS . The model took into consideration a limited number of ionic species (H + , K + , and Ca 2+ ) as well as ion transport mechanisms in the IMM including the PTP opening, among others. To further improve the previous model, presently we included a first-order differential equation describing the dynamics of the PTP opening probability as a function of the Ca 2+ and H + concentrations in the matrix. We shall first briefly clarify the differences between our current model, and previous models.

Ion transport across the IMM by ion exchangers and PTP
The ion transport through IMM is described by the Goldman equation: where p i is the permeability coefficient for the i-th ionic species, which includes the |e|/k B factor , T is the absolute temperature, k B is the Boltzmann constant, |e| is the absolute value of the electron charge, C i in/out is the i-th ionic species concentration inside and outside of the matrix, respectively, and z i is the relative charge of the i-th ionic species. Considering that PTP opening dynamics is described by the opening probability T op , Eq. (1) for ion transport through PTP is modified as follows: (2) Both Eqs. (1) and (2) were used for the description of ion transport V.I. Makarov, et al. Mitochondrion 50 (2020) 71-81 in our present model. Thus, the dynamic equation for the i-th ionic species transport may be represented as follows: where S(t) is the time-dependent mitochondrial surface area and V(t) is its time-dependent volume. Note that the permeability coefficient p i included in Eqs.
(1)-(3) was calculated as p i = p' i /S(0), where the p' i was the initial integral permeability coefficient of the IMM published earlier . This representation of the permeability coefficients allowed describing the ion transport rates in the function of the IMM surface area, which was changed by swelling. Regarding the description for C H+ in in our present model, Eq.
(3) was complemented with the W res term, describing H + generation due to respiration, and took the form: where respiration is described by the transformation of the respiration activator A. Respiration also creates a weak acid, which dissociates and alters pH. The respiration activator dynamics is described by: where C A in is the respiration activator concentration, W ox and W rem are activator generation rate and its transformation rate into the weak acid AH, respectively. The dissociation of AH is given by AH ↔ A − + H + .
Since ΔΨ m is directly coupled with the respiration dynamics, our present model described the W ox generation rate as where k A is a phenomenological rate constant, C S,0 in is the initial substrate concentration and δ(ΔΨ m ) = 43.42 mV is the characteristic parameter .

ΔΨ m dynamics
In the simplest case of modeling analysis, when the surface area of the ion exchanger channels and PTP were not included, and only average values of permeability coefficients were considered, the dynamic equation for ΔΨ m may be represented as follows: where |e| is the absolute electron charge, z i is the relative charge of i-th ionic species and C m is the ΔΨ m .

Dynamics of PTP opening probability
We used the following equation to describe the PTP opening probability: where p(ΔΨ m ) was given by: where k PTP m = 3.2 × 10 −5 mM, Ψ 0 = 110 mV, Ψ 1 = 20 mV, Ψ 2 = 70 mV, k op = 2 min −1 , and k cl = 100 min −1 . The upgraded model was analyzed numerically, and the results of the analysis are discussed in Section 3.

Nonlinear mitochondrion dynamics
Our previous model of MS used a nonlinear description of swelling based on the osmotic pressure created by ionic disbalance between the matrix and cytosol . The osmotic pressure was represented by the following equation: where N is the Avogadro number. The osmotic pressure inside the matrix induces MS, while the IMM strain compensates for the osmotic pressure in equilibrium conditions up to a certain limit, maintaining equilibrium: where ΔP IMM is created by the IMM strain, which was described by the rigidity tensor g. For mitochondria approximated by an ellipsoid of revolution ( Fig. 1), the pressure created by the IMM strain was described as follows: where the rigidity tensor components were defined as follows: where Δr(t) = r(t) − a and Δz(t) = z(t) − c. The parameters g 00 , g zz,0 , β 0 , β z , and n 1 describe the IMM rigidity and the surface area: and the volume of the IMM: Irreversible MS would occur at large IMM strains when the regidity tensor components vanish. The equation for the MS dynamics should include water fluxes into and out from the matrix. The water flux into the matrix is proportional to the osmotic pressure, while the flux out V.I. Makarov, et al. Mitochondrion 50 (2020) 71-81 from the matrix (IMS) is proportional to the pressure created by the IMM mechanical stress. The permeability constant p w for water transport is the same in both directions. Thus, the equation for the mitochondrion volume depending on the transported water volume may be represented as follows: Eq. (15) produces two coupled differential equations describing the dynamics of Δr and Δz. The equation for Δr may be rewritten in a spherically symmetric system as follows: where r 0 is the initial radius, Δr′ is the time-dependent IMM radial strain, g′ 00 is the rigidity constant of the spherical mitochondrion, β′ 0 and n 1 are the parameters describing the IMM strain nonlinearity. Numerical simulations were carried out for both ellipsoidal and spherical mitochondria.

Results
Here we simulated the experimental data reported earlier and carried out detailed numeric simulations of the system parameter dynamics at different initial conditions, comparing these results with the earlier reported numerical simulations .

Experimental data simulation
To obtain the IMM rigidity tensor components, we fitted the experimental data reported earlier by another group (Holmuhamedov et al., 1999) with the results of the present study (Fig. 2). Note that the biophysical model used for fitting is very limited as regards the number of variables. As the fitting procedure was applied simultaneously to two conjugated parameters (Ca 2+ concentrations in the sarcoplasmic/endoplasmic reticulum and in the matrix volume), and conscious of the model limitations, we optimized the fitting as much as possible. However, the fitted curve in Fig. 2b still deviates quite strongly from the experimental data. We expect that the fitting will produce better results with a more comprehensive model that would follow a larger number of ionic species. However, the presently used simple model estimated several of the essential parameters with an acceptable accuracy; these results were used in further numerical experiments.
In our simple modeling approach mitochondria were approximated by a sphere or an axi-symmetric ellipsoid. However, the real mitochondrial structure is much more complex, and the size and shape of mitochondria change dynamically in vivo. It should be noted that modeling analysis of any real biological systems requires certain approximations that do not take into consideration all of the available factors and parameters. As a result, a theoretical model only partially reflects the real biological system (event), in particular, our model does not take into account the dynamics of a real mitochondrial structure. Therefore, it is very important to enhance the accuracy level of the used approximations. We carried out fitting of the experimental data using axi-symmetric ellipsoidal and spherically symmetric systems, and our parameter values describing mechanical IMM properties do not dramatically differ from the real parameters. Consequently, we may expect that in conditions, when MS is developed, the used approximation remains acceptable for the analysis of reversible and irreversible MS, and can predict mitochondrial behavior in different conditions with a certain degree of accuracy. It should be noted that significant progress in adjusting the modeling analysis to reality may be obtained by the model upgrade that would take into account an extended set of factors affecting the system dynamics. Further studies should be extended to more complete biophysical models, where a mitochondrion is represented by a non-axi-symmetric ellipsoid, and an extended set of ions and processes is included. This would permit a better description of the real mitochondrion geometry. Thus, our current work is the first step to developing more realistic modeling tools.

Numerical analysis of mitochondrial dynamics
Thus, in the present model, we used the new values of the rigidity Here, V m,0 and V m (t) are the initial and current values of the mitochondrial volume, respectively. Our main objective is to improve our modeling tools; therefore, we analyzed the MS dynamics only for the initial conditions listed above. The currently reported modeling tools may be used in conjunction with more comprehensive biophysical models that would include transport of a larger set of ionic and neutral species across the IMM. On the other hand, we were testing the results of the current model against those of our previous developments, to evaluate the currently achieved improvements. Such an analysis provided indications for future improvements of the tools describing the nonlinear MS dynamics.
Since the 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 1. The numerical calculations were performed for [Ca 2+ ] = 1.0 and 500 μM in the sarco/endoplasmic reticulum, using the same initial concentrations of all of the other species. In the next section, we present the results of the numeric simulations.
It should be noted that we did not analyze the model of MS for different mitochondrial subpopulations in the current study. Note that there are three spatially distinct mitochondrial subpopulations in cells: subsarcolemmal, intermyofibrillar, and perinuclear mitochondria. Subsarcolemmal and intermyofibrillar mitochondria are mostly discussed in cardiac and skeletal muscle cells. These subpopulations demonstrate distinct functional and structural characteristics (Hollander et al., 2014;Kuznetsov and Margreiter, 2009) including different sensitivities to Ca 2+ (Palmer et al., 1986), among others. In addition, sex and age differences in the sensitivity to Ca 2+ were observed in cardiac mitochondria (Arieli et al., 2004). More detailed and complex models for different mitochondrial subpopulations depending on sex and age could be developed by upgrading the current simple model that would precisely predict MS dynamics of mitochondria in spatially distinct subcellular compartments in male and female mitochondria at a given point of the lifespan.

Dynamics of MS
We simulated the dynamics of ΔΨ m , pH, ζ and T op for an ellipsoidal mitochondrion using the same initial geometrical parameter values, where a = 0.5 μm, c = 1 μm (see Eqs. (13)) and the effective ratio of the cell volume to the mitochondrion volume is 2.86 (Song et al., 2013). We compared these results with the previously reported data, which led us to the following conclusions (Fig. 3): i) the pH dynamics has not changed significantly; ii) the dynamics of the PTP opening probability has changed weakly for [Ca 2+ ] ext 0 = 500 μM, disappearing at longer times; iii) the ΔΨ m dynamics has changed weakly; the difference observed between the two calculations for [Ca 2+ ] ext 0 = 500 μM apparently originates in the ΔΨ mdependent parameter, p(ΔΨ m ), and iv) the matrix volume dynamics for low [Ca 2+ ] has changed significantly, with the difference once more attributable to the p(ΔΨ m ), whereas no significant changes in the matrix volume dynamics were apparent at higher [Ca 2+ ]. We, therefore, conclude that the upgraded model describes the matrix swelling better than the earlier reported model , as it provides a more precise description of the rigidity tensor component values and the PTP opening probability dynamics.
To better understand the results obtained using the upgraded model, we represented the differences between the current and our earlier reported data  in Fig. 4. Fig. 4a demonstrates such difference for the matrix pH at [Ca 2+ ] ext 0 of 1.0 and 500 μM. For both Ca 2+ concentrations, the difference has a long period of oscillations (~500 s), which may be assigned to a better description of the mitochondrion state dynamics using Eq. (8). Notably, the amplitude of the difference in both cases is around 0.6, which is sufficient to affect the PTP opening dynamics. The difference in the PTP opening probability shown in Fig. 4b demonstrates decaying oscillations. Interesting effects can be seen in the differences for ΔΨ m (Fig. 4c), where stable oscillations with the oscillation period of about 64 s are observed for [Ca 2+ ] = 1 μM. This value is practically the same as observed for the PTP opening dynamics. In this case, we can suggest that in both cases the observed effects may be apparently assigned to the same reasons for low Ca 2+ (ΔΨ m ) and high Ca 2+ (PTP opening). However, at [Ca 2+ ] = 500 μM, ΔΨ m difference shows low-frequency decaying oscillations with an oscillation period of about 170 s. The volume dynamics differences show stable oscillations with the oscillation period of about 110 s at low Ca 2+ (Fig. 4d), while at high Ca 2+ , the difference demonstrates a peak with the maximum at around 277 s and the width of about 210 s. All of the observed effects in the parameter difference dynamics result from the better description of the PTP opening probability. Note that these effects are even observable at low Ca 2+ .
The differences between the current and the previous data on the PTP opening probability dynamics are shown in Fig. 3b. A detailed representation of the PTP opening probability dynamics shown in this figure for [Ca 2+ ] 0 out = 500 μl is given in Fig. 5. Note the regular oscillations with the period of ca. 37.9 s, which were quite unexpected in conditions leading to irreversible MS. The initial oscillation amplitude is about 0.031, and it decays with time, with the time constant of the exponential decay function of ca. 287 s. This effect should be caused by Ca 2+ uptake by the mitochondria, which occurred with a similar characteristic time, i.e. while the matrix Ca 2+ was below the critical level, the system was still controlling the PTP opening dynamics. However, the respective oscillation amplitude was only 3%, and difficult to detect experimentally. Therefore, the existence of such oscillations could be tested in numerical experiments, using comprehensive biophysical models with nonlinear MS description. Provided these models were able to reproduce such oscillations, we would gain a better understanding of the role of different interactions, leading to better modeling of the mitochondrial dynamics in different conditions.  . This should be caused by the irreversible swelling that eventually reduces the mitochondrial volume, apparently due to a strong efflux of water from the matrix through the PTP. The data analyzed in this section allow us to conclude that the behavior of the simulated mitochondrial parameters (matrix pH, ΔΨ m , PTP opening probability, and matrix volume) as a function of [Ca 2+ ] 0 out is quite reasonable. The results presented in Fig. 6a demonstrate that pH is about 7.17 and the relative mitochondrial volume is about 0.17 at high [Ca 2+ ]. The maximum matrix pH is around 9.4 and drops to about 7.17 at an infinite time for 500 μM Ca 2+ in the sarco/endoplasmic reticulum (Fig. 3a). Thus, matrix pH asymptotically approaches the value of 7.17 on this time scale. The same is apparent for the data shown in Fig. 6d, where the relative mitochondrial volume asymptotically approaches 0.17. To test these conclusions, we calculated the relative volume dynamics in the extended time range of [0; 6000] s. Both dependencies presented in Fig. 7a and b show that the system arrives to an almost steady state t > 3000 s, with the asymptotic values of 7.17 and 0.17 for pH and the relative mitochondrial volume, respectively.

Discussion
Currently, the existing experimental approaches have a limited ability to produce precise in silico models that describe the mechanisms of the MS dynamics in vivo. Besides, existing experimental data on the MS were obtained using a wide spectrum of techniques, experimental conditions, and tissue sources. This compromises a correct interpretation of the results and the development of a general model of MS. The changes in the influx and efflux mechanisms across the IMM for a broad range of ions and solutes, along with the physicochemical parameters of mitochondria such as ΔΨ m , PTP flickering, membrane rigidity, matrix colloidal osmotic pressure, and others, all should be monitored simultaneously. The latter approach can track more accurately the dynamics of the MS and the transition of mitochondria from a reversible to an irreversible state. A quantitative biophysical model would be very helpful in the interpretation of the experimental data, with the model parameters evaluated by simulation of the experimental results.
The MS generates mechanical stress in the IMM determined by the IMM strain and membrane rigidity, in response to the high colloid-osmotic pressure in the matrix. Large IMM strains during severe (pathological) MS stimulate mechanical degradation, inducing mitochondrial and cellular death. A wide range of biophysical modeling approaches was applied in previous studies that described the MS as a reversible linear process (Baranov et al., 2008;Bazil et al., 2010b;Massari, 1996;Pokhilko et al., 2006;Selivanov et al., 1998). The earlier discussed where λ is the IMM permeability for water, S m is the IMM outer surface area, R is the gas constant, T is the absolute temperature, ρ H2O is water density, C i in and C i out are the matrix and sarco/endoplasmic reticulum concentrations of the i-th solute species. The last relationship does not consider the IMM rigidity and the mechanical stress developed in the IMM when it is strained. According to Eq. (18), MS is only complete if On the other hand, such a model takes into account only the osmotic pressure factor. This approach is generally insufficient, as it enters in contradiction with the physical equilibrium concept requiring the osmotic pressure acting on an elastic membrane permeable to water to be compensated by the membrane stress forces. Let us consider an example where a mitochondrial matrix containing known concentrations of different ions is placed into a chamber with pure water, with the volume of water much larger than the mitochondrial volume. According to Eq. (18), we would expect MS continuing infinitely, as the ionic concentrations remain always higher inside the matrix as compared to outside. Thus, Eq. (18) will produce an apparently meaningless result. Therefore, we need a water outflow, which could be induced by negative osmotic pressure (negative terms under the sum in Eq. (18)), and/or by the pressure created by the mechanically strained IMM due to its nonzero rigidity. Assuming equilibrium conditions, the positive sum in Eq. (18) producing an osmotic pressure and a water inflow should be compensated by the water outflow induced by the IMM stress. Note that modeling analysis based on Eq. (18) is described as linear MS, where the rate of mitochondrial volume change is linearly dependent on the osmotic pressure. On the other hand, a linear MS model may be described by Eq. (16) for small IMM strains, where the pressure created by IMM strain is linearly dependent on the strain. It should also be noted that the model represented by Eq. (18) cannot reproduce irreversible MS or simulate the IMM degradation dynamics. As discussed above, the IMM rigidity tensor components defined for an ellipsoidal mitochondrion were represented by Eq. (13), where the nonlinearity of the IMM rigidity was modeled by the following expressions: where the parameter values were determined by fitting the model to the experimental data (see above). The functional form of Eq. (19) was chosen so that the rigidity tensor components would vanish asymptotically at large strains, modeling irreversible swelling. This functional form was chosen based on the ideas of the material resistance theory (Pelleg, 2013). Naturally, an asymptotically vanishing rigidity at large IMM strains implies physical disruption of the IMM, describing irreversible (pathological) MS. The analysis of the presently upgraded model improved the estimates of the parameters describing the IMM rigidity. Further studies are required to introduce an irreversible MS description into the general biophysical model that had been analyzed earlier (Bazil et al., 2010a,b;Pokhilko et al., 2006).
In the current study, we estimated the values of the phenomenological rigidity tensor components. It is important to examine the origins of these phenomenological parameters within the IMM structure. The IMM has a very high protein-to-phospholipid ratio (> 3:1 by weight, which is about 1 protein molecule for 15 phospholipids). Since the detailed structure of the IMM is still under intensive investigation, it is difficult to interpret a certain IMM structure in terms of the rigidity tensor component values. We assume that it is possible to simulate an optimized IMM structure using the Monte Carlo method. In this case, binary interaction potentials between phospholipids and proteins of IMM may be determined. Then, we may use the generalized interaction potential, obtained for the optimized IMM structure, in order to calculate the macroscopic rigidity parameters of the IMM. This analysis should take into consideration the well-known effects of different biomolecules on the IMM rigidity. For example, the cholesterol-to-phospholipid molar ratio affects the rigidity of the membrane. The presently reported upgraded model produced significant deviations in the ΔΨ m dynamics in comparison with the results obtained in the earlier version  at the sarco/endoplasmic reticulum [Ca 2+ ] 0 out = 500 μM. These deviations can be explained by the feedback coupling between the T op and ΔΨ m parameters. Apparently, similar reasoning may be used to explain significant differences in the mitochondrial volume dynamics observed between the two models ( Fig. 3e). Note that the effects of mechanic stress on the IMM behavior were ignored in the previous models, and different system parameters have been reported to oscillate in physiological conditions (Selivanov et al., 1998), although the respective oscillations were quite different from those presently observed. Apparently, such differences may be attributed to the diverse descriptions of the MS dynamics. Another reason could be the different forms of the master equation used in the respective models. Not that most every one of the earlier models used the dynamic equations in the form W X = J X , without expressing the flux J X of the species X in the function of S(t)/V(t). Here W X is the rate of the matrix concentration change of the species X, and S(t) and V(t) are the time-dependent mitochondrial surface area and volume, respectively. Such factorization was used in our analysis, that used time dependence of the permeability coefficients for the species X in the form p X (S(t)/S (t = 0)). Thus, although the reduced permeability coefficient p X /S (t = 0) remains constant, the species water flux in/out of the matrix varies due to changes in the IMM surface area, S(t).
Notably, in addition to the extent of osmosis in the matrix and cytoplasm, physical properties of the IMM including its fluidity, rigidity (stiffness), thickness, and permeability also play an important role in the MS process. There are several factors that may affect the physical properties of the membrane. The most important factor, as was already noted, is the cholesterol-to-phospholipid molar ratio (Shinitzky and Barenholz, 1978). Cholesterol increases the rigidity of the membrane fluid phase, thus the reduction of cholesterol content reduces the membrane rigidity. In contrast, phospholipids increase membrane fluidity and thus reduce its rigidity (Crockett et al., 2001;Wodtke, 1978). Therefore, an increased cholesterol-to-phospholipid ratio should increase the IMM resistance to swelling due to its higher rigidity. In contrast, decreased cholesterol to phospholipid ratio should reduce the membrane rigidity and facilitate its swelling. In addition, fatty acid composition and protein-to-lipid ratio may also affect the IMM physical properties and thus change its response to osmotic gradients. For instance, the specific molecular volume created by phospholipids may vary depending on the saturation and length of acyl groups, which can affect the close packing of phospholipids. Phospholipids containing longer acyl groups will enhance chain-to-chain interactions between fatty acids and thus reduce the membrane fluidity (Ladbrooke and Chapman, 1969). Thus, the IMM physical properties are dependent on a variety of factors, which warrants detailed experimental and theoretical analysis in order to relate the IMM composition and structure with its rigidity parameters.
It should be noted that the results of fitting analysis shown in Fig. 2b demonstrate high accuracy of the experimental data fitting. The average deviation of the fitting curve from the set of experimental data points is around 10%. We expect that fitting results can be improved by developing an upgraded biophysical model that would take into account an extended set of parameters. Future comprehensive modeling analysis of MS will be useful for understanding the molecular identity of the PTP complex and the dynamics of pore opening. In the present model, we made certain progress by taking into consideration the mechanical properties of the IMM. In particular, the estimated values of the IMM rigidity were obtained by fitting of the earlier reported experimental data. These estimated values give a key to the modeling of the IMM organization and to the structural identity of the PTP.

Conclusions
In the current study, we upgraded and further developed our previous biophysical model of MS . It should be noted that the upgraded model produced better estimates for the components of the IMM rigidity tensor, as evidenced by better fits of the present model to the previous experimental data. The new estimates of the IMM rigidity tensor components used in numerical experiments produced significant differences in the mitochondrial behavior in comparison to the earlier reported results. We conclude that the upgraded model produces a better description of the irreversible MS dynamics than the models developed previously. New tools developed and tested in the current study should be further applied for the analysis of irreversible MS. This will result in a more comprehensive model of MS reproducing various physiological and pathological conditions, by taking into consideration fluxes of a large spectrum of ionic species across the IMM and physicochemical parameters of mitochondria. Thus, the current study reports a significant step in the development and testing of the irreversible MS tools and their further application.

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. In addition to the limited number of parameters, the model did not take into consideration the differences between the distinct spatial subpopulations of mitochondria as well as sex-and agerelated differences. However, more complex and comprehensive models of MS may be developed in future on the basis of the current model, using more sophisticated modeling approaches.

Declaration of Competing Interests
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.