Contribution of the Collective Excitations to the Coupled Proton and Energy Transport along Mitochondrial Cristae Membrane in Oxidative Phosphorylation System

The results of many experimental and theoretical works indicate that after transport of protons across the mitochondrial inner membrane (MIM) in the oxidative phosphorylation (OXPHOS) system, they are retained on the membrane–water interface in nonequilibrium state with free energy excess due to low proton surface-to-bulk release. This well-established phenomenon suggests that proton trapping on the membrane interface ensures vectorial lateral transport of protons from proton pumps to ATP synthases (proton acceptors). Despite the key role of the proton transport in bioenergetics, the molecular mechanism of proton transfer in the OXPHOS system is not yet completely established. Here, we developed a dynamics model of long-range transport of energized protons along the MIM accompanied by collective excitation of localized waves propagating on the membrane surface. Our model is based on the new data on the macromolecular organization of the OXPHOS system showing the well-ordered structure of respirasomes and ATP synthases on the cristae membrane folds. We developed a two-component dynamics model of the proton transport considering two coupled subsystems: the ordered hydrogen bond (HB) chain of water molecules and lipid headgroups of MIM. We analytically obtained a two-component soliton solution in this model, which describes the motion of the proton kink, corresponding to successive proton hops in the HB chain, and coherent motion of a compression soliton in the chain of lipid headgroups. The local deformation in a soliton range facilitates proton jumps due to water molecules approaching each other in the HB chain. We suggested that the proton-conducting structures formed along the cristae membrane surface promote direct lateral proton transfer in the OXPHOS system. Collective excitations at the water–membrane interface in a form of two-component soliton ensure the coupled non-dissipative transport of charge carriers and elastic energy of MIM deformation to ATP synthases that may be utilized in ATP synthesis providing maximal efficiency in mitochondrial bioenergetics.


Introduction
One of the main bioenergetics issues is still not fully resolved-how the efficient energy transfer in the mitochondrial oxidative phosphorylation (OXPHOS) system between proton pumps and ATP-synthases is achieved. While Mitchell's chemiosmotic theory provides a general understanding of the OXPHOS system functioning, it ignores a large number of experimental facts showing that after transmembrane transfer or after H + dissociation at the membrane-water interface, protons do not immediately equilibrate with a water bulk phase but are retained at the membrane-water interface in nonequilibrium state (review in [1,2]). The retention of the protons on the membrane was shown in the model bilayer membranes using different techniques for proton release [3,4]. On the liposomes [5] and even on the octane-water interface [6], it was shown that surface protons drive ATP synthesis. The origins of proton affinity for the interfaces were studied both experimentally [7] and theoretically [8][9][10]. The tendency of the water bulk phase to electrical neutrality, which is an important organizing principle [11], provides the simplest explanation why uncompensated H + ions are retained at the interface. The retarded transfer of protons from membrane to the bulk phase was shown in photosynthetic systems [12], rhodopsin purple membranes [13,14] and mitochondrial inner membrane [15]. Not long ago, it was experimentally shown that, under operation of mitochondrial proton pumps, there is no significant acidification of the intermembrane space [16], and there are also indications that there is a lateral proton concentration gradient between proton pumps and ATP synthases [17].
Despite all of the above, the full theory of proton transport on the interfaces is not yet concluded, so there is still a need to develop a satisfactory physicochemical explanation of the proton transfer mechanism in the OXPHOS system on the membrane surface taking into account local environment. It remains unclear how exactly proton-membrane interaction occurs, how much energy can be transferred by nonequilibrium protons and in what form and how this energy is used by ATP synthase. The existing incompleteness of the conventional theory is evidenced by the emergence of suggestions for supplementing the chemiosmotic theory [11,18]. In this work, particular attention was given to the question of how fast proton motion along the membrane is ensured while the proton strongly interacts with the membrane, even changing the structure of the membrane itself [19]. At first glance, it should create high energy barriers for proton transfer or lead to significant energy dissipation. Indeed, in an artificial model system, when the distance between proton pump and ATP-synthase becomes more than 80 nm, energy dissipation led to the impossibility of ATP synthesis [5].
One of the molecular mechanisms underlying the effective proton transport over long distance was considered to be Grotthuss-like mechanism of proton transfer (structure diffusion) in the hydrogen bond (HB) chains formed by water molecules [20,21]. This mechanism was applied to the description of proton conductivity in various molecular systems, i.e., transmembrane channels of bioenergetics proteins in which proton hops from one water molecule or titratable group to the next ones [22][23][24]. It is reasonable to apply the same approach to the consideration of the direct motion of a proton along the membrane, taking into account the latest experimental data that the proton motion occurs at the interphase along bound water [7] but not along the titratable groups [25]. However, it is not excluded that membrane groups can modulate the water HB network formation.
It was recognized that the interaction of the ordered HB chains with their environment (amino groups at a surface of protein channels and lipid polar groups at the membrane interface) plays a significant role in maintaining effectiveness of proton transfer [1,7]. The proton transport along the proton-conducting HB structures requires the excitation of collective modes in these structures that cause in part self-assembling of water molecules into the ordered HB chains and its collective reorientation as a result of structural diffusion (hop-and-turn mechanism). Theoretical consideration of the proton transport as a collective phenomenon emerging in molecular systems has been first carried out in the framework of Davydov soliton theory [26]. Davydov soliton [27] or other similar self-localized states [28,29] and collective excitations [30,31] have long been and still are considered as promising candidates for the role of charge and energy carriers within proteins [32], along the lipid membrane surface [33,34], and on interfaces of biopolymers and artificial membranes [35]. The concept of solitary waves was developed to described the mechanism underlying highly effective transport of charge and energy in molecular systems in the form of collective excitations which are described by autolocalized solutions (solitons) of nonlinear wave equations [36]. It was shown that solitons can propagate in quasi-one-dimensional molecular systems over long distance conserving their shape and energy [37]. The further progress in this direction was made within different modifications of two-component soliton models, where interaction of proton motion with the hydroxyl ion vibration and dynamics of orientation defects in the HB chain were taken into consideration [38][39][40]. Note all these models were developed for the ice-like structures of the HB chains, and their application to biological structures requires further development and detailed consideration of the real structure of the proton environment affecting proton transport. The comprehensive model of the soliton transport was developed and applied to the proton transfer along an artificial assembled monolayer surface made up of sulfonic acid headgroups arranged in a regular hexagonal array [41]. In this model, relationship between the mobility of soliton and macroscopic parameters of the sulfonic acid array structure were established. The next step in this direction should be consideration of the dynamic behavior of the proton environment and its role in the effectiveness of proton transport in biological interface structures.
This paper aims at the elucidation of the molecular mechanism underlying collective dynamics in proton transfer along the HB chains formed by ordered water molecules interacting with the environment-mitochondrial inner membrane (MIM). We based on the new structural data on the OXPHOS system obtained by cryo-electron tomography [42], which allowed constructing a model reflecting structural, function, and physicochemical characteristics of the OXPHOS system and MIM. We suggested that the proton-conducting structures along the cristae membrane connect proton pumps (proton sources) with ATP synthases (proton sinks) and promote direct proton transfer in the OXPHOS system. We developed a dynamics model of long-range transport of energized protons along MIM, which is accompanied by collective excitation of localized waves propagating along the membrane and facilitating proton transfer. Coherent motion of proton in the HB chain and local deformation wave in MIM was described in a two-component soliton model. According to the developed model, the two-component soliton ensures coupled dissipationless transport of charge and elastic energy of MIM deformation to ATP synthases that may be utilized in ATP synthesis and control mitochondrial bioenergetics.

The OXPHOS System Organization in Mitochondrial Inner Membrane Folds
Development of the dynamics model of proton transport in the OXPHOS system should reflect the real structure of that system. In our model, we used the new data on the macromolecular organization of the OXPHOS system and complex cristae ultra-structure obtained by cryo-electron tomography (cryo-ET) [42,43] ( Figure 1A). MIM has many folds which form subcompartments called cristae. The most ordered known components of the OXPHOS system are located on the cristae folds and consist of clustered oligomeric structures of parallel rows of respirasomes and rows of ATP synthase dimers. Even in less compressed and ordered structures, ATP synthase dimers form rows at the edges of the cristae [44]. This makes possible to consider proton transport as unidirectional motion from the pump row to the ATP synthase row ( Figure 1B). The short distance of less than 80 nm between the observed rows of proton pumps and ATP synthase dimers provides the direct and fast transfer of proton to ATP synthases along the cristae membrane avoiding proton surface-to-bulk release [5]. The fact that hydrogen ions move laterally in a thin near-membrane layer of water [7] is also taken into account in our model. Schematic representation of proton transfer along the membrane from the respirasome to ATP synthase is shown on Figure 1C. It should be mentioned that apart from charge carrier transport, pumped protons can also transfer the excess free energy due to their nonequilibrium hydration shell formed on the membrane-water interface [2]. The excess free energy of the energized proton may cause local membrane deformation and induce collective excitation formation, the energy of which, along with the proton, is then utilized by ATP synthase.  [42]. Yellow-ATP synthase dimers, blue-complex I, purple-complex III dimers, green-complex IV, and grey-lipid membrane. (B) A dedicated direction of proton transfer between rows of proton pumps and ATP synthases. (C) Schematic reconstruction of the cluster in the OXPHOS system on the membrane fold and a pathway of the lateral transfer of protons from the respirasome to ATP synthase. The area of increased curvature of the membrane is enriched with CL molecules.

Role of Cardiolipin in Cristae Membrane Structure, Dynamics and Function
The folds of the cristae membrane are enriched with a lipid specific for MIMcardiolipin (CL). The mass fraction of mitochondrial CL in MIM is about 18%, and together with another cone-shaped phospholipid, phosphatidylethanolamine (PE,~34% in MIM), it segregates to the negatively curved monolayer leaflet facing the crista lumen, while the opposing, positively curved, matrix-facing monolayer leaflet contains predominantly phosphatidylcholine [45]. Due to its conical shape, CL especially accumulates in areas with a high curvature of the membrane [46], such as the cristae folds. CL molecules selectively interact with respiratory chain complexes [47,48], mitochondrial transporters [49] and ATP synthases [50,51]. Moreover, CL not only binds to components of the OXPHOS system, but also is essential for their functioning. When CL in submitochondrial particles is damaged by reactive oxygen species (ROS), complexes I, III, and IV become dysfunctional, while the addition of CL restores their functions [52][53][54]. CL is also involved in the functioning of ATP synthase [50,55], nucleotide translocator [56], and is also necessary for the assembly and functioning of respirasomes [57][58][59]. Thus, the areas of high curvature of MIM, where the clusters of the OXPHOS system are located, are enriched in CL, which is involved in the operation of nearly all components of the OXPHOS system [60].
CL possesses at least one or two strong acid phosphate groups [61], capable of hydrogen bond formation with water or hydronium molecules. Due to CL importance for the OXPHOS system and its ability to influence water orientation in the near-membrane zone, we considered its presence and special properties when constructing the model of proton transfer in the OXPHOS system. It is known that CL does not form covalent bonds with positively charged lysine residues of ATP synthase rotor, although it selectively interacts with them [62]. Thus, CL acts as a lubricant, preventing the ATP synthase rotor from strong chemical bonds formation with the environment. It is possible that CL can perform a similar role to ensure the motion of positively charged hydronium.

Grotthuss-Like Mechanism of Proton Transfer along the HB Chains
As was mentioned in the Introduction, one of the molecular mechanisms underlying the effective proton conductivity is Grotthuss-like mechanism of proton transfer in the HB chains which is applicable both to transmembrane channels and water-membrane interface systems.
In our model, we considered Grotthuss-like mechanism of proton motion in the HB chains formed by water molecules strongly hydrogen bonding with the phosphate oxygens at the water-lipid interface ( Figure 2). As mentioned above, MIM and, especially, the vicinity of the OXPHOS system are enriched with CL molecules. We proposed that the HB chain of water molecules may be assembled on the surface of CL domains, which were observed in MIM [63] and model membranes [46] and were investigated by computational modeling in curved two-component membranes [64]. The organized water network in the minimal case may be formed by quasilinear continuous chains of CL molecules. While proton transfer occurs along the water HB chains, the positioning of water molecules and distances between them are controlled by water-lipid interactions. Therefore, the water-lipid interaction on the interface makes it possible to consider the connection between proton motion and elastic deformation of the membrane. Note that, in chloroplasts, phosphatidylglycerol (PG) is the only abundant anionic phospholipid in thylakoid membranes, which is needed for photosystem complex functioning [65]. It may also stabilize water HB networks, as CL does in mitochondria.

A Two-Component Model of Lateral Proton Transport along Mitochondria Cristae Membrane
In this section, we propose a molecular mechanism of long-range proton transport along the mitochondrial inner membrane, which is facilitated by the acoustic soliton moving along the membrane. Assuming molecular structure of the MIM which was described above, we developed a two-component dynamics model of the proton transport The model includes two subsystems. The first one is the HB chain of water molecules bound to the second subsystem including the CL headgroup chains within CL-enriched domains. It should be noted that other phospholipids may also participate in the lipid subsystem, such as PE, which is also known to localize on the cristae folds and to be critical for functioning of the OXPHOS system [66]. Considering that PE headgroup is a zwitterion, we presumed that PE is much less suitable lipid to ensure strong electrostatic interaction with hydronium than CL, while it may participate in the process of lipid membrane deformation in mixture with CL.
For the first subsystem, we wrote the Hamiltonian of the proton of mass m in the HB chain in the following form: where p n characterizes the displacement of the nth proton relative to the middle of the nth HB. The point above p n stands for differentiation by time t. ω 2 1 = K/m ,where K is the stiffness constant of the proton-proton interaction in the HB chain. Consideration of the interaction between adjacent protons provides cooperative behavior and ordered states of the HB chain when protons are occupied either left or right well of the double-well potential U HB (p n , q n ). This proton-proton interaction defines two degenerate ground states, which are shown in Figure 2.
The HB potential U HB (p, q) is a two-dimensional potential surface depending on two variables, p and q, where q is the displacement of hydroxyl ions OH − from its equilibrium. It describes the energy potential of the proton, which is formed by the two neighbor hydroxyl ions OH − . A shape of the potential U HB (p, q) and its dependence on q are defined by a distance between two neighbor hydroxyl ions, R. At distance R in a range of 2.5-2.6Å, the potential U HB (p, q) is a double well, and it is transformed into a single well potential when R is reduced below the critical value R c = 2.5Å [67]. In our model, we used the following approximation of the potential U HB (p, q): where potential parameters ε, p 0 , and q c define its shape [26,39]. At q > −q c , the potential is double well and represents two proton degenerative equilibrium positions at p 00 = ±p 0 1 + q/q c separated by the potential barrier of height ε. At q ≤ −q c , the doublewell potential is transformed to a single well potential. This transformation occurs when the distance between two hydroxyl ions is reduced to less than critical displacement q c . Parameters of U HB (p, q) (Equation (2)) were determined based on the approximation of the quantum chemistry calculation of the HB potential by the ϕ 4 -function for the hydrogen bonding hydronium and water molecules at various distances between them [68]. The HB potential with parameters ε = 10.5 kJ/mol, p 0 = 0.26Å, and q c = 0.2Å describes well a potential change when distance between hydronium and water molecules changes from 2.3Å to 3Å. The HB potential surface U HB (p, q) (Equation (2)) at the defined parameters is shown in Figure 3. In the model, we considered Grotthuss-type conduction mechanism of proton in the HB chain using the double well potential (Equation (2)). Proton transfer is a fast exchange of proton between hydronium (donor) and water molecules (acceptor) that correspond to successive proton jumps from one equilibrium position to another generating hydroxyl and hydronium ions (Figure 2). Proton jumps occur through potential barrier, which depends on the displacement q between neighbour OH − ions according to Equation (2). The effective height of the barrier decreases with decreasing displacement q as at q < 0. The lowering the potential barrier facilitates proton transfer in a local region of compression in the HB chain. To consider the motion of OH − ions, we included in the model dynamic description of the CL molecules. The Hamiltonian of interacting CL in the harmonic approximation may be written as follows: where M is CL mass, and Ω 0 and Ω 1 are the characteristic frequencies of vibrations in the lipid system. The last term in Equation (4) takes into account dispersion of elastic waves in the CL chain. The Hamiltonian of the whole "HB chain + CL chain" system may be written in the continuum approximation as: where c 0 = aω 1 and V 0 = lΩ 1 .
The equations of motion corresponding to the Hamiltonian (5) take the form which lead to the following equations of motion for the coupled proton and CL systems: where Equations of motion (7) and (8) were derived by substitution of the following derivatives of the potential U HB (p, q) to Equation (6).
As a result, two coupled nonlinear differential Equations (7) and (8) were obtained where a link between two subsystems can be derived from the Hamiltonian H int of the interaction between the HB and CL chains: Turning to the new spatial variable ξ = x − Vt in Equations (7) and (8), the following set of two ordinary differential equations may be written: Of particular interest is the excitation moving along the chain at constant velocity V = V 0 . Then, variable q can be expressed from Equation (11): It allows expressing q through p in the form: where Equation (10) after substituting Equation (12) takes the form: where Equation (13) can be written in a form of the well-known equation of the ϕ 4 -theory where the following variable and parameters were introduced At a fixed value of velocity V = V 0 , Equation (14) has the exact soliton solution for proton displacement in the variables x and t: Substituting Equation (16) into Equation (12), we obtained a relationship for lipid headgroup displacement in the form of soliton solution: where q 0 = αp 0 and soliton half-width ∆ is defined as an inverse value of parameter σ 2 : The solution p(x, t) (Equation (16)) for Equation (7), shown in Figure 4A, corresponds to the soliton solution of the well-known equation of the ϕ 4 -theory (see, e.g., [36]) and defines a kink and antikink (upper and lower signs in Equation (16), correspondingly). The solution p(x, t) defines the motion of ionic defect (hydronium) and corresponds to domain wall motion, when all protons on the left side of the domain wall are in the right-hand well of the double-well potential (Equation (3)), and all protons on the right side are in the left-hand well (bottom scheme in Figure 2). The proton (hydronium ion) transport along the HB chain corresponds to the motion of this domain wall of width ∆, when proton moves from the left-hand well to the right-hand well of the double-well potential ( Figure 3). The solution q(x, t) (Equation (17)) represents the acoustic soliton ( Figure 4B) and describes the local deformation (compression) of the membrane region moving with velocity V in the CL chain accompanied with proton motion in the HB chain. We suggested that local compression of the lipid headgroups causes compression in the HB chain, hydrogen bonding to the phosphate oxygens of the CL headgroups. This leads to a decrease in the height of the HB potential barrier due to reduction of the distance between OH − groups less than critical value R c . According to Equation (2) and Figure 3, the HB potential at the distance R in the range of R c is locally transformed from double well potential to the low-barrier HB (LBHB) potential and then to the single well potential that facilitates proton transfer between donor and acceptor OH − groups.
In the model, we considered a special case of soliton motion with the fixed velocity V = V 0 which is the specific velocity defined by the dispersion of lipid molecule vibration in the MIM (Equation (4)). Soliton motion with other velocities needs further investigation of the model and computational solution of the nonlinear Equations (10) and (11). For similar two-component models allowing soliton solutions [26,40], dissipationless soliton has a admissible velocity range v < γc 0 , where parameter γ < 1 is defined by internal parameters of the systems. In our case, we estimated maximum soliton velocity as V 0 = 50-100 m/s using relationship V 0 = lΩ 1 , where the estimated distance between lipid molecules on the MIM surface is l = 0.5 nm for surface area of CL molecules of 0.85 nm 2 [69] and estimation of the characteristic frequency of the headgroup dipole oscillations being on the order of Ω 1 = 10 GHz-100 GHz. The estimated value of V 0 lays in the range of sound velocities 50-300 m/s measured in the lipid monolayers [70] and velocity of the solitary waves observed in liposomes [71]. The soliton half-width value ∆ ∼ = 8.0 nm was estimated in a dynamics model of soliton formation and propagation in domain structures in lipid membranes [34]. At the estimated value of width ∆, the soliton region covers approximately 10 lipid molecules. These values of soliton velocity V 0 = 50 m/s and width were used in calculation of soliton propagation shown in Figure 4.

Discussion
In our model, we considered the complex cristae ultra-structure and used new data on the macromolecular organization of the OXPHOS system obtained in our laboratory by cryo-electron tomography [42]. The ordered clustered oligomeric structure of parallel rows of respirasomes and ATP synthase dimers was observed to be located on the cristae folds of the MIM. The short distance between the observed rows of proton pumps and ATP synthase dimers provides the conditions for direct and fast transfer of protons to ATP synthases along the cristae membrane. Apart from the charge transport, pumped protons can also possess Gibbs energy excess due to their interaction with lipid and water molecules on the membrane-water interface. It is expected that this energy should not be wasted during efficient function of mitochondria, but could be consumed by ATP synthase. The interest to the kinetic coupling in mitochondria [16], which is most likely provided by fast lateral proton transport, is increasing, but the molecular mechanism of this process remains to be fully understood.
To describe the lateral transport of energized protons from proton pumps (sources) to ATP synthases (sinks), we developed a dynamics model of the long-range proton transport along the cristae membrane, which is accompanied by excitation of elastic waves moving along the membrane and facilitating proton transfer. We considered the Grotthusslike mechanism of proton motion in the HB chains formed by water molecules strongly hydrogen bonding with the phosphate oxygens at the water-lipid interface. The HB chain was suggested to be assembled on the surface of CL domains, which were observed in CL-enriched cristae membranes [63]. The developed two-component dynamics model of the proton transport combined two subsystems including the HB chain of water molecules and CL polar headgroup chain. We analytically obtained a two-component soliton solution in the model, which describes the motion of the proton kink, corresponding to successive proton hops from one quasi-equilibrium position to another in the HB chain. The proton kink motion in the HB chain is accompanied by the motion of compression soliton in the chain of CL polar headgroups. The local deformation of the membrane during the passage of the soliton facilitates proton hops due to water molecules approaching each other in the HB chain.
In our model, propagation of the proton-carrying soliton is a dissipationless process, in which a shape and energy of soliton are conserved when traveling from proton pumps to ATP synthases. We considered an extensive exchange of energy between two systems, i.e., the HB chain of water molecules and mitochondria membrane, but did not take into account interaction of the two-component soliton with the environment and processes of soliton thermalization through its interaction with heat bath. The problem of soliton stability and its thermalization were detailed investigated in various approximations and showed soliton stability at physiological temperature [36,40]. We suggested that adiabatic approximation in the model is obeyed, i.e., the two-component soliton travels sufficiently fast from the proton pumps to ATP synthases when compared with the rate of thermal relaxation.
Experimental observation of soliton-type excitations in lipid mono-and bilayers have been carried out in a number of experiments using different methods for excitation and registration of elastic pulses. In the experiments with optical generation of solitary waves, excitation of acoustic soliton-like pulses and their dissipationless propagation have been observed in lipid monolayers at surface pressure above a certain threshold value [72]. Soliton-type excitation of elastic pulse has also been observed in liposomes in temperature region of lipid phase transition and its estimated velocity was about 115 m/s [71]. When the soliton moves, the distance between the acid headgroups decreases, and compression of the lipid is up to 20%. It was proposed that nonlinearity of the membrane compressibility function (depends on temperature and pressure) in the vicinity of the melting transition provides the condition for soliton propagation [71]. The sharp decrease of proton lateral diffusion coefficient after lipid phase transition into the liquid-ordered state was also reported in the monolayer and was associated with hydrogen-bond network breakup [73]. Interestingly, we also registered earlier that effectiveness of mitochondrial ATP synthesis increased near the temperature of lipid phase transition into liquid-disordered phase [74,75]. Note not only the temperature and pressure can affect the phase state of the lipid membrane-it is also controlled by pH (by protons). In accordance to this, propagation of solitary waves was observed on a lipid interface at pulses generationed by local acidification [76]. The pH-induced perturbation of a lipid interface propagated at the velocities up to 1.4 m/s. The obtained correlation between the mechanical (monolayer compressibility), thermodynamic characteristics of the interface, and the pulse velocities confirmed the acoustic mechanism of these pulses and proton transfer along the lipid monolayer. Fichtl et al. proposed that proton soliton-like transport offers effective enzyme-enzyme communication in cells [76]. Note the defined velocity of solitary wave is much slower than the value of sound velocity ranging from 50 m/s to 300 m/s measured in the lipid monolayers [70]. The results of our modeling revealed the significant role of membrane deformation in the mechanism of lateral proton transfer. This finding agrees with experimental data which showed dependence of the proton transfer efficiency on the surface area of lipid headgroups and surface pressure applied to the membrane [73]. The fast proton transport was reported to occur in a defined region of surface pressure and was abolished beyond its boundaries. These experimental data also confirmed our model assumption on the tight coupling between the networks formed by the lipid headgroups and interface water molecules. In the model, strain in the water HB chain (Equation (7)) is defined by the local deformation in the CL chains (Equation (8)) that cause formation of the two-component soliton.
In our model, we considered a significant role of CL molecules in the assembling of proton-conducting structures in the MIM. The experimental and computational works proposed that CL molecules (being unique anionic lipid), besides their essential contribution to mitochondrial structure, functions, and stress response (see Section 2.2), could also play a central role in proton traps at the lipid-water interface and contribute to rapid lateral proton conduction along the membranes [77,78]. Proton-conducting structures in the MIM proposed in our model could be formed as a result of CL domain organization in the membranes which was obtained experimentally [63,79], as well as in computational simulation [69,80].
The developed model of the proton-carrying soliton described a stage of soliton propagation over membrane-water interface and did not consider the initial stage of the soliton formation. This stage should be investigated thoroughly by modeling of the local perturbation of the mitochondria membrane by proton or hydronium. At the current step of the model development, we proposed that proton-carrying soliton might be initially excited by protons pumped to the outer side of the MIM. The local decrease of pH and incompletely compensated highly localized proton/hydronium charge cause weakening of electrostatic repulsion interactions between lipid headgroups that lead to an increase of local lipid packing, accompanied by a decrease in lipid surface area and local deformation of the membrane [19,81,82]. We proposed that the mechanism of proton-induced local deformation of the membranes may be modeled based on the experimental data on the formation of inward plasma membrane curvature at decreasing local extracellular pH [83]. Consideration of this process would allow description of the full energy of proton/hydronium interaction with the membrane, which is stored not only in electrical field, H + concentration gradient, and hydration shell of the proton but also in elastic deformation of the MIM.
Within the developed model, we made several assumptions and did not take into account a number of important features of the MIM and OXPHOS structures. First, we did not consider electrochemical transmembrane potential driving ATP synthesis. Introduction of the electric transmembrane potential is proposed to lead to additional stress of the membrane [84]. The effect of electrical field on membrane deformation was observed, in part, by the method of inner field compensation on the bilayer membranes. Its application showed a change in the membrane thickness and capacitance when the membrane interacts with various compounds [85][86][87], including the membrane-bound proton fraction [3,4]. We proposed that a similar effect may also present in the MIM under proton pump functioning [2]. Consideration of the additional deformation of the MIM in the modified model would allow investigation of the peculiarities of the proton transport along the strained MIM.
Second, a significant factor defining lateral proton motion along the MIM is the local pH gradient of 0.3 units experimentally observed between proton pumps and ATP synthases [17]. Introduction of the corresponding electric field as lateral proton motive force into the model (in right-hand side of Equation (7)) would lead to unidirectional motion of proton-carrying soliton from the proton pumps to ATP synthases.
Third, the developed model did not take into account restoration of water molecule orientation in the HB chain after proton passage. This did not allow describing continuous proton transfer in the HB chain at the water-membrane interface. As seen in the bottom scheme in Figure 2, proton transfer along the HB chain reorients the chain configuration, which should be restored in order to reset transfer of the next proton. The restoration of the initial orientation of the HB chain significantly hinders proton transport following Grotthuss-like mechanism. As reported, fast proton hops between sites (1 ps) according to Grotthuss-like mechanism is followed by much slower reorientation of the HB chain (100 ps) [21,88]. To overcome this problem, additional mechanisms were suggested to reset the HB chain in the soliton models of proton transport and describe continuous proton transfer in the HB chain [89]. We proposed that the energy of water molecule reorientation might be provided by electrochemical potential across the MIM, as well as pH gradient along the MIM. These factors may accelerate relaxation of the HB chain to preferable orientation of electric dipoles of water molecules in the lateral electric field induced by the pH gradient along the MIM (Figure 2A). Consideration of these factors in the model at a further stage of its development will allow a description of the continuous proton conductivity through the HB chains of water molecules between proton pumps and ATP synthases along the MIM.

Conclusions
Based on the experimental data on cryo-electron tomography and nonlinear dynamics model, we developed the model-based approach to the investigation of the coupled proton and energy transport along mitochondria cristae membranes and proposed that proton conductive networks provide a direct connection between rows of proton pumps and ATP synthases assembled on the crista folds. The developed model describes collective excitation as two-component soliton that can be formed and propagate in the strongly linked systems of the hydrogen bond chains of ordered water molecules and membrane lipid polar headgroups. We suggested that the advantage of the proposed soliton-type mechanism lies in achieving coherent transport of proton and energy, which can be simultaneously utilized by ATP synthase. As we established in the model, the cristae membrane, its structure, topology, and lipid composition play a significant role in proton and energy transport and communication between key components of the OXPHOS system. The further investigation of the molecular mechanism of the direct proton and energy transport in mitochondria crista needs detailed consideration of the proton environment including a unique system of the water-membrane interface in the OXPHOS system. Experimental verification of the proposed soliton-type communication mechanism in mitochondria could be carried out by investigation of the effect of proton environment modulation on the effectiveness of ATP synthesis. Another way to verify the proposed mechanism of the coupled proton and energy transport along the mitochondria membrane is the use of bioinspired biotechnology in the investigation of ATP synthesis effectiveness in giant liposomes or in other artificial systems encapsulating ATP synthases and mimicking mitochondrial environment.

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