Cheng et al. (2011) studied the effect of exogenous NO on VDAC purified from rat heart mitochondria. They reported a biphasic effect of PPN on VDAC reconstituted in planar lipid bilayers. Peak inhibition was reported at 25 μM PPN (PPN25) which corresponds to a NO concentration of approximately 350 nM. This is within the range of NO that has been shown to be anti-apoptotic. At higher levels of NO, this inhibitory effect was attenuated and eventually enhanced VDAC conductance at 100 μM PPN (PPN100; corresponding to approximately 900 nM of NO). All these experiments were performed at room temperature and command voltages of −10 mV. For detailed methodology for these protocols see Cheng et al. (2011). Here data reported by Cheng et al. (2011) are analyzed to characterize the effect of NO (possibly via nitrosation) on VDAC purified from primary tissue.

VDAC full length cDNA clone from rat heart with the corresponding GenBank accession number BC072484 was purchased from the Mammalian Gene Collection of Open Biosystems (Huntsville, AL). The coding sequence of VDAC was amplified and cloned by standard PCR methods in the BamHI/Not1 sites of Escherichia coli expression vector pET-21a (Novagen, Darmstadt, Germany). Colonies were screened by restriction digestion and the VDAC sequence together with its His-tag at the C-terminus was verified by sequencing. Construction of the VDAC phosphomimetic, where Ser 137 was mutated to Glu, was achieved by overlapping PCR amplification utilizing the QuickChange Site-Directed Mutagenesis Kit (Stratagene, La Jolla, CA). Under physiological conditions, this Ser residue has been shown to be phosphorylated in rat liver VDAC (Distler et al., 2007), and we have also confirmed its basal phosphorylation in rat heart by tandem mass spectrometry using the Thermofinnigan LTQ XL ion-trap mass spectrometer (MS) in the Biotechnology and Bioengineering Center at the Medical College of Wisconsin and analyzed by multistage-fragmentation analysis (Chesnik et al., 2011).

Induction of VDAC expression in BL21(DE3) cells was achieved by incubating the overnight culture for an additional 4 hours after addition of 1 mM isopropyl β-D-10 thiogalactopyranoside. Cells were then pelleted and lysed, and the VDAC containing inclusion bodies were subsequently isolated by centrifugation. The inclusion bodies were dissolved in solubilization buffer and then loaded onto nickel-nitrilotriacetic acid (Ni-NTA) His•Bind Resin (Novagen) columns. Bound proteins were dislodged from the matrix with a Ni-elution buffer. Subsequently, the purified VDAC was refolded by drop-wise dilution in a 1:10 ratio of elution buffer to the refolding buffer. The un-phosphorylated/recombinant VDAC (rVDAC) and phosphomimetic VDAC (S137E) proteins were reconstituted into planar lipid bilayers to measure their respective channel current activities as previously described (Cheng et al., 2011). The term un-phosphorylated is used throughout the text to indicate the recombinant VDAC.

Briefly, a mixture of phospholipids in chloroform were mixed and dried under a stream of nitrogen and dissolved in n-decane with a final lipid composition of phosphatidylethanolamine, phosphatidylserine and phosphatidylcholine (Avanti Polar Lipids, Alataster, AL) at a ratio of 5:4:1 (v/v). The lipid bilayers were formed in symmetrical solutions, with both the cis and trans chambers containing (in mM): 10 HEPES, 500 NaCl, 1 CaCl₂, 0.2 MgATP, with pH 7.4. The cis chamber was held at virtual ground, and the trans chamber was held at the command voltages. The VDAC protein was incorporated from the cis side. Currents were digitized at 5 kHz and low pass filtered at 1 kHz using a voltage clamp amplifier (Axopatch 200B, Molecular Devices, Sunnyvale, CA) via a digitizer (DigiData 1332, Molecular Devices), and recorded in 1 minute segments at a command voltage of −10 mV. Experiments were performed at room temperature. The pClamp software (version 10, Molecular Devices) was used for data acquisition and analysis is done using a MATLAB based script implementing the developed MCMC algorithm. Although decane-containing membranes have been commonly used in recording the biophysical characteristics of VDAC reconstituted in planar lipid bilayers, it should be noted that the solvent could increase the thickness of the membrane by approximately 1 nm and potentially alter the properties of the channel (Tillman and Cascio, 2003). Despite such limitations, recording VDAC in planar lipid bilayers is a preferred method due to the significant difficulty in measuring VDAC activity in its native environment.

The reconstituted VDAC activity recorded from the planar lipid bilayer experiments is sampled with a time-interval of fixed length (τ) and represented as a sequence of discrete-time channel events I = (I_k)^N_{k = 1}, where N is the number of data-points in the recording and I_k is the single channel current at time t_k. A typical VDAC current recording has multiple open and closed (or minimally open) states with a dominant sub-state that can be regarded as half-open or a half-closed state. Thresholding of VDAC current is employed to identify channel openings and closures which are classified into three distinct conductance states: open (O), sub-state (O_S), or closed (C) event: (1)Ek=E(tk)={O, if |Ik|>|IO|,C, if |Ik|<|IC|,OS, otherwise, where |I_O| > |I_C|. The thresholds I_O and I_C are set at fixed values to identify distinct VDAC conductance states. These threshold values are different for WT and PTMed VDAC activity, however were chosen to be same within a given data type (i.e., WT, S137E, etc.). The threshold values used are tabulated in Table 1. The event sequence (E_k)^N_{k = 1} obtained in Equation 1 can be represented by a Markov model. A Markov model is a symmetric directed graph. The set of vertices contain distinct Markov states {S₁, S₂, … S_n} and the set of edges contain the non-negative constants q_ij (i,j: 1, 2, …, n) that govern the transition rate from state S_i to state S_j. Figure 1 shows different Markov models.

In general, a given channel event can be represented by more than one Markov state. For example, S₅, S₁, S₄ represent the same event O_S for the Markov model shown in Figure 1A. Therefore different Markov sequences (M_k) can yield same event sequence (E_k). To circumvent this problem and identify the best Markov model, Siekmann et al. (2011) construct a probability distribution P((M_k)|(E_k), Q), where Q is known as the infinitesimal generator of a continuous-time Markov model. It is a matrix consisting of entries given by q_ij (i ≠ j) and zero when two states are not connected. The diagonal entries of the matrix Q are given by: (2)qii=−∑i ≠ jqij,i = 1,…,n.

Let p₀ be vector containing the initial probability distributions for the set of Markov states. Then the time evolution of the Markov state probabilities, p(t), can be calculated by the following differential equation: (3)dp(t)dt=Qp(t)⇒p(t)=exp(Qt)p0, where exp denotes matrix exponential which is computed using the Pade approximation (Arioli et al., 1996). Using Equation 3 the transition probability from state S_i to S_j in a given sampling interval (τ) is computed: (4)Aτ=exp(Qτ)=aij,i,j:1,…,n.

To account for different conducting states of VDAC activity the Siekmann et al. (Siekmann et al., 2011) method is modified to introduce three projection vectors θ_O, θ_OS, θ _C which correspond to the three possible events i.e., O, O_S, C. The purpose of these projection vectors is to restrict transition from state S_i(t_k−1) to S_j(t_k) if S_j(t_k) ∉ E_k. This is done by modifying the transition probability matrix A_τ according to an event at a given time t_k, such that: (5)Aτ|t=tk=AτθO if E(t=tk)=O.

Such projection matrices are useful in developing Markov models of channels with multiple sub-conducting states such as, VDAC. Here, the Metropolis–Hastings (MH) algorithm is used for sampling different Markov models (Brooks, 1998; Siekmann et al., 2011). MH algorithm is basically a two-step process:

Independent runs of the algorithm tested different competing Markov models (shown in Figure 1). The best model is chosen as the one that explains all data from all of the experiments on the four different form of VDAC analyzed and based on the uncertainty in the associated parameter distributions. The model that is selected (shown in Figure 1A) invokes a total of eight parameters which are found to be different for WT, nitrosated, un-phosphorylated and phosphorylated VDAC recordings.

To validate the robustness of the Markov model, open-time, half-open-time, and close-time distribution of the VDAC activity are computed and are compared with the open-time, half-open-time, and close-time distributions derived from the experimental data. These probability density functions govern the time for which the channel remains in open state, half-open state or closed state. Colquhoun et al. (Colquhoun and Hawkes, 1981, 1982) derived expressions for these probability distributions making use of the infinitesimal generator Q. In general, the probability density function of open times can be written as: (7)f(t)=Φ0exp(QAAt)(−QAA)uA;Φ0=πBQBA/πBQBAuA.

Here, “A” is the subset of Markov states representing only the VDAC open states, Q_AA = {q_ij : i, j ∈ A} and u_A is a column vector with n_A entries which are all equal to one. Φ is the initial (row) vector that defines the relative probability of an opening starting in a given Markov open state. π_B is the stationary distribution of the Markov states which do not represent the VDAC open state and Q_BA = {q_ij : i ∈ B, j ∈ A}. In the event when there is just one open state (such as shown in Figure 1A) then both Φ₀ and u_A are unity and can be omitted, and the distribution reduces to a simple exponential distribution. For a more comprehensive detail of this approach the reader is referred to Colquhoun and Hawkes (1981, 1982).

In the current recordings analyzed here, the VDAC activity has at least three conducting states which are identified using the thresholds (I_O and I_C). For example, the WT VDAC data shown in Figures 2A,B has three different conducting events which for the Markov model shown in Figure 1A correspond to: open state S₃ if |I(t_k)| ≥ (|I_O| = 20) or closed state S₂ if |I(t_k)| ≤ (|I_C| = 15) or else half-open state(s) S₁, S₄, S₅. Similarly, distinct conducting events for nitrosated VDAC data are also identified using the threshold values listed in Table 1. It is apparent from the purified WT VDAC data (see Figure 2A) that the channel openings and closings are variable in amplitude but model simulations have fixed discrete conductance (see Figure 2A). The simulations assume three distinct conductance states (closed, half-open or open) to match VDAC opening or closing identified using the threshold values.

It is apparent from Figure 2A that the model is successfully able to capture the relative frequency of transitions of WT VDAC to the open and closed state which is apparent from the short 2-second long trace of WT VDAC data and model simulation (see Figure 2B). Nitrosation of VDAC by PPN induces robust changes in its conductance and kinetics which are also successfully captured by the Markov model (see Figures 2C,D). The marked reduction and increase in basal VDAC conductance are captured by changing model VDAC conductance in an ad-hoc manner as described above. The model simulations are able to qualitatively reproduce purified WT and nitrosated VDAC kinetics. The burst-like behavior with PPN 25 μM (see Figure 2C) and WT-like behavior with PPN 100 μM (see Figure 2D) are both successfully emulated by the Markov model. The mean values of the estimated parameters for the two data types are listed in Table 1 and shown in Figure 7A.

The rVDAC activity behaves differently than the purified or WT VDAC activity. Compare the data shown in Figures 2A, 3A. In particular, the rVDAC tends to close more rapidly than the purified WT VDAC. To ensure that the model parameter estimates of the rVDAC data are minimally affected by background noise, the thresholds for closed and open states were slightly increased (compared to purified WT). While determining model parameter estimates of the S137E data, the threshold values were not changed as the background noise level is low (see data shown in Figure 3B). The model simulations of un-phosphorylated (see Figure 3A) and phosphorylated data (see Figure 3B) match the experimental data reproducing the marked increase and decrease in VDAC activity. The rVDAC data and simulation show a high probability of closing while the S137E data and simulation show a low probability of closing. Again, VDAC conductance levels (openings and closings) during model simulations are assigned in a manner similar to the purified WT and nitrosated data. The mean values of the estimated parameters for the two data types are listed in Table 1 and shown in Figure 7B.

The parameter convergence plots are shown in Figure 4A which are obtained from fitting the WT VDAC data, with the Markov model shown in Figure 1A, using the MCMC approach described earlier. The initial estimate of the parameters was obtained by multiple runs of the MCMC method for 10,000 iterations with variable values of δ ∈ [0.0005, 0.05]. This approach helped in obtaining a reasonable initial estimate of model parameters. Posterior probability distributions from the convergence plots are also computed. While computing the posterior distributions of the parameters, first 10,000 iterations are discarded as the burn-in period. Posterior distributions of the parameters provide valuable information regarding the Markov model e.g., model complexity and goodness of the model fits. Figure 4B shows the posterior distributions computed from the convergence plots obtained by fitting purified WT VDAC activity. All of the rate constant histograms have one pronounced peak which represents the best model parameter values. WT VDAC is known to reside predominantly in an open sub-state (S₄ in our model), therefore transitions to a higher conductance state (S₃) or a lower conductance state (S₂) seldom happen. The model captures this aspect of WT VDAC.

The open-time, half-open-time, and close-time distributions are computed from the experimental data (shown in Figure 5). These dwell-time distributions are represented by (square-root transformed) frequency, similar to (Colquhoun et al., 1996), of observing an event of given duration. These dwell-time distributions are compared with the theoretical distribution (η · f (t)) inferred from the Markov model. Here, η is an adjustable parameter accounting for total number of a given event, since very brief events may escape detection due to limitations in the experimental data time-resolution (Lauger, 1988). Dwell-time distributions calculated for purified (Figure 5A) and recombinant VDAC (Figure 5B) confirm the kinetic differences suggested from their single channel recordings. Indeed, rVDAC closes more frequently and has longer lived event durations. Notice the differences in the ordinate limits in Figures 5A,B. The theoretical dwell-time distribution (dashed line) matches well with the dwell-time distributions inferred from the experimental data (see Figures 5A,B).

It is apparent from Figure 2 that PPN 25 μM induced significant changes to the VDAC gating kinetics. In order to quantify these changes we calculated the stationary distributions of the Markov model following each condition WT and PTMs. The stationary distribution suggested that the most dominant Markov state of VDAC is S₄ irrespective of the type of experimental data (WT or PTM). Therefore, it is excluded from the distributions shown in Figure 6 to understand relative dwell-time changes among states: (S₁, S₂, S₃, S₅). Stationary distributions, shown in Figure 6, are the solution of πQ = 0, Q is the infinitesimal generator of Markov chains describing WT and PTM VDAC activities. WT VDAC (Figure 6A) and PPN100 VDAC (Figure 6C) distributions suggest minimal changes in the dwell-times of the four states. On the other hand, there is a contrasting change in the dwell-time distribution for the PPN25 VDAC gating kinetics (see Figure 6B). The dwell-time distribution suggests that VDAC stays ~99% of the time in its closed state (S₂) as compared to other states (S₁, S₃, S₅). Un-phosphorylated and phosphorylated VDAC also share significant differences in their stationary distributions. The S137E data based stationary distributions (see Figure 6E) indicate that there is an appreciable increase in the number of open events, evident from the increase in dwell-time probability of the open state as compared to the control that remains predominantly closed (see Figure 6D). Apparently, nitrosation and phosphorylation of VDAC induce significant changes in its gating kinetics, both of which may have profound effects on mitochondrial (dys)function.

In Figures 7A,B, the mean (with standard error) taken over the parameters estimated from each individual electrophysiological data for a WT, nitrosated, un-phosphorylated or phosphorylated VDAC activity is shown. As purified and recombinant VDAC activities are found to be noticeably different (see Figures 5A,B), purified and recombinant data are shown separately for comparison purposes (Figures 7A,B, respectively). Figure 7A shows purified VDAC associated parameters: WT (Wild-type; n = 4), P25 (PPN 25 μM; n = 5) and P100 (PPN 100 μM; n = 5), and Figure 7B shows recombinant VDAC associated parameters: rVDAC (un-phosphorylated; n = 3) and S137E (phosphorylated VDAC; n = 3). Overall, the estimated parameters are relatively tight. Most significant changes in the estimated parameters of nitrosated VDAC occurred for the PPN 25 μM data. In particular, q₁₃ and q₅₁ have significantly increased, and q₃₁, q₂₄ and q₄₂ have significantly decreased. As a result of these changes the PPN25 VDAC gating kinetics are much faster, as compared to purified WT and PPN 100 μM VDAC data. The phosphorylated VDAC (S137E) closes much less as compared to the recombinant VDAC. This behavior is captured by the significant differences in the q₃₁, q₃₄, q₂₄, and q₄₂ parameters of both the data types (see Figure 7B).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.