Assembly of Spinach Chloroplast ATP Synthase Rotor Ring Protein-Lipid Complex

Rotor ATPases are large multisubunit membrane protein complexes found in all kingdoms of life. The membrane parts of these ATPases include a ring-like assembly, so-called c-ring, consisting of several subunits c, plugged by a patch of phospholipids. In this report, we use a nature-inspired approach to model the assembly of the spinach (Spinacia oleracea) c14 ring protein-lipid complex, where partially assembled oligomers are pulled toward each other using a biasing potential. The resulting assemblies contain 23 to 26 encapsulated plug lipids, general position of which corresponds well to experimental maps. However, best fit to experimental data is achieved with 15 to 17 lipids inside the c-ring. In all of the simulations, the lipids from one leaflet (loop side of the c subunit) are ordered and static, whereas the lipids from the other leaflet are disordered and dynamic. Spontaneous permeation of water molecules toward Glu61 at the active site is also observed. The presented assembly approach is expected to be generalizable to other protein complexes with encapsulated lipid patches.


INTRODUCTION
Rotor ATPases are large multisubunit membrane protein complexes found in all kingdoms of life that convert the energy stored in transmembrane potential into the energy of a covalent bond in adenosine triphosphate (ATP), or vice versa (Walker, 2013;Junge and Nelson, 2015). Generally, ATPases consist of two large parts, one of which is soluble and another one is embedded in the membrane. The membrane part includes a so-called c-ring-a symmetric homo-or heterooligomer of subunits c, which is rotating in the membrane during the protein operation ( Figure 1A). Currently, c-rings consisting of 8 to 17 protomers are known (Kühlbrandt and Davies, 2016;Schulz et al., 2017). The number of c subunits defines the energetics of ATP synthase because it corresponds to the number of ions (H + or Na + ) transported per a full 360 • rotation, which invariably results in synthesis or hydrolysis of 3 ATP molecules. Most of the subunits c consist of two transmembrane helices with the connecting loop directed toward the soluble part of ATP synthase, whereas others have four transmembrane helices, which presumably results from gene duplication. Some of the c-rings can self-assemble (Arechaga et al., 2002) whereas others require accessory subunits such as UncI (Ozaki et al., 2008). The inside pore of the c-rings in some cases is plugged by phospholipids Oberfeld et al., 2006).
Highly detailed atomistic structures of isolated c-rings and whole ATP synthases have recently became available (Kühlbrandt and Davies, 2016;Guo and Rubinstein, 2018;Murphy et al., 2019). High resolution crystallographic structures reveal both the overall fold as well as the geometry of the active site comprising a carboxylate amino acid, which is de-and re-protonated during FIGURE 1 | Assembly of the spinach ATP synthase rotor ring protein-lipid complex. (A) Structure of the spinach ATP synthase (Hahn et al., 2018). The rotor ring is highlighted in green and the approximate membrane boundaries are shown with black lines. (B) The rotor ring assembly approach. c 7 oligomers are separated by 5 nm and pulled toward the experimentally determined positions using harmonic restraints. (C) Assembled c 14 ring. (D) Convergence of the c 7 oligomers during pulling simulations. Yellow and magenta traces converge to 3 Å; the green simulation is trapped in a high-RMSD conformation (4.6 Å); the blue simulation doesn't converge in 1,500 ns. (E) RMSD of C α positions in unrestrained atomistic simulation of the assembled c-ring with 17 lipids inside relative to the starting experimental structure (PDB ID 6FKF). the protein operation. Structures obtained using single particle cryo-electron microscopy (EM) currently have lesser resolution, but are better at revealing the overall complex structure, often in several conformations. Interestingly, in some recent EM structures of yeast V-ATPase, additional transmembrane (TM) helices have been observed inside the c-ring (Mazhab-Jafari et al., 2016;Roh et al., 2018). While the protein components of ATP synthases are generally well resolved in the structures, the densities for the surrounding lipid molecules, including those presumably trapped inside the c-ring, are of lesser quality. In the best available structure, that of Enterococcus hirae V-ATPase K 10 ring, electron densities of 20 diacylglycerol/phosphatidylglycerol molecules are resolved in the immediate vicinity of the protein inside the ring, but not in its center or outside the ring (Murata et al., 2005). Most probably this is due to the great flexibility and dynamic nature of lipids (Buslaev et al., 2016) and lack of specific lipid-protein interactions, or due to detergent extraction used for sample preparation. Consequently, the identities of lipids, and their exact positions, are not well resolved in the available experimental density maps, and corresponding moieties are not modeled in experimental structures.
The lipids in and around c-rings were simulated in several previous studies, for example (Aksimentiev et al., 2004;Sengupta et al., 2009;Krah et al., 2010;Pogoryelov et al., 2010). In the smaller c 10 rings, usually, four lipids were modeled inside (Aksimentiev et al., 2004;Sengupta et al., 2009;Zhou et al., 2017a). Larger c 14 and c 15 rings were simulated with 12 and 15 lipids in the central plug, respectively . Finally, an even larger K 10 ring (40 TM helices) should contain even more lipids (Murata et al., 2005;Leone et al., 2015). Molecular dynamics (MD) simulations have also been used to probe the specific interactions of cardiolipin with the outer surface of metazoan c-rings (Duncan et al., 2016).
Recently, a Cryo-EM structure of spinach chloroplast ATP synthase with recognizable densities in the membrane region has become available (Hahn et al., 2018), providing an exciting opportunity to compare the modeled lipid positions with the experimental data. Neither this structure, nor previously determined crystallographic structures of similar chloroplast c 14 rings (Vollmar et al., 2009;Saroussi et al., 2012;Balakrishna et al., 2014), contain lipid molecules. Here, we use an approach inspired by the natural assembly process to model the lipids inside the c-ring. It builds upon the existing approaches for insertion of preformed oligomers (Kandt et al., 2007;Lindahl and Sansom, 2008;Wolf et al., 2010;Biggin and Bond, 2015;Javanainen and Martinez-Seara, 2016), and uses a biasing force to assemble the whole ring, essentially by incorporating experimental restraints into the simulation. Using the approach, we observed trapping of up to 26 palmitoyl-oleoyl-phosphatidylcholine (POPC) lipids inside the c-ring, and positions of these lipids correspond well to the experimentally determined density maps. Then, we compare simulations with different numbers of lipids inside the c-ring to experimental maps and find the best fitting ones.

RESULTS AND DISCUSSION
Our approach to assembly of the c-ring protein-lipid complex is inspired by the natural assembly process. In nature, probably, subunits c are synthesized by ribosome independently and inserted into the membrane one-by-one. After that, the protomers assemble into oligomers, sometimes observed in atomic force microscopy images Pogoryelov et al., 2005), and eventually form complete rings. Some of the c-rings can assemble in liposomes in vitro while others require presence of accessory proteins (Rühle and Leister, 2015). On the contrary, the existing computational methods for insertion of membrane proteins into membranes involve quite unnatural assembly steps: most of the approaches artificially insert fully assembled proteins and protein complexes into pre-equilibrated membrane patches, whereas others rely on the self-assembly of the membrane around the membrane proteins or their complexes  7, 4.6, 4.5, 4.6, 4.5 10, 9, 9, 10, 9 15, 14, 13, 15, 13 25, 23, 22, 24, 22 DPPC 0.01 5 (10) 33 1.0, 1.4, 2.2, 0.8, 0.6 4.5, 3.0, 4.5, 3.0, 3.4 11, 10, 10, 10, 10 14, 14, 15, 14, 14 25, 24, 25, 24, 24 Convergence times, RMSD and numbers of lipids are indicated only for the trajectories where complete assembly has been observed. ( Kandt et al., 2007;Lindahl and Sansom, 2008;Wolf et al., 2010;Biggin and Bond, 2015;Javanainen and Martinez-Seara, 2016). In most cases, these unnatural approaches work very well, but in special cases, where some of the lipids are encapsulated within the membrane protein or protein complex, special care may be required.
Since the c subunits are relatively simple hairpin-shaped proteins, and there are no lipid molecules trapped in the c:c interfaces, we skipped the initial oligomerization steps in our simulations and started directly from partially assembled rings. Partially assembled rings have been simulated before, when the true stoichiometry of the ring was not known (Schlegel et al., 2012) or in a permeability study (Zhou et al., 2017b). Spinach c-ring contains 14 protomers, which we separated into two heptamers for simplicity of the assembly protocol. While it is likely that in nature the assembly proceeds via oligomers of different stoichiometry, we believe that the resulting arrangements of encapsulated lipids should be similar for different assembly pathways. The c 7 heptamers, separated by 5 nm, were inserted into the membranes consisting of the commonly used model lipid POPC and were pulled toward the experimental positions by a biasing potential realized as a harmonic restraint to the experimental structure (Figures 1B,C). While theoretically the halves could assemble into a complete rotor without any biasing force, in practice this might require a large amount of computational time. Using a force too weak might also result in slow assembly. On the contrary, applying a force too strong could result in unnatural deformation of the bilayer and trapping of a system in unnatural conformation. Therefore, we tested several harmonic biasing potentials with spring constants of 0.002, 0.005, 0.01, and 0.02 kJ/(mol·nm 2 ) and followed root-mean-square deviation (RMSD) of atomic positions relative to the reference experimental structure (Hahn et al., 2018). Using a coarse grained (CG) force field allowed for much faster assembly simulations, due to both computational efficiency and faster dynamics in CG simulations (Marrink and Tieleman, 2013; Buslaev and Gushchin, 2017). In many of the simulations, the structures converged into an assembled ring in <1.5 µs. The numbers of lipids trapped at the loop side and the N-and C-termini (NC)-side were consistent in different simulation runs: 9 to 11 and 13 to 15, respectively ( Table 1). Using stronger restraints often resulted in trapping of the complex in a high-RMSD conformation (∼4.6 Å), whereas using weaker restraints resulted in most cases in absence of convergence ( Figure 1D and Table 1).
After the heptamers have approached each other, we followed the initial weak restraint simulation step by a series of three consecutive 200 ps simulations with strong restraint spring constants of 10, 100, and 1,000 kJ/(mol·nm 2 ), after which the RMSD relative to the experimental structure dropped below 0.8 Å. After that, the CG positions of plug lipids were converted into all atom (AA) representation using backward (Wassenaar et al., 2014), and the CG protein model was replaced by the experimental atomic structure. The resulting complex was solvated by water molecules and simulated without restrains. First, we conducted an AA simulation with an average amount of plug lipids (25). The c-ring remained stable but some lipids from each side of the ring were pushed out of the c-ring core (Figure 2A). Therefore, we have attempted a second simulation, starting from the converged structure with 23 plug lipids. During the strong restraint preparatory CG simulation, three lipids were also pushed out of the c-ring ( Figure 2B). These short simulations suggested that our assembly simulations overestimated the number of lipids inside the spinach chloroplast c-ring, and prompted us to conduct a quantitative comparison of the lipid positions obtained in simulations to experiment. The original structure (Hahn et al., 2018) is not symmetrical due to presence of ATP synthase subunits other than the c-ring. Yet, the ring itself in the experimental structure is symmetric, and a free c-ring in a lipid membrane is expected to be symmetric as well. Thus, for analysis, we averaged the experimental Cryo-EM map according to the 14-fold symmetry of the ring (Figure 3). We have conducted a number of 1 µs long CG simulations of c-ring assemblies where lipids were removed from each side one-by-one, while the protein particles were restrained to the experimental positions. The densities obtained in simulations were averaged over the trajectory length and also over the 14-fold symmetry of the ring. The best fit is observed visually for the system with 6 or 7 lipids at the loop side and 9 or 10 lipids at the NC side ( Figure 3D). Real space correlation coefficients (RSCC), used routinely to compare EM maps to each other and to evaluate model-density fits (Afonine et al., 2018), highlight the 6:9 system as the best fitting one with the RSCC value of ∼0.75 ( Table 2).  (Hahn et al., 2018) contoured at the level of 3σ (black). Excellent correspondence between positions of protein atoms (green) and experimental map is observed. The lipids (magenta) fit very well into the densities in the membrane region. (B) Density distributions of outer (left) and inner (right) POPC phosphate (red) and diacylglycerol (blue) moieties overlaid with three exemplary trajectories of permeating water molecules. Average position of Glu61 centers of mass is shown using the solid line. (C) Contact maps between water molecules and Glu61 side chains. Presence of water molecules whose oxygen atoms are within 4 Å of Glu61 O ε1 and O ε2 atoms is indicated using gray (1 molecule) and black (2 molecules).
The molecular entities with RSCC values above 0.7 are generally recognized as consistent with corresponding maps (Neumann et al., 2018). Overall, the 6:9 system reproduces the densities in the acyl chain region very well, but clearly lacks needed density in the head group region.
After the CG simulations, we have also conducted short atomistic simulations of 6:9 and 7:10 systems. Both revealed similar correlation with the experimental map (Table 2), with the 6:9 system severely lacking densities in the head group region, and the 7:10 system having better correspondence there ( Figure 3B). The reason for the discrepancy in the head group region might be that we used the model lipid POPC in our simulations, whereas the natural thylakoid membranes contain a complex mixture of phosphatidylglycerol, monogalactosyldiacylglycerol, digalactosyldiacylglycerol and sulfoquinovosyldiacylglycerol lipids (Webb and Green, 1991). Phosphate and choline moieties of POPC have joint molecular weight of about 200 Da, whereas one galactose moiety weighs around 180 Da, two galactoses weigh ∼360 Da, and sulfoquinovose weighs ∼250 Da. Therefore, native lipids have head groups with higher molecular mass and bigger volume, and this might account for the additional density observed in the lipid head group region of experimentally determined density maps.
To gain insight into atomistic details and to probe the system stability, we have simulated the 7:10 system without any applied restraints for 500 ns. Overall, the system was stable with the protein C α RMSD of ∼1.1 Å ( Figure 1E) and the lipids remaining close to the original position. We observe a reasonable correspondence of modeled lipid positions to the experimental densities ( Figure 4A). Interestingly, the lipids encapsulated within the c-ring are displaced relative to the surrounding bilayer (Figure 4B), as it was observed before in experimental and modeled structures of other c-rings (Murata et al., 2005;Krah et al., 2010). At the loop side of the c subunit, the lipids are straightened and immobile. At the NC side, the lipids are more disordered and dynamic, forming a conical shape. Oleoyl chains of three POPC molecules at the NC side are bent so that some acyl chains are oriented along the membrane plane and perpendicular to other palmitoyl and oleoyl moieties. Finally, while no water molecules were initially placed within the membrane, we observed diffusion of several molecules into the membrane toward the active site glutamate (Figures 4B,C).
Based on the reported simulations and comparison of density distributions in the acyl chain region to the experimental data (Figure 3 and Table 2), we conclude that the experimental maps are best described by systems with 15 to 17 lipid molecules. These numbers are in contrast with 23 to 26 lipids observed in biasing potential-driven assembly simulations. Is it possible that trapping of the excessive number of lipids is a consequence of using a biasing potential? Our analyses indicate that this likely might not the case. First, most of the assemblies with the biasing potential coefficient below 0.01 kJ/(mol·nm 2 ) do not converge in reasonable time (Table 1), and thus the biasing potentials are not too excessive and affect the system moderately. Second, some lipids are known to interact tightly with corresponding membrane proteins, and such interaction might prevent reaching equilibrium number of lipids inside the c-ring in accelerated assembly simulation. However, analysis of the residence times of POPC lipids in the inner volume of c 7 half-rings separated by 5 nm reveals that lipids in most cases do not stay within the cavity for longer than ∼100 ns (Figure 5). Given that the typical assembly takes 0.5-1.5 µs, the lipids have an opportunity to leave the cavity during the assembly simulation. Third, lipids are known to co-diffuse with membrane proteins (Niemelä et al., 2010), and thus might be dragged by the half-rings pulled by the biasing potential. In assembly simulations with biasing potential coefficients of 0.01 kJ/(mol·nm 2 ) most of the lipids that are eventually trapped inside the c-ring are initially positioned in between the c 7 halfrings ( Figure 6A). On the contrary, in the 0.002 kJ/(mol·nm 2 ) assembly simulation, the lipids are evidently free to diffuse inside and outside the space between the approaching half-rings, and most of the lipids that are eventually encapsulated are initially positioned outside of the starting positions of c 7 half-rings (Figure 6B), yet still 23 lipids are trapped. Finally, assembly simulations with dioleoyl-phosphatidylcholine (DOPC), a lipid with two unsaturated, and strongly disordered, oleoyl chains, and dipalmitoyl-phosphatidylcholine (DPPC), a lipid with two saturated, and ordered, acyl chains result in very similar numbers of trapped lipids ( Table 1). To sum up, it appears that the effects arising due to using a biasing potential, if they are present, are moderate and do not explain the difference between the numbers of 23-26 lipids trapped in the assembly simulations and 15-17 lipids observed to fit best the experimental maps. Therefore, it is tempting to speculate that during the natural assembly process the c-ring might initially include a larger number of lipids, and then consequently some lipids are extruded, as depicted in Figure 2, and either can diffuse away themselves or might be carried away by some soluble transport proteins.
Preliminary data shows that the presented approach works well also for larger K 10 (40 TM helices) and smaller c 13 (26 TM helices) rings (data not shown). We believe that it can be used for assembly of other protein complexes with encapsulated lipids such as the multidrug exporter AcrB (Qiu et al., 2018) or oligomeric light-driven pumps, microbial rhodopsins (Gushchin and Gordeliy, 2018;Bratanov et al., 2019;Kovalev et al., 2019) as well. As lipids are often critical for correct folding and operation of membrane proteins (Cournia et al., 2015;Hedger and Sansom, 2016;Gupta et al., 2017), reproducing their correct positions in simulations is of utmost importance.

Preparation of Starting Models
Starting coordinates for the spinach c-ring were taken from PDB ID 6FKF. The symmetry axis of the ring was oriented along the z axis. For assembly simulations, half-ring oligomers were separated by 5 nm along the x axis. Atomic structures were converted into Martini representation using martinize.py (http://cgmartini.nl/index.php/tools2/proteins-and-bilayers/ 204-martinize). CG lipids, ions and water particles were added using insane (Wassenaar et al., 2015). The resulting CG systems contained 1308 POPC molecules and 8, 15 and 25776 Na + , Cl − and water particles, respectively, in a rectangular unit cell with dimensions of ∼24 × 18 × 12 nm 3 . For atomistic simulations, the systems were reassembled using the internal lipid positions taken from CG assembly simulations (Wassenaar et al., 2014), the experimental c-ring structure, and the external lipid positions generated using CHARMM-GUI (Jo et al., 2007). The resulting AA system contained 17 lipids inside the c-ring, 242 lipids outside of the c-ring, and 7 Cl − and 28348 water molecules, in a hexagonal unit cell with dimensions of ∼11.5 × 11.5 × 13 nm 3 . All simulations were performed using periodic boundary conditions. pKa values of titratable amino acids were determined using PROPKA3 (Olsson et al., 2011). Active site Glu61 were protonated. Glu37 of even chains were protonated because of borderline pKa value of ∼6.8. Other titratable amino acids and N-and C-termini were charged. All simulations were performed using GROMACS 5.1 (Abraham et al., 2015).

General Details of the Simulations
CG and AA simulations were conducted using leapfrog integrator with time steps of 20 and 2 fs, at a reference temperature of 350 and 303 K, respectively, and at a reference pressure of 1 bar. Protein, lipid and water molecules were coupled to the temperature bath separately. Temperature was coupled using velocity rescale (Bussi et al., 2007) and Nosé-Hoover (Nosé, 1984) thermostats with coupling constant of 1 ps −1 , respectively. Pressure was coupled with semiisotropic Parrinello-Rahman barostat (Parrinello and Rahman, 1981) with relaxation time of 12 or 5 ps, respectively.

Details of CG Simulations
CG simulations were conducted using Martini 2.2 (Marrink and Tieleman, 2013) and ELNEDYN (Periole et al., 2009) force fields. The center of mass of the reference structure was scaled with the scaling matrix of the pressure coupling. The non-bonded pair list was updated every 10 steps with cutoff of 1.1 nm. Potentials shifted to zero at the cutoff of 1.1 nm and a reactionfield potential with ε rf = ∞ were used for treatment of the van der Waals and electrostatics interaction as recommended (de Jong et al., 2016).

Details of AA Simulations
AA simulations were conducted using CHARMM36 (Klauda et al., 2010;Best et al., 2012) force field. The covalent bonds to hydrogen were constrained using SHAKE algorithm (Ryckaert et al., 1977). The non-bonded pair list was updated every 20 steps with cutoff of 1.2 nm. Force-based switching function with the switching range of 1.0-1.2 nm and particle mesh Ewald (PME) method with 0.12 nm Fourier grid spacing and 1.2 nm cutoff were used for treatment of the van der Waals and electrostatics interactions.

Details of RSCC Calculations
RSCC values have been calculated using Chimera (Pettersen et al., 2004). The experimental map sampling density is 1.053 Å in all directions, as used in the original report (Hahn et al., 2018). The simulated maps were obtained using the volmap utility in VMD (Humphrey et al., 1996) with the sampling density of 1 Å. The symmetry axes were determined using Chimera, and the maps were averaged with their symmetry-related images rotated by increments of 360 • /14≈25.7 • . The map regions where the density of the simulated maps was above 0.8 (the volume inside the cring) have been used for RSCC calculation. For reference, the maps shown in Figures 3D,E are drawn at the level of 1.3.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
ON and PB conducted the molecular dynamics simulations, analyzed the results and helped with manuscript preparation. IG designed and supervised the project, analyzed the results, and prepared the manuscript.

FUNDING
This study was funded by the Russian Foundation for Basic Research according to the research project No. 18-34-00986. We are grateful to the Center for Scientific Computing (CSC-IT Center for Science, Espoo, Finland) for computational resources.