Molecular dynamics simulation of proton-transfer coupled rotations in ATP synthase FO motor

The FO motor in FOF1 ATP synthase rotates its rotor driven by the proton motive force. While earlier studies elucidated basic mechanisms therein, recent advances in high-resolution cryo-electron microscopy enabled to investigate proton-transfer coupled FO rotary dynamics at structural details. Here, developing a hybrid Monte Carlo/molecular dynamics simulation method, we studied reversible dynamics of a yeast mitochondrial FO. We obtained the 36°-stepwise rotations of FO per one proton transfer in the ATP synthesis mode and the proton pumping in the ATP hydrolysis mode. In both modes, the most prominent path alternatively sampled states with two and three deprotonated glutamates in c-ring, by which the c-ring rotates one step. The free energy transduction efficiency in the model FO motor reaches ~ 90% in optimal conditions. Moreover, mutations in key glutamate and a highly conserved arginine increased proton leakage and markedly decreased the coupling, in harmony with previous experiments.

. Simulation setup. (A) Side view of the yeast F O F 1 ATP synthase whole complex (PDB ID: 6CP6). The ac 10 part of F O region is depicted as the transparent surface, while the F 1 region and the b-subunit are in the ribbon diagram. The rectangular gray region represents the membrane. The beads inside F O region represent cE59, aE223, aE162, and aR176 (blue). Putatively protonated and deprotonated glutamates are drawn in green and red, respectively. (B) Top view of our simulation system, ac 10 in the F O . We assign integers 1 to 10 to identify each c-subunit in c 10 -ring in clockwise order starting from the nearest subunit to aE223 at the initial configuration. The aE223 site is connected to the IMS side, while the aE162 is connected to the matrix side. (C) Schematic picture of our model. EH and Erepresent protonated and deprotonated glutamates, respectively. Membrane drawn in gray is modeled implicitly. Protons can hop between cE59 and glutamates in a-subunit, aE223 and aE162. Also, aE223 and aE162 exchanges their protons with the IMS and the matrix aqueous environment, respectively. Arrows in green indicate the net proton flow in the ATP synthesis mode. We set the rotational axis of the c 10 -ring as the z-axis and the position of aR176 on the x-axis. (D) Our simulation protocol. We performed short MD and MC simulations in turn. MD simulation updates the protein structure and MC simulation updates the protonation states of 12 glutamates.

Results f o rotation driven by the proton-transfer across the membrane.
We begin with simulations of rotary motions of F O motor driven by a proton-motive force, which corresponds to the ATP synthesis mode. Without external torque applied to the c 10 -ring, we performed proton-transfer coupled MD simulations of the F O motor ac 10 with the pH difference between the IMS and matrix sides being 1 and the membrane potential 150 mV ∆Ψ = . Using the hybrid MC/MD method (see Materials and Methods), we carried out 10 independent simulations with stochastic forces. Figure 2A shows cumulative rotational angles of c 10 -ring in each trajectory together with the number of protons that are transferred from the IMS to the matrix throughout (a representative trajectory in red corresponds to Movie S1). Hereinafter, unless otherwise denoted, we use the rotational angle of c 10 -ring, relative to its initial angle with the counterclockwise direction being positive when it viewed from the matrix side. Clearly, in Fig. 2A, we see that F O rotates stochastically and unidirectionally, and that protons were transferred as F O rotates. Note that here we counted the number of protons that transferred all the way from the IMS to the matrix and thus the protons that were bound in F O in the initial state were not included.
The averaged net number of protons that passed from the IMS to matrix sides, the proton flux, per c 10 -ring revolution was ~10, which clearly suggests the stoichiometric tight coupling between the proton transfer and the rotation (Table S1). www.nature.com/scientificreports www.nature.com/scientificreports/ f o rotates stepwise by 36° coupled with a proton transfer from IMS to matrix. Since the architecture of the c 10 -ring has a 10-fold symmetry, it is thought that F O rotates in °36 steps, and the existence of such steps has been confirmed for E. coli F O F 1 36 . In Fig. 2B, we find the time when the c 10 -ring rotary angle comes over every °36 and the time at which the protonation state changes are synchronized although the rotation steps are not very clear in our simulations. Then, we plotted the population of c 10 -ring rotational angle from 10 trajectories (Fig. 2C), finding that there are 10 low-population gaps in symmetric locations that separate 10 states. This is in a good agreement with the 10 symmetric states observed in the experiment 36 . Notably, the rotational angle distribution of c 10 -ring has relatively broad modes separated by narrow unstable gaps. From this figure, we find the angles at −°6 and °30 correspond to low-population angles (the angle zero corresponds to the initial state). Thus, the c 10 -ring does rotate stepwise by ~ °36 . The °36 step is mainly caused by the low-population bottleneck angles.
Additionally, in the experiment 36 , °72 steps were observed. Our result suggests that the thermal fluctuation can occasionally lead F O to rotate faster and the successive two transitions can be interpreted as one °72 step. Thus, one possible explanation about why they observed 72° step can be the low time resolution compared to the F O rotation.
Pathways in the proton-motive force driven F o rotation. Next, we analyzed detailed pathways of proton-transfer coupled rotary motions in the above simulations. Classifying the pathways, we obtained a few major routes in the ATP synthesis mode (Fig. 2D). Taking advantage of the 10-fold symmetry, here we focus on one step of proton-transfer coupled rotation by °36 . Transitions from the previous step occur around θ −°18 , relative to the angle in the cryo-EM. Given that, the transition pathway starts with the state {DEP(c 1 , c 2 ), θ −°18 }, in which only c 1 and c 2 in c-ring are deprotonated. In the most prominent path (Fig. 2D left), firstly, upon a counterclockwise rotary fluctuation up to θ −°2 from ~18 θ −°, the protonated c 3 gives a proton to the matrix site aE162 and becomes deprotonated, {DEP(c 1 , c 2 , c 3 ), θ −°2 }. Subsequently, it rotates up to θ°3 , and then, the deprotonated c 1 receives a proton from the IMS site aE223 and becomes protonated {DEP(c 2 , c 3 ), θ°3 }. Then, the c 10 -ring further rotates counterclockwise, resulting in {DEP(c 2 , c 3 ), θ°18 }. This final state is equivalent to the initial state with a shift in the chain identify. We call this pathway as Path 1 hereafter.
We also found sub-dominant pathways, Path 2 ( Fig. 2D right) in which the protonation of c 1 occurs before the deprotonation of c 3 .
Additionally, we found the Path 3 which is a small variant of Path 1; the second half of F O rotation occurs before the protonation of the c 1 .
Path 3: DEP(c 1 , c 2 ) 18 the external torque can pump the proton across the membrane. F O F 1 ATP synthase can work as a proton pump driven by the ATP hydrolysis. ATP hydrolysis in F 1 leads to the clockwise rotation of the central stalk, which transmits the torque to F O c 10 -ring. The rotation in F O c 10 -ring results in pumping protons from the matrix side to the IMS side (Note that, for yeast mitochondrial F O F 1 , the proton pump occurs only in vitro, but we use the terms, the matrix side and the IMS side for simplicity). Here, we investigated the proton-pump function of F O driven by an external torque applied to c 10 -ring in the ATP hydrolysis direction. In this study that contains only the F O part, we applied an external torque to the residues on the top of c 10 -ring, cS38, in the direction of ATP hydrolysis, by which we mimic the work by the ATP hydrolysis in F 1 . We turned off the proton-motive force by setting pH in both sides equal and the membrane potential ∆Ψ = 0mV. Figure 3A shows 10 trajectories of rotation angles and the number of protons pumped from the matrix to the IMS environment under . ⋅ 86 2 pN nm torque, which corresponds to the F 1 torque observed in experiments 37 (the proton flux from the matrix to IMS sides is defined as the negative value)(a representative trajectory in red corresponds to the Movie S2). Clearly, protons are pumped up stochastically and unidirectionally across the membrane driven by mechanical forces. At this condition, the average number of pumped protons per revolution was ~10 (Table S1).
Compared to the ATP synthesis process shown in Fig. 2B, there is a clearer stepwise motion in °36 (Fig. 3B,C). It is because the applied torque pushes the c 10 -ring near to the unstable gap found in the ATP synthesis mode, resulting in a sharp and high population there.
Next, we move to a tug-of-war situation: We applied the . ⋅ 86 2 pN nm external torque to F O with the standard proton-motive force; the pH difference being 1 and the membrane potential ∆Ψ = 150 mV. At this condition, the external torque is expected to be stronger than the torque generated by the proton-motive force. Indeed, the simulation resulted in proton pumping coupled with the clockwise rotation of c-ring. Here, the average number of protons pumped per revolution was about 9 with some proton leakage (Fig. 3E). We describe the distinction between the functional proton pumping and proton leakage. In the functional pumping process, a proton first exits from c 10 -ring to the IMS side. Subsequently, a proton enters from the matrix side into c 10 -ring and the proton rotates clockwise with c 10 -ring (the trajectory of each proton is drawn in Fig. 3F bottom left). Rarely, however, after a proton leaves from c 10 -ring to the IMS side, another proton comes in to c 10 -ring from the IMS side. This proton occasionally reaches to the matrix side via a clockwise rotation of c 10 -ring, resulting in the leakage from IMS to matrix side ( Fig. 3F bottom center). Compared with the functional proton pumping, the leakage occurs Pathways in the proton pump motion. Figure 3D shows major pathways observed in the simulations of the ATP hydrolysis mode in the absence of the proton-motive force. For simplicity, we start from the state in which c 2 and c 3 are deprotonated {DEP(c 2 , c 3 ), θ°34 }. Notably, even though the protonation state DEP(c 2 , c 3 ) The positive values mean that the transport from IMS to matrix side and vice versa. "transport", "leak", and "idle"; the definitions are described in (F). (F) Schematic trajectories that protons pass in each process; in the top (bottom) three cases, the c-ring rotates counterclockwise (clockwise). In the productive transport, protons rotate ~°× 36 8 in c 10 -ring; in the leakage process, protons rotate ~ ×  36 2 in c 10 -ring; in the idle process, protons rotate ~ 360 . The yellow-colored processes found in the hydrolysis mode are discussed. (2020) 10:8225 | https://doi.org/10.1038/s41598-020-65004-1 www.nature.com/scientificreports www.nature.com/scientificreports/ has its intrinsic free-energy minimum at θ°36 , the external torque skews the free energy surface so that the minimum shifts to θ°0 . In this initial state (Fig. 3D left), c 1 frequently exchange its proton with the IMS site aE223. If c 3 receives a proton from matrix side aE162 while c 1 is protonated, the protonated c 3 can immediately enter into the membrane, driven by the torque applied to c 10 -ring, and F O rotates in the clockwise direction by °36 . Then the protonation state returns to the state equivalent to the initial state with a shift in the chain identify. We call this Path 4. DEP(c 1 , c 2 )θ −°2 Also, a sub-dominant pathway Path 5 in which c 3 is protonated before c 1 is deprotonated was observed (Fig. 3D center).
The Path 4 is more frequent than the Path 5 because it is easier to receive a proton from the IMS side than from the matrix side due to the difference in the local pH and the membrane potential.
Next, we investigate the pathways in the simulations of the ATP hydrolysis mode in the presence of the proton-motive force, i.e., the tug-of-war situation. Results in Fig. S2 shows that the proton pump activity is much slowed down (Fig. 3E). As described above, we found some proton leak, as well. Since the proton exchange between c 10 -ring and the IMS site is fast, F O occasionally has a proton in c 1 when c 3 is protonated. After the protonation at c 3 , F O can quickly be rotated °36 in the clockwise direction by the applied torque, by which a proton escapes from c 1 to the matrix side immediately ( Fig. 3F the bottom center). In this pathway Path 6, the proton does not go around c 10 -ring and causes proton leakage.
In summary, given the two dominant paths, the Path 1 and the Path 4, there is a tendency that deprotonation of the c 10 -ring occurs before protonation in both ATP synthesis and hydrolysis processes and subsequently F O rotates either in counterclockwise (in the ATP synthesis mode) or clockwise direction (in the ATP hydrolysis mode).
Mutations in f o . So far, we have shown that our computational model shows 10-fold symmetric and stepwise motion consistent with the F O motor, and revealed detailed proton-transfer pathways coupled with rotations that was previously unclear. To validate/support these pathways, we will show that our computational model reproduces the characteristic F O motor behavior in addition to the 10-fold symmetrical motion. Now, we address effects of some mutations on critical residues; cE59 and aR176. It has been experimentally confirmed that a bacterial F O F 1 loses its proton pump function when the corresponding arginine is mutated to alanine 21 . In addition, it has also been confirmed that the ATP synthesis activity of a bacterial F O F 1 drastically decreases when the acidic residue equivalent to cE59A in one of c-subunits is mutated 17 . In order to test whether our model can reproduce these behaviors and in order to investigate the mechanism how these mutations lose the functionality, we carried out MC/MD simulations for these mutants.
The aR176A mutant rotates with some protons leak in the ATP synthesis mode. First, we modeled the mutant in which aR176 is substituted with alanine (aR176A) and performed the hybrid MC/MD simulation with no external torque in the same way as above. Results show that, while the unidirectional rotation does occur, the rotation for aR176A was clearly slower than that of WT (Fig. 4A). Also, the number of transported protons per revolution was about 10 (Fig. 4C).
Interestingly, in the aR176A, the most frequent proton transportation pathway was Path 2, rather than Path 1 that was the dominant pathway in WT (Fig. S3). Namely, from DEP(c 1 , c 2 ), the protonation at c 1 occurs before deprotonation at c 3 . The resulting DEP(c 2 ) state, which corresponds to an intermediate state of Path 2, becomes a broader distribution in rotary angles, which results in increase in the population of this intermediate state and the Path 2. These changes can be understood in the following manner. Compared to the WT case, because the attractive interaction between deprotonated cE59 and aR176 disappears in the aR176A case, the energy barrier on protonation of c 1 , the bottleneck in Path2, disappears. This results in the acceleration of Path 2. We also note that, compared to the WT, the angle distribution in the initial DEP(c 1 , c 2 ) state of the aR176A moved downward, i.e., clockwise rotated ( Fig. S3), which increases the energy barrier for the first step of Path 1. On top, the lack of the attractive interaction slows down the rotation of F O due to the lack of driving force for the counterclockwise rotation.
Broadening of the angle distribution in DEP(c 2 ) of the aR176A mutant increases the proton leakage via a reverse rotation as well. We denote the proton leakage path observed in aR176A without external torque as Path 7 ( Fig. S3B magenta path), In the path, the lack of electrostatic attraction of aR176 also contributes to the clockwise rotation DEP(c 2 ) 18 rot θ −°→ DEP(c 2 ) 34 θ −°. Next, we tested a sensitivity of the mutant with respect to a weak opposing torque, which could mimic the load by F 1 . With a constant and weak torque . ⋅ 8 6 pN nm, we repeated the otherwise same MC/MD simulations. The results (Fig. 4D) show qualitatively similar behavior with a more enhanced difference. While the WT F O rotated nearly at the same angle velocity as the case of no torque, the aR176A mutant exhibited more stochastic and, on (2020) 10:8225 | https://doi.org/10.1038/s41598-020-65004-1 www.nature.com/scientificreports www.nature.com/scientificreports/ average, slower rotation ( Fig. 4D,F). Therefore, the aR176A mutant would not provide sufficient torque to turn γ-subunit in F 1 .
The aR176A does not have the proton pump function. Next, we applied a strong external torque, which corresponds to the ATP hydrolysis mode, to the aR176A mutant F O . While the torque-driven rotation velocity is nearly the same between WT and the aR176A, there is a critical difference in the proton pump function between the WT and the aR176A (Fig. 4G). The proton leakage drastically increased in the aR176A compared to the WT, and the aR176A pumps almost no proton in total (Fig. 4G,I). It has already been assumed that the mutation at the corresponding arginine in a bacterial F O F 1 lacks proton pump activity 38 and our result supports this assumption. From our result, it turned out that actually aR176A transports some protons from the matrix to the IMS sides, but it leaks almost the same number of protons towards the opposite direction, resulting in almost no net proton pumping.
Since the lack of aR176, which means that no barrier between IMS and c-subunit that locates near matrix side, protons can enter into c 2 from IMS side directly before c 3 receives a proton from matrix. Because c 2 locates near at the matrix side, a proton they have can leave from c10-ring and exit to the matrix side. the c 2 E59A mutant rotates with reduce velocity and power. We also performed hybrid MC/MD simulations for the mutant in which one of the cE59 is substituted with alanine (c 2 E59A). Without the external torque, the c 2 E59A mutant can rotate with a reduced velocity (Fig. 4B). The plot shows that the rotation pauses every °360 , apparently corresponding to the mutant stoichiometry of the c 10 -ring. Additionally, the number of protons transported per revolution is about 9 (Fig. 4C), which is a reasonable number because the c 2 E59A lacks one protonatable site out of 10.
Although it is experimentally confirmed that the cE59A mutant loses its ATP synthesis activity, it successfully transports protons from IMS to matrix and rotates in counterclockwise direction in our model. It is because our model lacks F 1 . In reality, F O protein is connected to the central rotor to conduct its torque to F 1 , hence it receives some resisting power to induce conformational changes in F 1 . In order to take this into account, we applied a small torque, . ⋅ 8 6 pN nm and repeated the otherwise same simulations (Fig. 4E). Figure 4E shows that the www.nature.com/scientificreports www.nature.com/scientificreports/ c 2 E59A mutant hardly rotated, although WT rotated nearly equally to the case of no load (gray in Fig. 4B-E). This result suggests that the c 2 E59A mutant is less powerful than the WT, i.e. the mutation makes it impossible to generate a torque strong enough to drive F 1 to synthesize ATP. the c 2 E59A mutant works as a proton pump, though it is less efficient. At last, we applied a strong torque . ⋅ 86 2 pN nm to the c 2 E59A that corresponds to the ATP hydrolysis mode. The cumulative rotational angles are plotted in Fig. 4H. As shown also in Fig. 4H,I, it works as a proton pump though it is less efficient. Since it loses one component of c 10 -ring, the number of protons transported decreases compared to the WT.
The dependency of the rotary mechanism on the computational parameters, pKa. We have shown that our model is able to reproduce the behavior of F O motors characterized experimentally. Next, we address what changes in behavior occurred when pKa, pH, and ∆Ψ values were changed; these values have only been predicted or rather arbitrarily set in the previous reports.
So far, we used one set of pKa values in 12 glutamates based on previous suggestions 20 . However, these values are either merely assumed or inferred from indirect data; there is no direct measurement on these pKa values. Here, to investigate how F O functions can differ depending on pKa values in key glutamates. Here, we changed pKa values of aE162 as in Table S2 and repeated otherwise the same simulations; we performed simulations with no load (as the ATP synthesis mode) and with a strong load, . ⋅ 86 2 pN nm (as the ATP hydrolysis mode). Results in Fig. 5 show that the pKa value of the aE162 does not strongly affect the rotational velocity and the rate of proton transfer in the ATP synthesis mode (Fig. 5A-C). With a smaller pKa value of aE162, the protonated aE162 is less stable, which slows down the proton transfer from cE59 to aE162, but accelerates the subsequent release to matrix aqueous environment by the same factor; the two effects cancel each other. The similar argument could be applied to the pKa value of aE223.
On the other hand, in the ATP hydrolysis mode, the pKa value of aE162 has sharp effects. A smaller pKa value, 8 or 7, decreases the rotational velocity and sharply reduces the rate of proton pump (Fig. 5D,E). In the ATP hydrolysis, a proton must be taken in from the matrix side via aE162. A smaller pKa value of aE162 makes the proton uptake from the matrix side more difficult, thus the rotational velocity decreases. This process is a kinetic bottleneck in the ATP hydrolysis mode.
In the ATP hydrolysis mode, a proton entered into c 10 -ring from IMS side, instead of the matrix side, causes proton leakage (see Path 6). Generally, proton concentration in the matrix side is lower than the IMS side so that proton uptake from the matrix side is slower than that from the IMS side. Therefore, it is hard to strictly stop the proton leakage while keeping the ATP synthesis ability. Actually, as the pKa of aE162 decreases, the probability of aE162 having a proton also decreases and it increases the proton leakage from IMS to matrix (Fig. 5F). With the value 7.0, the number of protons leaked from IMS to matrix is about the same as the number of protons pumped from matrix to IMS. Thus, with that value, F O can virtually pump up proton. Those results show that the range of parameters that allows F O to work as both ATP synthase and proton pump is narrow and our parameter set used here is inside of it. www.nature.com/scientificreports www.nature.com/scientificreports/ The effect of pH. Next, we address how pH in the IMS and matrix environments affect the F O function; notably, the difference in pH values, but not individual pH values, modulates the dynamics. Here, we varied the pH value in the matrix side (see Table S2) and performed otherwise the same simulations as above. Results in Fig. 5 show that a larger value of the pH in the matrix side makes the rotational velocity and the rate of proton transfer larger in the ATP synthesis mode (Fig. 5A-C), as expected from experimental results 39 . The large pH of matrix side facilitates proton release from aE162 and thus the proton transfer from c 10 -ring to aE162 is more efficient. As discussed before, the proton transfer at the matrix side often becomes a rate-limiting process; hence it is reasonable that a large pH in the matrix side increases the rotational velocity in the ATP synthesis mode.
On the contrary, in the ATP hydrolysis mode, a smaller pH value in the matrix side makes the rotational velocity and the rate of proton pumps larger (Fig. 5D-F). The reason coincides with the case of ATP synthesis process. When the matrix side pH becomes larger, the probability of aE162 being protonated becomes smaller, thus the proton transfer from aE162 to c 10 -ring is enhanced. At the same time, proton leakage also occurred via the same mechanism as the case of a large pKa value at aE162 because the probability of being protonated is controlled by the difference between pKa and pH.
The effect of ΔΨ. Finally, we investigate the effect of the membrane potential; Either deleting or doubling the membrane potential (Table S2), we repeated the simulations as before. In the ATP synthesis mode, a larger membrane potential makes the rotational velocity larger in the ATP synthesis process (Fig. 5A-C), as expected from experimental results 39 . It is because the protons are stabilized by being passed to aE223 from IMS and to matrix from aE126.
On the other hand, in the ATP hydrolysis mode, a smaller membrane potential makes the rate of proton pumps larger (Fig. 5D-F). This is because the energy cost to make aE162 protonated decreases. At the same time, the possibility of aE223 being protonated decreases and thus the proton leakage become less probable. It drastically accelerates the proton pumping.
Notably, within our simulation method, both pH and ∆Ψ values appear only in the same formula that represent the probability of protonation at the glutamates in the two half channels (see Materials and Methods: Simulating Proton-Transfer Coupled Protein Dynamics). Thus, the total proton motive force, but not pH and ∆Ψ separately, controls the entire behaviors, which is experimentally well examined 39 .
Protonation state dependent free energy surfaces. In order to investigate why the Path1 is preferred than the Path 2 and the Path 3, we estimated free energy surfaces along the rotation angle for different protonation states (see more detailed method in Section S1). In this analysis, to compare the free energies, we need to keep the total number of protons in the system the same. Thus, for the Path 1/3 and Path2, respectively, we include 9 and 10 protons in the system and explored the free energy surfaces separately (Fig. 6A,B). Figure 6A shows free energy surfaces without the external torque. The first transition in the protonation states in Path1 drastically decreases the free energy compared to that of Path2 (Fig. 6B). Additionally, after the protonation of c 1 , the most stable angle is almost kept, suggesting that the significant rotation occurs after exchanging protons with bulk. The significant rotation might possible before the protonation at c 1 , but it has large energy barrier. This is why Path 1 is more preferred than Path 3. On the other hand, the first protonation state transition in Path 2 slightly increases the energy. This is the reason why Path1 is the dominant pathway compared to Path 2.
Estimation of the torque generated by the F o motor. In order to synthesize ATP using proton motive force, F O needs to generate a torque strong enough to induce conformational change in F 1 . Therefore, how much torque c 10 -ring generates by the proton motive force is important.
To estimate the torque generated by the F O , mimicking experimental setup, we applied an external torque in the direction of ATP hydrolysis and tested if F O can rotate in the direction of ATP synthesis against the external torque, with varying proton motive forces (Note that, alternatively, the torque could be estimated via the free energy surface that takes into accounts the proton motive force). The . ⋅ 86 2 pN nm of torque resulted in the reverse rotation (ATP hydrolysis direction), coupled with the proton pumping, as described. Thus, here we test weaker torques varying between 0 and . ⋅ 86 2 pN nm. F O rotated in the ATP synthesis direction up to . ⋅ 43 0 pN nm of torque transporting protons from IMS to matrix. Whereas, with the . ⋅ 51 6 pN nm of torque, we observed, as the average, the reverse rotation (Fig. S4). An external torque at which the net rotation becomes zero, which corresponds to the torque generated by F O , is estimated as ~⋅ 50 pN nm. Interestingly, our result coincides with the previous theoretical prediction by Oster 40 , from 40 to ⋅ 50 pN nm. From this result, we can estimate the efficiency of the energy conversion from proton gradient to the torque by our F O model. Since c 10 -ring rotates in °36 per one proton, the energy converted to torque per one proton is On the other hand, the free energy loss due to one proton transfer from the IMS to the matrix side is a sum of the energy owing to the difference in pH and the membrane potential. The former is www.nature.com/scientificreports www.nature.com/scientificreports/ The effect of membrane region. As discussed so far, the location of the protein-membrane interface near the matrix channel has an important role for the rotary movement in both ATP-synthesis and proton pumping activity. It is known that mutations in aromatic residues around the matrix channel causes mitochondrial diseases 41 , possibly because the mutation makes a-subunit be distant from c-ring and the lipids intrudes into the gap. Here we investigated this effect by changing the position of the surface of the implicit membrane potential V mem (Fig. S5A). We performed the hybrid MC/MD simulations in the ATP synthesis mode, with varying angles of borders between lipid phase and a-subunit-facing phase; the border is changed from the angle of WT to the angle at aE162 (Fig. S5A).
Interestingly, the rotational velocity of the c-ring decreases monotonically (Fig. S5B) and the number of protons transferred per rotation also decreases (Fig. S5C). This is caused by the wider membrane region. Since the protonated state is more preferable in the membrane, in order to pass a proton to the matrix channel, cE59 needs to escape from the membrane region solely by the diffusion. This is similar to the case of c 2 E59A and the number of transported protons is almost the same as that of c 2 E59A (Fig. 5C).
Especially, when the membrane region locates at the aE162, no protons are transported but protons leak without rotating the c-ring (Fig. S5C). In this case, passing a proton increases the energy when the cE59 approaching from the membrane region but it does not when the cE59 approaches in the opposite direction. As the proton transfer becomes difficult, the leaky pathway becomes more popular. The mechanism is similar to the case of the aR176A mutant, in which protons leak without rotating c-ring.
The result explains why the mutations in the aromatic residues around the matrix side result in loss of its activity. As the a-subunit become distant from c-ring, the membrane-protein interface invades toward the matrix channel and it affects the activity in similar ways to the cases of aR176A and c 2 E59A.

Discussions and conclusion
The current model F O motor reproduced the rotations in the ATP synthesis direction coupled with proton transfers across the membrane driven by the proton motive force as well as the proton pump functions coupled with the reverse rotations when the torque was applied in the ATP hydrolysis direction. Furthermore, the inactivation of ATP synthesis ability and/or proton pump ability for mutants was also reproduced. By visualizing the www.nature.com/scientificreports www.nature.com/scientificreports/ movement of every protons in time domain, we found that protons-transfer from c 10 -ring to a-subunit before proton-transfer from a-subunit to c 10 -ring is the dominant pathway in both synthesis and hydrolysis modes.
In this study, we employed a simple coarse-grained molecular model and combined it to a coarse-grained proton hopping model, which needs some discussions. First, our model contains only Cα atom positions which makes it difficult to estimate the side-chain electrostatic interaction and the steric effect. To circumvent this issue, we used the residue-dependent excluded volume interaction and introduced a direction-dependent kinetic rate which models the relative orientation of the side chains. Second, while we treated only 12 sites as proton-relaying sites in MC protocols, there are other residues which, in principle, can mediate the proton transport. This choice is based on previous studies; c 10 -ring glutamates are mainly discussed in the literatures 20,23,27,40,42 , and aE223 and aE162 are considered as the main proton donor and acceptor to the c 10 -ring 20 . Since we are mainly focusing on the long-time dynamics, we limit ourselves to these representative protonatable sites assuming that the average location of proton in MD steps can be represented by dominant locations. Third, compared to related other approaches, such as the empirical valence bond-based approach 31,32 , our modeling is relatively simple. We considered barriers by a single kinetic rate that depends on the distance and relative direction of side chains and did not include any intermediate protonation states. Clearly, there is much room of improvement. Notably, although both protein structure and proton translocation are simply coarse-grained, still our model reproduced the mechano-chemical features, such as the strength of the generated torque in ATP synthesis mode, external torque driven pumping, and the effect of the mutation. It is worth noting that, not only the strength of the torque, our model reproduced the experimentally-consistent behavior of some mutants supporting that our model captures the important nature of the F O system.
In the ATP synthesis mode, the de-protonation of c 10 -ring cE59 plays a central role in rotational movement, while in the ATP hydrolysis, the rate-limiting event is the protonation of c 10 -ring cE59. It is clear that proton-transfers occur efficiently when the donor and the acceptor are close enough, but, in addition, the range at which c-ring faces the lipid membrane is a crucial factor to determine the pathways. Within the lipid-facing membrane region, cE59 is hardly deprotonated, due to the electrostatics in the hydrophobic lipid regions. Therefore, it may be interesting to investigate the detailed pathways of proton-transfer coupled rotations for F O motors from different species with different numbers of c-subunits.
It has been well established that the c-ring rotation is realized by the basic mechanism that c 10 -ring glutamate (cE59) enters and exits the membrane lipid-facing region through protonation and deprotonation, respectively, and that protons are imported/exported via two half-channels separated by a-subunit arginine (aR176). However, experimentally, it is difficult to directly observe the temporal movement of protons, and thus, it has been unclear in what orders the three events; protonation, deprotonation, and rotation occur. Also, theoretically, mathematical models were not directly verified on the basis of the structure, so it was obscure how the rotation occurs coupled with the number of deprotonated glutamate. In this study, we found that three E59 in c-subunits are deprotonated as the major intermediate state during the rotation of the F O motor. Starting from the c-ring with two E59 deprotonated, the c-subunit near the half-channel connected to the matrix side is deprotonated to reach the state with three E59 deprotonated. Only in the region facing a-subunit, c 10 -ring can be deprotonated because it is not under lipid environment. In other words, since deprotonated c 10 -ring cannot enter the lipid membrane environment, the c 10 -ring cannot rotate in the reverse direction via slippery diffusion. Subsequently, the c subunit near the half-channel connected to the IMS side is protonated. Coupled with that, the c 10 -ring rotates counterclockwise. This dominant proton-transfer pathway in the ATP synthesis mode is an optimal reaction path for efficiently synthesizing ATP because there is no intermediate state that allows slippery and backward diffusive motions for the given F O structure.
Not only functional proton translocation, we also observed that some proton leakage without contributing to the c-ring rotation, especially in the ATP hydrolysis mode. Currently, it is not clear whether this leakage phenomenon may actually occur in vitro or it is an artefact of our simulation setup. The mitochondrial F O F 1 ATP synthase primarily works for and thus may be optimized for ATP synthesis, which might be related to leakage. Since the observed leakage mechanism is composed only of the acceptable proton jumps and forward rotation, it is virtually unavoidable. However, without the electrostatic attraction between the protonated cE59 and aR176, such leakage pathways should have relatively high free energies. Our model might under-estimate the interaction, resulting in an overestimate of the leakage probability.
In the previous work 27 , the authors estimated that the state where only one c-subunit is deprotonated is preferred compared to the state where two c-subunits are deprotonated, based on the free energy estimate. This is primarily by the lack of detailed atomic structure of a-subunit at that time. Very recently, with the atomic structure of a-subunit, the same group found the possible proton pathway includes the states where two c-subunits are deprotonated 42 . Our free energy evaluation (Fig. 6) confirms that the state where only one c-subunit is deprotonated has a higher free energy compared to the state where two or three sites are deprotonated, mainly because of the electrostatic stabilization by the a-subunit.
Limitations of this work. Our current F O model is a minimal model to address proton-transfer coupled F O rotations at structural details. First, for computational ease, we utilized a coarse-grained model throughout, which is lower resolution and unavoidably less accurate than the gold-standard atomistic molecular mechanics. The similar proton-transfer coupled atomistic MD is highly desired in future. Second, we observed characteristic rotational behaviors, however, these motions might be sensitive with pKa values and the treatment of V mem . Third, our model contains only a-subunit and c 10 -ring. However, in the actual F O motor, additional subunits such as b-subunit exist. If these subunits play important roles other than just stabilization of the relative positioning, the ignorance of them could affect the motions/pathways. In addition, because F 1 has 3-fold symmetry, being mismatched to the 10-fold symmetry in c 10 -ring, the whole F O F 1 ATP synthase must show more complex and asymmetric behaviors, which can be a future direction to investigate. Thus, the current study can be viewed as the Scientific RepoRtS | (2020) 10:8225 | https://doi.org/10.1038/s41598-020-65004-1 www.nature.com/scientificreports www.nature.com/scientificreports/ first step of chemical-reaction coupled molecular dynamics simulation, providing a new framework, which has many rooms of improvements

Materials and Methods
Experimental design. Figure 1 shows the overview of the simulation setup used in our research. Fig. 1BC shows the F O ac 10 complex used in our simulation. Figure 1D schematically describes the hybrid MC/MD simulation that alternatively update the protein configuration by MD and the protonation state by MC. More precise methods are described in the following three sections.
Protonation-state dependent energy function. We first describe a general formulation that concisely represents the protonation-state dependent protein energy. We assume that the target protein contains n p amino acids that can take both protonated and deprotonated states, which we call the protonatable sites. We introduce an n p -dimensional vector = { } i , of which each element + h i represents the state of the protonatable site and takes either 1 when protonated, or 0 when deprotonated. Then, we introduce a protonation-state dependent energy function for proteins in general, total n on es es pKa where R collectively represents coordinates of the protein structure. The first term in the right-hand side is the protein interaction energy excluding the electrostatic interactions. The second term represents the electrostatic interaction that depends both on the structure R and on the protonation state + H . The last term is for the intrinsic energy for the protonation of individual protonatable amino acids, which depend on their intrinsic pKa values at given environment.
Protein modeling. Our simulation system contains an a-subunit and 10 c-subunits that form the c 10 -ring.
The structure model is based on the cryo-EM structure of a yeast mitochondrial ATP synthase (PDB ID: 6CP6 20 ) (Fig. 1A). For the c 10 -ring structure, we first trimmed one C-terminal residue in two chains to make the chain length uniform (chains M and R in the PDB file). While the c 10 -ring structure in the cryo-EM model is asymmetric due to the interaction with other subunits, especially with a-subunit, we need a conformation of the c-subunit alone. As a proxy, we obtained the c-subunit structure averaged over all the 10 c-subunits in 6CP6, which was used as the reference structure of the c 10 -ring. For all the proteins, each amino acid is coarse-grained to a single bead located at the position of the Cα atom of the residue. Membrane is treated implicitly via effective potentials and the effect of water solvent is also treated implicitly via Langevin dynamics. . The proteins have both fixed charges and charges on the protonatable sites. The protonatable sites contain 12 glutamic acids according to Srivastava et al. 20 (Fig. 1B,C); cE59 in each of 10 c-subunits suggested as the responsible residue for the proton transport, aE223 and aE162 in the a-subunit which are considered as the proton donor and acceptor, respectively, to the c-subunits along the proton transport pathway 20 . The aE223 resides at the interface to c-ring and connects to the IMS side via the IMS half channel, whereas aE162 sits at the interface to c-ring connecting to the matrix side. We assign −1 charge for the deprotonated state and 0 for the protonated state of these 12 glutamic acids. For other fixed charged residues, we put +1 charge on arginine and lysine residues and −1 charge on aspartic acid and glutamic acid residues. We set the dielectric constant of the system to 10 as an average value in the hydrophobic membrane environment 43,44 . For mutants, we changed the residue type and the charge on it.
The membrane-environment term approximates a free energy cost when a charged amino acid is embedded into the central hydrophobic layer of lipid membrane. Specifically, cE59 sits in the middle of the transmembrane helix, facing outward of the ring form. When the deprotonated cE59 faces directly to lipid, we assign a constant free energy cost, mem  . On the other hand, when the deprotonated cE59 is at the interface to a-subunit, no free energy cost is applied. Thus, the free energy value for deprotonated cE59 differs between at the interface to lipid and at the interface to a-subunit. We model this difference in the free energy as a smooth and implicit membrane-environment potential applied to charged residues.
represents the charge of cE59 and the potential r v ( ) mem i is a smooth function that takes 1 when the residue i faces directly to membrane lipid and switches off at the interface to a-subunit (Fig. S1A). We note that the protonated, and thus uncharged, cE59 has no such energy cost. Specifically, to define r v ( ) mem i , we first define the two border planes (dashed lines in Fig. 1C with the angle 51.2 and −69.6 degree from x-axis) that separate the lipid-facing region (gray shaded in Fig. 1C) and the a-subunit facing region (white in Fig. 1C). For each of the borders, we then define a coordinate r x( ) normal to the border plane and that takes the positive value in the lipid membrane side. We employed a smoothly switching function, (2020) 10 In the full F O F 1 system, the c 10 -ring would be tightly attached to the central rotor of the F 1 subunit to transmit the torque generated by the rotation of c 10 -ring. These interactions together with the peripheral stalk should keep the c 10 -ring close to the a-subunit and, at the same time, also allow it to rotate. Whereas, in our simulation setup, we only have a-and c-subunits. Thus, in order to prevent the c 10 -ring from tilting and diffusing away from a-subunit, we added an extra-restraint that keeps the rotational axis of the ring to the z-axis (Fig. 1C). To define the rotational axis, we choose two residues from each c-subunit that locates the first and third quarter in z coordinate and fit a circle to them. Then we applied harmonic restraints to the residues to keep the distance between them and the center of the corresponding circle. We also applied 3 harmonic restraints to the a-subunit to prevent free diffusion and rotation after aligning aR176 on the x-axis.
The last term in eq. (1), , simply reflects the intrinsic preference of the protonation at each site and can be written as, where pK a i , takes into account the intrinsic pK 0 of glutamic acid as well as its basal environment in Fo motor. Specifically, we set = . cE pK ( 59) 8 0 a , = . aE pK ( 223) 6 0 a , and = . aE pK ( 162) 9 0 a , as the default values, while we will vary these values later to test roles of these values. We chose those values according to the discussion made by Srivastava and coworkers 20 . We note that, as shown in Fig. 1A as red and green beads, all the 12 protonatable sites locate almost the same z-coordinate values in the membrane. Thus, we set the membrane potential of all the protonatable sites identical and thus the membrane potential does not appear in the energy function. Instead, the membrane potential ∆Ψ appears in the boundary condition, as described in the next sub-section.
To identify each subunit in the c 10 -ring, we assign an integer "1" to "10" to each subunit in the clockwise manner looking from F 1 -side (Fig. 1B). Note that these integers are associated with molecules; when the c 10 -ring rotates, these integers rotate together.
Simulating proton-transfer coupled protein dynamics. Given the protonation-state dependent energy function, we now describe the dynamic simulation method. In the current modeling, a proton transfer can be modeled as a change of a pair of charges on two protonatable sites; the proton donor changes its charge from 0 to −1, while the proton acceptor changes its charge from −1 to 0. Because traditional MD simulation is not capable to treat these charge changes, we combined the MD simulation with the MC transitions, which represent proton transfers, to realize the proton-transfer coupled rotation of F O in silico. Short-time MD stages and MC stages are conducted alternatively as in 45 , which is also of some similarity to constant pH MD 29,30 . In the MD stage the protein structure moves with a fixed protonation state, whereas in the MC stage the protonation state changes with a fixed protein structure (Fig. 1D).
In each MD stage, with the given protonation state, we perform a MD run for τ = 10 5 MD steps (Fig. 1D). The MD simulations were conducted by the underdamped Langevin dynamics at a temperature of 323 K. We used default parameters in CafeMol except the friction coefficient, which is set as 2.0 to model the effect of membrane environment. All MD simulations are performed using the CafeMol package version 2.1 46 .
In the MC stage, we first determine the protonation states of the IMS-connected site aE223 and the matrix-connected site aE162. We assume that the aE223 and aE162 sites are able to exchange protons with the IMS and matrix aqueous environments, respectively, and that the proton exchange between these sites (aE223 and aE162) and the aqueous environment is sufficiently fast compared to the timescales of F O rotational motions. Hence, the protonation states on these sites are well equilibrated with the respective environment based on the pKa, the environmental pH, and the membrane potential ∆Ψ. Unless otherwise denoted, we set the pH values of IMS and matrix as 7.0 and 8.0, respectively, and the mitochondria membrane potential as ∆Ψ = 150mV; the membrane potential at the matrix environment is −150 mV relative to the IMS environment. Since the 12 protonatable sites are located near the center of the membrane, we set the membrane potential of all these sites as ∆Ψ/2. Given that, we obtain the local equilibrium probabilities of protonation at aE223 and aE162 as, There protonation probabilities serve as the boundary conditions. Each MC stage begins with sampling the protonation state of these two sites, aE223 and aE162, from these equilibrium probabilities.
Once the protonation states of aE223 and aE162 are set in each MC stage, we then simulate proton transfers among the 12 protonatable sites. Each proton in the protonated site can transfer to a deprotonated site, keeping the total number of protons unchanged. We limit the transfer to be between c-subunit sites and a-subunit sites. For the i-th proton, its transfer probability to the j-th site depends on the kinetic weight factor ← w j i and the energy change upon the proton transfer ∆E, where + H after ( ) and + H before ( ) represent the protonation state after and before the trial of a proton transfer, respectively, and k is the overall rate constant and τ = 10 5 MD steps.
The kinetic weight factor ← w j i reflects the geometry of the donor and acceptor residues. The weight is calculated by a function that depends on the distance r between the donor and acceptor sites as well as an angle θ between two vectors that emulates the conformation of side chain (see Fig.S1B). This modeling is similar to the recent coarse-grained modeling in DNA and RNA base-base interactions 47,48 . The first vector connects the protonation site on c 10 -ring with the site on a-subunit. The second vector connects the middle point of the adjacent residues on the c 10 -ring and the site on c 10 -ring, which indirectly monitors the side chain orientations (Fig.S1B). The functional form is given as, The f is an exponentially decaying function of the distance with the threshold distance = . r 08nm 0 at which the tip of the side chain of glutamic acids are nearly in the direct contact and ∆ = .
r n m 0 4 . The g is a Gaussian function that has the peak at θ 0 , the corresponding angle in the reference cryo-EM structure, and the width σ =°36 . Note that θ is the angle between vectors defined by local conformation as mentioned before and its breadth σ = 36 is unrelated to the rotation angle of c 10 -ring. Notably, the above defined geometric factor ← w j i is symmetric with respect to + H before ( ) and + H after ( ) so that the MC algorithm, together with the Metropolis criterion, satisfies the detailed balance. Given that we do not know the absolute overall rate k, we simply chose it as τ = k 1. The trial proton transfer is accepted by the Metropolis criterion, would not contribute in the energy change.
Additionally, we introduced an extra geometric factor h R176 . It is known that the arginine in a-subunit, aR176, serve to block proton leaks from the IMS to the matrix environment without transferring into c-subunits, which thus separates two half-channels of the a-subunit 21,49,50 . In order to take this effect into account, we put a simple factor h R176 , which is normally unity, but takes zero when the y-coordinates of the aR176 is in between the y-coordinates of the proton donor and acceptor sites (see Fig. 1AB for the definition of the coordinates). For the simulation of the mutant on aR176, we removed this extra factor.
In each MC stage, we repeat these MC trials for every proton in protonatable sites in random order. At the first stage of the simulation, we conducted one round of a MD stage for the equilibration and then we performed 3000 rounds of the hybrid MC/MD simulations.
As the initial condition, we chose a protonation state where the following three sites are deprotonated; the two glutamates of c-subunits, which are facing to a-subunits and the aR162 connected to the matrix side where the proton concentration is lower than in the IMS side. However, the choice of this initial states does not affect the overall behavior we discussed in this paper since the protonation state changes much more rapidly relative to the whole rotation process, and thus is quickly equilibrated.