Brownian dynamics simulation of cytochrome c diffusion and binding with cytochrome c1 in mitochondrial crista

Cytochrome c (Cc) protein shuttles electrons from respiratory chain complex III — from cytochrome c1 (Cc1) subunit — to complex IV during oxidative phosphorylation, in intermembrane space of mitochondria and cristae lumen. With Leigh syndrome (LS), the crista lumen width (CLW) increases, and ATP production declines. One of the questions raised by this situation is to find out how ATP production impairs at LS. Using the simulation of Brownian dynamics, we tested whether the increase in CLW declines respiration at the stage of electron transport of Cc to Cc1. We designed a Brownian dynamics model of horse Cc diffusion and binding with bovine Cc1 in solution by the ProKSim software. The values of the model parameters were estimated to obtain the same dependence of the second-order association rate constant on the ionic strength as in the experiment [1]. Estimated values of the model parameters were used in the model of the reaction in the cristae lumen. The model scene was a parallelepiped. The distance between the two surfaces simulated crystal membranes varied. We received increasing of half-life time of Cc diffusion and binding with Cc1 at increasing CLW. For membrane surface 90Åx100Å (close to the membrane size of complex III), the half-life time of the process changed from 0.098 to 0.22 μs with increasing cristae lumen width from 120 to 160 Å. But due to the half-life time of electron transfer between proteins in the complex, estimated in [1], is higher (100.5μs), the overall time shouldn’t change. To simulate impair of ATP production in the model with an increase in the crista lumen width, we probably need to add to the model IV complex and take into account the dimerization defect of ATP synthase.


Introduction
For most cell types in the human body, mitochondria are the central sources of ATP production. In aerobic respiration, membrane-bound protein complexes pass electrons from NADH or succinate to oxygen, pumping protons across the membrane, and creating a proton gradient used by ATP synthase. The respiratory chain complexes and ATP synthase together form the oxidative phosphorylation (OXPHOS) system -figure 1. Respiratory chain complexes consist of four multi-subunit complexes named complex I -complex IV. The electron transfer from NADH to O 2 releases a lot of energy, which is trapped by the I, III, and IV complexes in transporting protons across the inner mitochondrial membrane for ATP synthase to synthesize ATP. Oxygenic phosphorylation entails the transfer of electrons from NADH to Figure 1. Scheme of oxidative phosphorylation. I -complex I (NADH dehydrogenase), II -complex I (succinate dehydrogenase), Cc -cytochrome c, III -complex III (cytochrome bc1 complex), IVcomplex IV (cytochrome c oxidase), V -complex V (ATP synthase), Q -ubiquinone, QH2 -ubiquinool molecular oxygen. In the respiratory chain, electrons pass from complex I and II to complex III via the hydrophobic electron carrier ubiquinol and from complex III (from subunit cytochrome c1) to complex IV via cytochrome c (Cc). Cc is a 12-kDa soluble protein, diffusing in intermembrane space of mitochondria and cristae lumen. The mitochondrial respiratory chain complexes reside in the inner mitochondrial membrane or cristae. Cristae are the folds of the inner mitochondrial membrane. The respiratory chain complexes form supercomplexes (SC) in the inner mitochondrial membrane [2,3].
Many studies refer to an altered mitochondrial ultrastructure as an indicator of mitochondria dysfunction. Changes in cristae number and shape define the respiratory capacity as well as cell viability [4]. Leigh syndrome (LS) is a mitochondrial disease affecting mitochondrial energy production via OXPHOS. It caused by several mutations and is a progressive impairment of cognitive and motor functions and premature mortality [5]. In [6] using cryoelectron tomography, were found deep ultrastructural defects in mitochondria of patients, including bloated balloon-like cristae. Crista lumen width CLW (the distance between the two opposite sides of the cristae lumen) was significantly increased: 164±59 Å in patient versus 120±32 Å in healthy control.

Description of the model
We designed a Brownian dynamics computer model of diffusion and binding of oxidized horse Cc molecules with water-soluble part of reduced bovine cytochrome c1 (Cc1) in solution by the ProKSim software [7]. We used structures of Cc1 with PDB ID 1BGY [8] and Cc with PDB ID 3O1Y [9]. In the model, when two molecules approached each other, they were directed by the electrostatic field and could occupy a favorable position for docking. In the right orientation, proteins form the so-called transient, or encounter, complex [10]. Further conformational changes in protein molecules result in the formation of the final complex, in which the reaction of electron transfer takes place. In our model, we consider only electrostatic interactions and do not explicitly simulate the final complex formation and electron transfer. We assume that the transient complex transforms into the final one faster than the transient complex forms. Dark grey (blue) corresponds to potential + 6.5 mV, light gray (red) to -6.5 mV. Heme is shown by sticks, the Fe atom -as a sphere.

Electrostatic potentials of proteins
Charged amino acid residues and partial charges of the proteins produce a heterogeneous electrostatic field. If a protein is far away from other proteins, it moves under the effect of the Brownian force. When two protein molecule comes close to each other, their behavior is directed by the electric field of both particles, so that a favorable position for the formation of a transient complex can be achieved. The Poisson-Boltzmann formalism was used to determine the electrostatic potential field generated around the proteins. The step of the grid for calculation potential was 1 Å, the dielectric constant of the cells inside the protein was assigned 2; inside the solution -80; at the boundary protein solution -40; and the value of the ionic strength inside the protein was assigned 0. The electrostatic charge distributions on protein atoms are according to the CHARMM27 field [11] and those on the heme, according to [12]. pH in the model was 7. Protein potentials and the molecules were visualized in PyMol (http://pymol.org/). Figure 2 presents equipotential surfaces of Cc, Cc1 at pH =7 and I=130 mM. Dark grey corresponds to positively charged regions, light grey to negatively charged ones. The Cc1 molecule is negatively charged (net charge -14) and its part, facing the cristae space in mitochondria and binding Cc molecule, has a region of negative potential ( figure 2A). On its other part, the molecule has small regions of positive potential. The Cc molecule is positively charged (net charge +8), and its cofactor heme is located in the region of positive potential ( figure 2B). The molecule has two negatively charged regions. This inhomogeneity of the equipotential surface should facilitate the proper orientation of Cc relative to its protein partners. The reaction rate decreases with increasing ionic strength, indicating that the interaction has an important electrostatic component [1].

The model scene in solution
We estimated values of the model parameters (the docking distance and the docking energy) to get the same dependence of the association second-order rate constant on the ionic strength as in experiment at pH 7 [1]. If at some point in time the distance between Fe-atoms of hemes in cytochromes was less than the docking distance, and the electrostatic energy of The paper [1] shows values for the second-order reaction between bovine III complex and horse photooxidized Ru labeled Cc above 130 mM ionic strength since at lower values of ionic strength kinetic curves are determined by intracomplex electron transport. In computer experiments in cristae we took I=130 mM for approximation the reaction at physiological ionic strength (i.e., I =100-150 mM [13]).

The model scene in cristae
To study kinetical characteristics of Cc and Cc1 complex formation in crista lumen, we used the parameter values estimated for the buffer solution. The model scene was a parallelepiped with soft boundary conditions. CLW -the distance between 2 surfaces, simulated crista membranes, varied. Cc1 molecules were fixed at their initial positions near the center of one crista membrane -figure 4.
Cc molecules diffused under the action of the random Brownian force and electrostatic force and formed a transient complex with Cc1. We repeated simulation 20000 times. We calculated the number of complexes formed at each time step, normalized kinetic curves obtained are shown in figure 5A. From kinetic curves, we estimated the half-life time of the process t 1/2 -figure 5B.

Results
In one part of model experiments in the cristae, the size of the cristae membranes had a fixed value of 90Åx100Å (close to the membrane size of complex III) and CLW changed from 92 to 300 Å. We got approximately linear increasing of t 1/2 on CLW -circles in figure 5B. That is because the second-order association rate constant is independent of protein concentration, and half-life time is inversely proportional to the product of the association rate constant by concentration and so proportional to CLW. The half-life time of the process changed from 0.098 µs to 0.22 µs with increasing cristae lumen width from 120 to 160 Å (as at LS [6]). For CLW 300 Å the t 1/2 was 1.7 µs. The half-life time of electron transfer between proteins estimated in [1] as ln2/(6900 s −1 )=100.5µs is higher than that time. That means that the  electron transfer is the rate-limiting step of Cc-Cc1 interaction. So, the electron transfer rate from Cc1 to Cc, for this size of cristae partition, is determined by the rate of electron transport, and should not depend on CLW. In this case, reaction volume changed linearly on CLW.
To exclude the effect of varying concentration at changing CLW, in another part of model experiments, we changed the size of cristae membranes with changing CLW to get the volume as at case with dimensions 90x100x120 Å 3 . In this case, the dependence of t 1/2 on CLW deviates from linear in contrast to the case with the changing volume -cross in figure 5B. In the third part of model experiments, the surface size of the cristae membranes was constant and bigger than dimensions of the respirasome (190Åx300Å [3]) -900Åx1000Å. In the last case, the half-life time of Cc diffusion and binding with Cc1 increased from 42 µs to 56 µs with increasing cristae lumen width from 120 to 160 Å. So, the reaction became non-limited by electron transport but less sensitive to CLW.

Discussion and conclusions
On the Brownian dynamics model of Cc diffusion and binding with Cc1 in the cristae partition of complex III size, we didn't get the impairment of the overall electron transport in the respiratory chain with the increase in CLW. The impair of ATP production may be caused by: 1) the dearth of cristae [1]; 2) the disruption of ATP synthase dimerization [1]; 3) the defective assembly of respiratory chain complexes in SC that may decrease the respiratory chain efficiency [14]. In future studies for the simulation of the drop of ATP production under the condition of the increased crista lumen width, we will add to the model complex IV and the defect of ATP-synthase dimerization.