Understanding the Dynamics of the Transient and Permanent Opening Events of the Mitochondrial Permeability Transition Pore with a Novel Stochastic Model

The mitochondrial permeability transition pore (mPTP) is a non-selective pore in the inner mitochondrial membrane (IMM) which causes depolarization when it opens under conditions of oxidative stress and high concentrations of Ca2+. In this study, a stochastic computational model was developed to better understand the dynamics of mPTP opening and closing associated with elevated reactive oxygen species (ROS) in cardiomyocytes. The data modeled are from “photon stress” experiments in which the fluorescent dye TMRM (tetramethylrhodamine methyl ester) is both the source of ROS (induced by laser light) and sensor of the electrical potential difference across the IMM. Monte Carlo methods were applied to describe opening and closing of the pore along with the Hill Equation to account for the effect of ROS levels on the transition probabilities. The amplitude distribution of transient mPTP opening events, the number of transient mPTP opening events per minute in a cell, the time it takes for recovery after transient depolarizations in the mitochondria, and the change in TMRM fluorescence during the transition from transient to permanent mPTP opening events were analyzed. The model suggests that mPTP transient open times have an exponential distribution that are reflected in TMRM fluorescence. A second multiple pore model in which individual channels have no permanent open state suggests that 5–10 mPTP per mitochondria would be needed for sustained mitochondrial depolarization at elevated ROS with at least 1 mPTP in the transient open state.


Introduction
The mitochondrial permeability transition pore (mPTP) is a non-specific, transmembrane structure in the inner mitochondrial membrane (IMM) that normally resides in the closed state. Several external factors can cause the pore to switch from the closed state to the open state. The opening can be transient or permanent. Previous studies in isolated mitochondria have suggested that conditions of oxidative stress, depletion of adenine nucleotides, and/or high concentrations of matrix Ca 2+ ([Ca 2+ ] m ) underlie the mPTP opening events [1]. Oxidative stress occurs when the rate of production of reactive oxygen species (ROS) is greater than the rate at which ROS is detoxified, resulting in elevated ROS levels in the mitochondria. This leads to mPTP opening events [2].
The mPTP opening events are considered to have major implications for several biological processes and conditions. For example, permanent mPTP opening has been associated with cell death [3,4]. There are two types of mPTP opening events: transient and permanent. In transient opening events, the pore switches from the closed state to the transient open state and then back to the closed state, while permanent opening events switch the pore to the permanent open state. On a cellular level, permanent mPTP opening events depolarize the IMM, inhibit ATP synthesis and possibly induce cell death. Previous research has shown that when mPTP is in the permanent open state, mitochondria swell, eventually causing necrotic cell death [5]. Conversely, closed mPTPs maintain the membrane potential and integrity of the IMM, preserving homeostatic conditions in the mitochondrial matrix [6]. Once the mPTP opens permanently, the mitochondrion is no longer intact and the IMM is irreversibly depolarized. Since the membrane potential of the IMM drives ATP synthesis, the depolarization caused by permanent mPTP opening stops any and all ATP synthesis.
Although several implications of mPTP have been reported, the behavior of mPTP and the transition between transient and permanent opening events still pose unanswered questions. In this study, a computational model was developed to describe the opening and closing events of the mitochondrial permeability transition pore using recent experimental data from isolated heart muscle cells (cardiac ventricular myocytes) in Boyman et al. (2019) [7]. These data were chosen to constrain the model since it is one of the few studies of mPTP in functioning mitochondria in living cells and the experimental manipulations were minimal so as not to perturb normal function. The model includes TMRM (tetramethylrhodamine methyl ester), a potentiometric fluorescent dye, as an indicator of mPTP opening events, according to the experimental protocol of Boyman et al. (2019). Along with better understanding the correlation between single and multiple mPTP opening events, another major purpose of this model was to augment the interpretation of the mPTP dynamics beyond what was experimentally measurable. Particularly, the model would verify how the sampling rate in the experimental study impacted the experimental data, including whether mPTP transient opening events truly have different amplitudes. The application of Monte Carlo methods simulates the data observed in the experimental study. The usage of the Hill Equation characterizes any cooperativity in the impact of ROS on the opening and closing events of mPTP. Furthermore, Euler's method was used to linearly approximate the reaction rates of the transitions from the closed state to the transient open state, the  transient open state to the closed state, and the transient open state to the permanent open state. In addition, the model incorporated the addition of two different types of inhibitors of mPTP opening, cyclosporine-A (CsA) and n-acetyl cysteine (NAC).

mPTP Model
Data from the experimental study of Boyman et al. [7] constrained the computational model. Figure 1 represents the different factors that drive mPTPs to transition from one state to another. Reactive oxygen species (ROS) from various sources and increased concentrations of Ca 2+ 2+ , high levels of laser illumination can increase mPTP activation because fluorescent dyes such as TMRM, used to measure the potential across the inner mitochondrial membrane (IMM), produce ROS by photooxidation [7]. In order to attenuate the effect of the light intensity on mPTP activation, the experimental study used low levels of illumination, which lowered the amount of ROS produced by "photon stress" [7]. The model included two inhibitors used in the experimental study, both of which lower the probability of mPTP transitioning from the closed state to the transient open state: cyclosporin A (CsA), which reduces sensitivity of mPTP to the driving factors and N-acetyl cysteine (NAC), a ROS scavenger that lowers cellular ROS levels.   (2)). In the equation, V max represents the maximum reaction velocity, [S] is the substrate concentration, n is the Hill coefficient, and K is the Michaelis constant.
The mPTP has been shown to increase its open probability in response to increased ROS [8]. Moreover, antioxidant compounds or ROS scavengers reduce mPTP opening [7]. The model reflects this by the Hill equation (in µM) with parameters (n-cooperativity and K-binding constant) adjusted to reproduce the Boyman The values of n = 4 and K = 0.1 µM were fitted to simulate the experimentally measured transient opening. The V max is omitted here, but it is included in the final transition probability in Equation (6). Transitions to the permanent open state were also assumed to be ROS-dependent and were also modeled with a Hill equation, with the values for the cooperativity (n = 4), Michaelis constant (K = 0.5 µM) and maximal opening rate (V max = 0.001 s −1 ) fitted to match the experimentally measured permanent opening events.

ROS
Experimental studies conclude that the mPTP has an increased open probability in response to elevated extramitochondrial [Ca 2+ ] e . However, experimental quantitative data measuring this is limited. In experiments, Clarke et al. used changes in mitochondrial volume (measured by light scattering) to assess the extent of mPTP opening (both transient and sustained) [9]. Within the context of the experiments, we assume that the rate of shrinkage of preswollen mitochondria is proportional to the fraction of mitochondria with at least 1 open pore. In the current study, the Ca 2+ is fixed so this representation will suffice until better data becomes available. This and other data used for these modeling studies were digitized using the Engauge Digitizer (https://markummitchell.github.io/engaugedigitizer/-last accessed 20 December 2021) and normalized. The binding constant and cooperativity were fitted with the following Hill equation (in µM) using the online curvefitting tool at https://mycurvefit.com/ (last accessed on 20 December 2021) ( Figure 2): The V max is omitted in the normalized equation but will be accounted for in Equation (6). Calcium dependence of mPTP opening based on the assumption that the normalized rate of mPTP opening is linearly proportional to the rate of experimentally measured volume change (blue) [9]. Model fit to data (red).
The opening transition probability P (Equation (6)) is the product of the Ca 2+ -dependent opening rate and the ROS-dependent opening rate multiplied by the integration time step dt. The maximal rate 8.25 × 10 −7 fits the opening rate to the experimental data and captures both the rates due to ROS and Ca 2+ . The value of [Ca 2+ ] e was assumed to be at the physiological concentration of 0.1 µM. Multiplying the closing rate of 0.1 s −1 , which matches the experimental average open time of 10 s, by dt yields the closing transition probability Q (Equation (7)). The transition probability to the permanent open state Z was also ROS dependent to match the experimental data on permanent openings (Equation (8)

Mitochondrial Membrane Potential (ψ)
In the model, the mitochondrial membrane potential (ψ) is set to −180 mV at rest. With mPTP opening, the membrane potential depolarizes. This is based on the modeling studies by Nguyen et al., where the opening and closing of an ion channel, with a large ion conductance similar to mPTP, results in step-like changes in the mitochondrial membrane potential [10,11].

TMRM
TMRM is a positively charged fluorescent dye that crosses the polarized mitochondrial membrane in a concentration dependent manner [12]: For simplicity, the assumption of mass action kinetics yields the following dynamic equation to describe the change in concentration of TMRM inside mitochondria ([TMRM] in ): where [TMRM] Total is the total TMRM concentration, k in is the rate constant for TMRM entry (influx) into the mitochondria, and k out is the rate constant for TMRM efflux. The influx rate constant is k in = 0.0002 × abs(ψ) (11) where abs(ψ) is the absolute value of the mitochondrial membrane potential ψ. The remaining parameters are given in Table 1.

ROS Production
Photoexcitation of TMRM has been shown to produce singlet oxygen which can quickly be converted to any of the common ROS species [13][14][15][16][17]. A minimal model was developed to describe ROS dynamics. The ROS concentration ([ROS]) is assumed to be 0.01 µM under resting conditions. This study assumes that once the experiment starts, ROS will be produced by the interaction of the applied imaging laser and TMRM in a continuous manner with rate k ros. The removal of ROS by a buffering system is assumed to be fast. It is modeled using the rapid buffering approximation where [BUF] Total is the total 'concentration' or capacity of the buffering system and K BUF is the equilibrium constant for the interaction of ROS with the ROS buffering system [18]. The dynamic equation for [ROS] is described by

Numerical Methods
The differential equation describing TMRM dynamics was solved using explicit Euler methods. The model was written in FORTRAN and the figures generated were produced using Excel. The Monte Carlo simulation methods are described below. Model code is provided in Appendix A.

Markov Chain Monte Carlo Method
The Markov chain Monte Carlo method is used to reproduce the stochastic transition between states in the mPTP model. The opening and closing of mPTP is assumed to obey the Markov property-that the reachable states and the probabilities of reaching those states are only dependent on the current state and no other previous states. In this study, this method was used as a means of representing the randomly determined transitions between the closed state, the transient open state, and the permanent open state of the mPTP. A uniform pseudorandom number on the interval [0, 1) was generated on the computer and compared to the transition probabilities. For example, when the pore is in the transient open state, there is a transition probability Q that it will transition into the closed state (Equation (7) above) and a transition probability Z (the difference between Q + Z and Q, as shown by the bracket in Figure 3) that it will transition into the permanent open state (Equation (8)  to the closed state. The same approach is used for when the pore is in the closed state, with the probability that it will transition to the transient open state, P, given by Equation (6) (above). The probability that the pore will transition out of the permanent open state is considered zero.

Results
In the experiments of Boyman et al., laser irradiation of myocytes containing TMRM resulted in efflux of the dye from mitochondria, consistent with depolarization caused by ROS-induced opening of the mPTP [7]. As shown in Figure 4, the experimental data (blue line) showed a maximum decrease in TMRM fluorescence of 54% on average for 478 transient opening events over 180 s, which was matched well by the model for 500 simulated openings (red line). Data was fit using an average transient opening event that started at the 19-s mark and lasted for a duration of 14 s. During this event, depolarization of the mitochondrial membrane potential was presumed to drop from −180 mV to 0 mV.  [7]. Model fit to data with 500 transient opening events (red). Mitochondrial membrane potential during an average transient opening event (gray-right axis) was assumed to exactly follow the timing and duration of an average transient opening event. Figure 5A is a histogram of the amplitude distribution of transient opening events. The amplitudes ranged from 0 to 1, with decreasing numbers of events for increasing amplitude. The model distribution (for 800 events) closely aligns with the experimental data for 478 events with amplitudes greater than 0.3 but contains many more events at lower amplitudes [7]. When the 10 s sampling rate used experimentally was changed to less than 0.1 s the amplitude histogram for all transient opening events changed considerably. As shown in Figure 5B     In addition to CsA, the effect of the ROS scavenger NAC on the observed frequency of transient opening events [7] also was replicated in modeling. The addition of NAC was mimicked by changing the rate at which ROS levels increased, since NAC scavenges ROS. As depicted in Figure 8, the total number of transient opening events over 30 min decreased from 44 events without the addition of NAC to 3 events with the presence of NAC, an 80% decrease. ROS levels were essentially flat over time resulting in the low frequency of transient opening events throughout the experiment.  Figure 9A shows the transition from transient to permanent opening events over time as photon stress increases in different cells, which models the experimental data represented in Figure 9B Figure 10A shows the transition from transient to permanent opening events over time for models as photon stress increases in cells from different groups: control, CsA, and NAC, similar to the experimental data shown in Figure 10B [7]. Over time, the average TMRM fluorescence decreased in cells from all groups. However, a clear distinction exists between the control group and the two experimental groups (CsA and NAC). Whereas the average TMRM fluorescence decreased by almost 40% over the course of 60 min in the control group, the decrease was significantly less in the groups with CsA or with NAC. Figure 10C shows that the addition of CsA or NAC significantly decreases the number of simulated transient and permanent opening events with increased photon stress. In the control, as more mPTP pores enter the permanent open state, the number of openings decreases since there are fewer pores available.    Figure 11C. In this figure, no transition can be seen at all, as all 100 pores continue to stay in the closed state over all 60 min, suggesting that NAC has significantly decreased ROS levels, and subsequently inhibited mPTP openings to an even greater extent than CsA.  where the transition probabilities P and Z are defined similarly as above with the factor 1/n added to the opening rate to keep the number of opening events relatively consistent as the number of channels changes P = The modeling studies next explore the transition from transient to permanent events. The first set of simulations with this model explored how the number of channels affected simulation results. Figure 12A Figure 12B). There is variability in the transition similar to experiment and the previous model. The only observable differences are that there seems to be more positive inflection (caused by all the mPTPs in a mitochondrion being closed) than in experiment and the previous model.   Figure 13A. By contrast, in a mitochondrion with only 1 pore as shown in Figure 13B, there can be times when there is no pore in the open state and the mitochondrion membrane potential can recover until around 50 min. Overall, such simulations indicate that once the number of pores is five or higher, the chance of all pores being closed is unlikely, in effect leading to a transition to a permanently depolarized state.

Discussion
This study examines the behavior of the mPTP through the development of a computational model replicating the results of an earlier experimental study by Boyman et al. [7]. The computational model used Markov Chain Monte Carlo simulation to simulate mPTP dynamics and ordinary differential equations to describe the ROS and TMRM dynamics. A Hill equation describes the dependence of the mPTP opening rates on ROS and Ca 2+ with the parameters fit to match the experimental data. The main purpose was to clarify the interpretation of the experimental data and verify the accuracy of the conclusions drawn in the experimental study. The first model based on observed phenomena was able to simulate the experiments.
The model assumes that ROS and Ca 2+ are independent factors. Some recent studies have suggested that ROS is responsible for mPTP formation and Ca 2+ is primarily responsible for opening [19]. Others have suggested that Ca 2+ -dependent opening of mPTP is modulated by ROS [20]. An alternate formulation of the opening transition is possible. For example, if the level of ROS modulates the Ca 2+ -dependent opening of mPTP the opening rate can be written as (17) where V max is the maximal rate, K ROS is the EC 50 for ROS, K Ca is the EC¬ 50 for calcium, and n and m are Hill coefficients for the ROS and Ca 2+ dependence, respectively. However, in the current study, which focuses on ROS dependent dynamics, the Ca 2+ levels are not varied. Once sufficient additional quantitative data are gathered on the interplay of mPTP opening with Ca 2+ and ROS levels, such a model can be developed. The mPTP is generally thought to be a large conductance channel (on the order of 1 nS) with multiple subconductance states [21]. However, it is also possible that the smaller and larger conductances represent separate physical processes, as suggested recently [22]. Previous studies have suggested that fA currents are sufficient to depolarize mitochondria [23]. The current model does not make any assumptions about pore conductance, but does assume that, once a mPTP opens, the mitochondrion depolarizes. Given the small volume of the mitochondrion, even a single pore opening of 100 pS conductance is sufficient to quickly depolarize the mitochondria, as concluded in our previous studies [11]. Furthermore these previous studies suggest that such mitochondrial depolarizations would produce TMRM fluorescence profiles similar to those observed experimentally [11].
With regards to the idea that transient and permanent changes in fluorescence might represent separate processes, an alternative model can be developed involving different channels described as: The formulations of P, Q, and Z would be similar to those in the current model with some adjustments to match the data. Such a model most likely could reproduce the experimental data. However, at this time we prefer the current, simpler model that explains the progressive nature of the data in terms of a transition (transient to permanent openings) in a single pore entity with increasing [ROS].
The simulations suggest that the sampling rate did have an impact on the amplitudes of the transient events in the experimental study, as illustrated in Figure 5. When a large sampling interval of 10 s was used in the computational model ( Figure 5A), simulation results were similar to those of the experimental study. On the other hand, when the sampling rate was decreased to less than 0.1 s ( Figure 5B), results were clearly different.
Whereas the data in Figure 5A has the appearance of a roughly normal distribution, that in Figure 5B is skewed resembling exponential decay as expected for stochastic ion channel gating [24]. The experimental data miss a majority of the events, particularly those with shorter openings and lower magnitude depolarizations, as indicated by reduction in TMRM fluorescence. Thus, a higher sampling rate would be needed to determine whether events with shorter amplitudes and time spans occur and (as predicted by the modeling) account for most of the total transient opening events. However, while a faster sampling rate is possible, it would, it would entail more frequent laser irradiation and change the dynamics of the system through increase ROS production. The net result would be more sustained depolarization events. The model also predicted that, as the number of mitochondria with mPTP in the permanent open state increases, the rate of transient openings should decrease due to fewer mitochondria with closed mPTP and intact membrane potential ( Figures 9C and 10A).
The computational model replicated the experimental study accurately in several other aspects. As shown in Figure 4 [7], albeit at different calcium levels, would help to model the calcium dependence of mPTP behavior.
In order to get a better understanding of mPTP dynamics, a second model in which there were no permanent open states was presented in Figures 12 and 13. As [ROS] rises the opening rate of the pore increases so that the probability of at least one pore being open is high. In this case the mitochondrion will remain depolarized. If there are only a few channels it is possible that all the channels will close resulting in repolarization and an increasing TMRM signal. For five or more channels this is unlikely. Under the assumption that the depolarization does become permanent without an increase in flickering before this transition suggests that there are at least 5 mPTP per mitochondria. There are probably less than 10 pores since a larger number of pores would lead to excessive flickering before transition to the permanent state. This estimate of 5-10 mPTPs per mitochondrion is consistent with estimates in the experimental literature [25]. Furthermore, experimental data tracking single mitochondria would be helpful to resolve this issue. The model assumes that all mPTP openings are independent, but there is probably some coupling due to the fact that as the mitochondrion depolarizes, the mechanism to replenish ROS removal is attenuated. This coupling would make the transition to the permanently depolarized state more switch-like. Future studies could explore this situation with a more detailed study of ROS dynamics. The The model assumes that the action of CsA simply reduces the opening rate for each channel. The model does not explicitly use CsA binding sites. A more detailed model is possible, which includes the K 0.5 and number of binding sites using mass action kinetics or steady-state approximation. Such a model can be constrained using data from the experimental literature [26]. However, in this research study, the current simplified approximation is sufficient.
There are ongoing studies to determine the molecular identity of mPTP, with various indications that ANT and ATP synthase contribute to the observed phenomena, with regulation by cyclophilin D, a protein "foldase" regulated by CsA [27]. The modeling in this study was agnostic to competing theories about molecular composition of mPTP, but rather modeled the overall behavior of mPTP in terms of the stochastic behavior of ion channels [28]. An expectation would be that with transient depolarizations caused by channel openings, ATP production would decrease until repolarization occurs.

Conclusions
Combining computational models with experimental results greatly helps to clarify how different molecular mechanisms affect the data recorded experimentally. Here, the newly developed computational model for mPTP, based on minimal assumptions about channel behavior, accurately simulated experimental data. The model replicated many features of mPTP openings observed in cardiac myocytes triggered by "photon stress." Importantly, it suggests that mPTP openings should follow an exponential open time distribution predicted by stochastic channel theory. Furthermore, the model suggests that for 5-10 mPTPs in each mitochondrion, transient openings can produce mitochondrial depolarization that mimics a single permanent opening. Finally, the model suggests how to vary conditions employed in experiments to garner further insights into mPTP dynamics.

Conflicts of Interest:
The authors declare no conflict of interest.