Affinity of rhodopsin to raft enables the aligned oligomer formation from dimers: Coarse-grained molecular dynamics simulation of disk membranes

The visual photopigment protein rhodopsin (Rh) is a typical G protein-coupled receptor (GPCR) that initiates the phototransduction cascade in retinal disk membrane of rod-photoreceptor cells. Rh molecule has a tendency to form dimer, and the dimer tends to form rows, which is suggested to heighten phototransduction efficiency in single-photon regime. In addition, the dimerization confers Rh an affinity for lipid raft, i.e. raftophilicity. However, the mechanism by which Rh-dimer raftophilicity contributes to the organization of the higher order structure remains unknown. In this study, we performed coarse-grained molecular dynamics simulations of a disk membrane model containing unsaturated lipids, saturated lipids with cholesterol, and Rh-dimers. We described the Rh-dimers by two-dimensional particle populations where the palmitoyl moieties of each Rh exhibits raftophilicity. We simulated the structuring of Rh in a disk for two types of Rh-dimer, i.e., the most and second most stable Rh dimers, which exposes the raftophilic regions at the dimerization-interface (H1/H8 dimer) and two edges away from the interface (H4/H5 dimer), respectively. Our simulations revealed that only the H1/H8 dimer could form a row structure. A small number of raftophilic lipids recruited to and intercalated in a narrow space between H1/H8 dimers stabilize the side-by-side interaction between dimers in a row. Our results implicate that the nano-sized lipid raft domains act as a “glue” to organize the long row structures of Rh-dimers.


Introduction
The visual pigment rhodopsin (Rh) initiates the phototransduction cascade in vertebrate disk membranes of rod photoreceptor cells. Rh is a prototypical seven-transmembrane G proteincoupled receptor (GPCR) and is highly concentrated in the disk membrane, occupying approximately 30% of the total disk membrane area [1,2]. PLOS  suggested that approximately 70-80% of Rh molecules form dimers in the disk membrane [13,18,25]. Furthermore, Rh-dimers are known to locate in the disk membranes where saturated and unsaturated lipids are essential components [13,18]. Therefore, for simplicity we focused on the dynamics of these two types of lipids and Rh-dimers in the current model in order to advance our insight into the mechanisms of row structure formation by Rh-dimers. The disk membrane is a closed lipid bilayer membrane incorporating a tremendous amount of Rh. The dominant movement of Rh and lipids is 2-D diffusion along the membrane plane. The lipid-lipid and lipid-Rh positional exchanges along the membrane plane are the predominant limiting processes. Again for simplicity, we constructed the disk membrane  [30]. The ribbon model of H1/H8 dimer viewed from the eighth helix (H8) side. The yellow alpha-helix indicates H8 in each Rh-monomer, which was assumed to be palmitoylated. The red alpha-helices indicate H4 and H5. (b-c) Basic two-dimensional (2-D) structure models. (b) H1/H8 dimer 2-D structure constructed based on the Rh-dimer contour profile in (a). (c) H4/H5 dimer 2-D structure. The yellow particles indicate regions of H8 assumed to be raftophilic. The red particles are assumed to be parts of H4 and H5. The other particles (blue or green) represent other regions of the Rh dimer. (d) Schematic illustrations of unsaturated and saturated lipids (lower figures) and their representative 2-D particle models (upper figures). The affinity between particles of saturated lipids and the yellow particles in the Rh-dimer 2-D model was assumed (see panels b and c). The tails of each saturated lipid and cholesterol were assumed to frequently associate with each other, making the saturated lipids rigid compared to the unsaturated lipids.
https://doi.org/10.1371/journal.pone.0226123.g001 model using a complete 2-D particles system that described the structures and dynamics of lipids in one monolayer and Rh-dimers on the plane using circular particles (Fig 1). The influences of water and lipids in the opposing disk membrane monolayer were modeled as drag forces and noise.
The 2-D structure model of each Rh-dimer was constructed as a 2-D elastic network using 16 circular particles. For the H1/H8 dimer model, the particles occupied positions to imitate the contour profile of the basic Rh-dimer structure viewed from the H8 side based on X-ray crystal structure analysis of the activated form of Rh-dimers (PDB ID: 3CAP) [30] (Fig 1A and  1B and S1A Fig). The radii of the particles were assumed to be similar to those of alpha helices and the nearby particles were connected by springs where the natural length of each spring connecting two particles was equal to the distance between them according to the basic 2-D Rh-dimer structure model (S1B Fig and S1 Table). As there are no reference structures based on experiments such as X-ray crystal structure analysis, the H4/H5-dimer model was constructed by exchanging the positions of the particles on the left half side of the H1/H8 dimer model corresponding to one Rh-monomer with those on the opposite side ( Fig 1C, S1C Fig,  S1D Fig and S2 Table). We note that the results presented in this paper were independent of the detailed shape of the 2-D model of Rh-dimers and were qualitatively unchanged.
The structures and physical properties of the unsaturated and saturated lipids in the 2-D model were as described. Both types of lipids were assumed as circular particles. The radii of the circular particles were assumed as those when lipids are approximated as cylinders ( Fig  1D) = 0.43 nm, which is comparable to the scale of the section of the phospholipid glycerolhead estimated by recent studies [31,32]. No specific attractive interactions were assumed among the lipids. Notably, disk membranes have been known to contain sufficient concentration of cholesterol compared to that of saturated lipids; the frequent interactions of cholesterols with saturated lipids are expected to make the saturated lipids rigid and raftophilic compared to that of unsaturated lipids, where the stoichiometric ratio between total lipids and saturated lipids is~100 : 8, while that between total lipids and cholesterol is~100 : 11, in disk membrane [1]. Therefore, we assumed the cholesterols were always associated with saturated lipids and the lipids were always more rigid than unsaturated lipids. These assumptions resulted in the saturated lipid membranes being more rigid and less elastic than that of unsaturated lipid membranes, which is consistent with experimental findings [33]. Furthermore, as mentioned below, the simulations under these assumptions demonstrated the following results, i) the diffusion of saturated lipids was slower than that of unsaturated lipids and ii) the saturated lipid domains that may correspond to the ordered lipid raft domains appeared through the phase separation between saturated and unsaturated lipids (see the "Parameters for simulations" section in the Results). These results were consistent with previously reported experimental findings [29,34,35].
Additionally, two H8 helices located near the center of the H1/H8 dimer and near the edges of the H4/H5 dimer were expected to be partially raftophilic since acyl chains in H8 are known to be palmitoylated [4,12,36]. Therefore, in the models of these dimers, we assumed an affinity between saturated lipids and the tip of H8 (yellow particles in Fig 1C), except for some modified models as detailed below.

Model implementation
In the model used in this study, we described the central region of the disk membrane to consist of Rh and two types of lipids using 2-D particles moving on a 2-D plane. Since the lipids and Rh would be subjected to water and lipids in the opposing monolayer of the membrane, we assumed that the motion of each particle, the lipids, or parts of Rh obeyed overdamped Langevin equation, as follows: where x i = (x i ,y i ) was the position of the i-th particle and V indicated potential of forces working on particles. The γ i and R i (t) indicated the coefficient of drag force and the random force working on i-th particle by water and lipids in the opposite monolayer, respectively. The term R i (t) was given as Gaussian white noise and satisfied hR i (t)i = 0, and hR i (t)R j (s)i = 2γ i k B Tδ ij δ(t −s) where k B was the Boltzmann constant, T was the temperature, δ ij indicated the Kronecker delta, and δ() indicated the Dirac delta function. The first term of the right-hand side of Eq (1) indicated interactions among particles provided by the potential of the system V according to the following equation: where V collision was the potential of the excluded volume effects among the particles, V bond was the interaction potential among the particles forming Rh-dimers to sustain the shape of each dimer, and V raft was the affinity for the interactions between Rh dimers and saturated lipids. The potential of the excluded volume effects among particles (V collision ) was denoted by the following: where k c ij ¼ kq i q j was the elastic constant when the i-th and j-th particles contacted each other, q i indicated the nondimensional parameter of rigidity of the i-th particle, r i indicated the radius of the i-th particle, and θ was the Heaviside step function defined by the following: ( We assumed q i depended on the type of particles, where q i of saturated lipids was assumed to be appropriately larger than those of unsaturated lipids. In this case, the lipid type-dependent diffusion constant and the saturated-unsaturated lipid phase separations were consistently obtained based on recent experiments as described below in the Results section and as previously described [29,34,35]. The variables r h and r l were the radii of cylinders, which approximated an alpha-helix of Rh and the lipids in the disk membrane, respectively. We assumed r i = r h when the i-th particle was a component of Rh and r i = r l when the i-th particle was a lipid. The term ∑ indicated the sum of the i-th and j-th particles that did not belong to the same Rh-dimer.
The interaction potential among the particles forming Rh-dimers was defined as V bond , and was calculated as follows: Where d Rh ij was the distance between the i-th and j-th particles of the basic 2-D Rh-dimer structure (S1 Fig and S1 and S2 Tables) and k b was the elastic constant for sustaining the shape of the basic structure of each Rh-dimer. The ∑ in this equation indicated the sum of the i-th and j-th particles that belonged to the same Rh-dimer and that were in close proximity to each other (S1 Fig).
The interaction potential among the palmitoylated H8s of the Rh dimers and the saturated lipids was indicated as V Raft and was calculated as follows: where ε raft indicated the affinity between a palmitoylated H8 and a saturated lipid. The S indicated the sum of particles pairs of a particle occupied in H8 of Rh and a saturated lipid (Fig 1).

Simulation parameters
As far as possible, the parameters of the model employed in the simulations of the current study were based on experimental findings. The parameters included radii of lipids (r l ) and alpha helix (r h ) of Rh and were assumed to be 0.43 nm and 0.88 nm, respectively. To simulate molecular dynamics similar to the situation at the central region of the disk membrane, we used a 40 nm × 40 nm square box with periodic boundary conditions, which included 14 Rh-dimers (14 × 2 × 8 = 224 helices), 156 saturated lipids, and 1,704 unsaturated lipids as entire the simulation space. The area occupancy rates of the Rh-dimer, saturated lipids, and unsaturated lipids in the present model were given as approximately 27%, 5%, and 58%, respectively, which matched the rates for the entire disk membrane obtained from previous experiments [1,2]. The parameters for the excluded volume effect of each particle (k and q i ) were assumed according to k/k B T = 7.5 (nm −2 ), q i = 1.5 for unsaturated lipids, q i = 4 for saturated lipids, and q i = 9 for particles in Rh. To sufficiently sustain the Rh-dimer structure, we assigned k b /k B T~18870 (nm −2 ). The drag coefficient γ i was assumed at γ i = 6πηr i . In the current model, η was expected to be relatively large compared to cytoplasm as we did not know the precise value from previous experiments because of the influences of the lipids at the opposing lipid bilayer. Thus, we assumed η/k B T = 6×10 −6 (mm −3 s) in the current arguments, while η/k B T~1.55×10 −7 (mm −3 s) is estimated in cytoplasm [37]. In this case, the additional simulation of the model using these values for the parameters without any affinity among the lipids and Rh-dimer (such a dimer was termed a raftophobic H1/H8 dimer and is defined in the next section), we obtained diffusion coefficients for unsaturated lipids (in an unsaturated lipid domain), saturated lipids (in a saturated lipid domain), and Rh-dimer (in an unsaturated lipid domain) of approximately 27 (μm 2 / s), 3.0 (μm 2 /s), and 0.5 (μm 2 /s), respectively. These values were estimated by mean square displacement of the particles (S2 Fig and S3 Table). The diffusion coefficients were consistent with those observed in recent in vitro experiments using artificial model membranes [11,13,38,39].
The parameter used to estimate the affinity between the particle at H8 of the Rh-dimer and the saturated lipid was assumed according to ε raft /k B T = 62.5, which was based on the results of the following additional simulation of a model membrane containing only one Rh-dimer and an equal number of saturated and unsaturated lipids with the abovementioned parameters (S3 Fig). Using this simulation, we obtained the probability of contact between the Rh-dimer and saturated lipids as approximately 0.627. This result seemed consistent with our in vitro study using artificial lipid bilayers containing saturated and unsaturated lipid domain where approximately 62% of the Rh-dimers with H8 palmitoylation were found to have saturated lipid domains at equilibrium (Tanimoto, et al., submitted). However, we note that the detailed values ε raft /k B T were not essential and similar contact probability values could be obtained when ε raft /k B T�1.

Analysis methodology
To simulate the present model, the time integral of Langevin Eq (1) was calculated numerically using the Eular-Maruyama method with unit MD step = 0.1ps [40] (S1 File). To estimate the degree of order of the spatial distribution of Rh-dimers, we defined the degree of row structure S(t) as follows: where the position of the center of mass of the m-th Rh-dimer (R m and θ mn ) and the angle between the orientation of the m-th and n-th Rh dimers as shown in S4 Fig. The contribution of m-th and n-th Rh-dimers to S(t) was large when the centers of masses were close to each other and when they were facing the same direction. The value of cos(θ mn ) was determined from the inner product between the vector from one Rh-Rh interface particle to another of the m-th Rh-dimer and that of the n-th Rh-dimer (S4 Fig), where 2(cos 2 (θ mn )−1) was used as the degree of plane orientation in liquid crystal [41]. Since the distance between Rh-dimers observed in cryo-electron tomography is approximately 5.5 nm (Gunkel), we assumed S(t) drastically decreased when the distance between two Rh-dimers was greater than 5.5 nm. Ideally, S(t) = 0 when Rh-dimers were distributed randomly.

H1/H8 dimers could form row structures of Rh-dimers by their raftophilicities
The model simulations of disk membrane containing H1/H8 dimers were performed in the current study. We assumed all molecules were randomly distributed under the initial conditions (Fig 2A). Results of the simulations after extended periods of time from the initial condition revealed that the model with H1/H8 dimers exhibited some row structures of Rh-dimers (Fig 2A, S1 Movie), consistent with findings observed in AFM and cryo-electron tomography [15][16][17][18]. On the other hand, no Rh-dimer row structures were observed in models using raftophobic dimers (Fig 2B), where raftophobic H1/H8 dimers indicates H1/H8 dimers but assumed ε raft = 0. According to the measurement of S(t) and their averages over 10 simulations (hS(t)i), the model using H1/H8 dimers exhibited larger S(t) and hS(t)i values and greater averages after extended time from the initial conditions compared to those using raftophobic dimers (Figs 2D and 3).

H4/H5 dimers accumulated but failed to form row structures of rhodopsin (Rh)-dimers
The simulation of the model containing H4/H5 dimers showed Rh-dimer accumulation with some branching structures, but no row structures were observed (Figs 2C, 2D and 3).
In models using either dimers, H1/H8 or H4/H5, associations among the Rh-dimers occurred through the small regions of saturated lipid accumulation near the areas corresponding to raftophilic H8s of each Rh-dimer. Notably, for H1/H8 dimers, two dimers associated with each other when the saturated lipid domains formed at the regions between two the central portions of these dimers. Each domain connected only to the central regions of two Rh dimers that corresponded to their raftophilic H8s. This fact enabled the H1/H8 dimers to form long rows of Rh-dimer structures (Figs 2B and 4A); the lengths of the rows were dependent determined using the number of lipid rafts formed in the membrane. In contrast, H4/H5 dimers could associate with more than two dimers through the domain of saturated lipids that formed at the edges of each dimer. This made it difficult for H4/H5 dimers to form row structures.

Saturated lipid domains stabilized the row structure of Rh-dimers
Saturated lipid domains were formed due to the affinity of the saturated lipids to raftophilic H8 regions of dimers and the phase separation of saturated and unsaturated lipids (S5 Fig). The stability of the row structures of H1/H8 dimers was influenced by the size of the saturated lipids domains (i.e., the number of saturated lipids in each domain) between two H8 regions of neighboring dimer pairs. The contact between the two Rh-dimers became unstable with decreases in domain size. We determined that hS(t)i of the model system containing only two H1/H8 dimers decreased with the decrease in the number of saturated lipids (N s ) between the two dimers when N s was less than 5 (Fig 4 and S6 Fig). However, hS(t)i t for all simulation results (results of ten simulations) exhibited values near the mean when N s �6, even though hS(t)i t was often near 0 when N s �5. This finding indicated that lipid domains containing more than 5 or 6 saturated lipids were essential forming and stabilizing Rh-dimer row structures, not single saturated lipid molecules.

Discussion
In this study, we developed a coarse-grained model of retinal disk membrane consisting of saturated lipids, unsaturated lipids, and Rh-dimers used to simulated row structure formation of Rh-dimers observed in recent experiments. Using this model, we clarified that H1/H8 dimer is the primary element constructing the row structures of Rh dimers. As far as we know, the present study was the first to provide a scenario to explain the mechanism of formation and stabilization of row structure by Rh-dimers in the disk membrane.
Recent experimental evidence have suggested that more than 70% of Rh is in dimeric or higher oligomeric state(s) existing in a dynamic equilibrium [13,18,25]. Based on these findings and our current results, the H1/H8 dimer forms transient row structure through lipid raft-based interaction between the dimers, which appears consistent with results from recent studies [13,42,43]. In addition, we found that the nano-sized domains of saturated lipids with cholesterols (raftophilic lipids, i.e., lipid rafts) were also key factors in forming Rh-dimer row structures with H1/H8 dimers. Such domains acted as "glue" to connect two raftophilic H8 regions of two Rh-dimers to construct the row structures.
Due to the hybridization in the presented results and the recent idea of a dynamic scaffolding mechanism for the Rh-G t (G protein transducin) interaction [27], the dynamics of G t on the disk membrane were predicted as follows. Our in vitro experiments suggested that G t was basically raftophobic (Tanimoto et al., submitted). However, the volume fraction of the raft forming saturated lipids on disk membrane is only~5%, and our simulation results suggested almost all raft lipids were confined between each pair of dimers, stabilizing the structure of the rows of Rh-dimers. Therefore, we expected G t to diffuse freely on the inter-Rh-dimers rows region with a few saturated lipids, or to prebind to Rh-dimers and hop on the regions of stable Rh-dimers rows, similar to the previously proposed model [17,20,27], without any limitations caused by lipid raft domains (Fig 5).
In our model, the saturated lipids were assumed to always associate with cholesterols, which allowed the saturated lipids to take rigid structures compared to that of unsaturated lipids. As a result of this assumption, the lipid raft could form by phase separation between saturated and unsaturated lipids. It should be noted that our simulation strongly suggested that lipid raft formation was essential for the formation of the row of Rh dimers, but did not suggest much about the molecular mechanism of lipid raft formation. In order to reveal the detailed mechanism of lipid raft domain formations, we need to consider various other physicochemical properties of lipids, such as the electrostatic properties of the polar head [44,45], as well as, the effects of cholesterols by the considering molecular models with microscopic details; this is one of our future challenges. On the other hand, we expect to obtain results that are similar to Coarse-grained molecular dynamics simulation of rhodopsin aligned oligomer formation the ones presented, if such physicochemical properties of lipids play weak or positive roles in raft formation.
In this study, we carefully chose the parameters of the model. However, we found that Rhdimers could form the row structures if the ratio of rigidity (q i ) of saturated and unsaturated lipids and the raftophilicity of Rh-dimers were sufficiently large. This fact provides additional evidence confirming the abovementioned fact that only the H1/H8 dimers and saturated lipid domain formations were essential for Rh-dimer row formation.
In our model, the model of H1/H8 dimer was constructed based on the crystal structure of the activated form of the Rh-dimer [30]. Recently, photoactivation of Rh was known to induce a structural rearrangement, for example, an increase in the distance between specific helices of Rh [43,46,47]. However, such local deformations did not affect the global geometric shape of Rh-dimer, and the structures around H1/H8 interface were highly conserved. Therefore, we can construct the model of the deactivated form of the Rh-dimer using the present model with a few modifications and obtain almost the same results.
The current model was unable to reproduce the repetitions of formation and collapse of row structures of Rh-dimers that have 100 ms range of life time [13,[48][49][50] because of difficulties such as simulation costs for example. The limitation and difficulties will be addressed in future studies by modifying the model to allow for accelerated calculation.

S1 Fig. Detailed construction of a two-dimensional (2-D) elastic network model of rhodopsin (Rh) H1/H8 and H4/H5 dimers. (a)
Relationship between H1/H8 dimer structure obtained by X-ray crystal structure analysis (ribbon model) seen from H8 side (PDB ID: 3CAP) and particle positions of the 2-D model of H1/H8 dimer (black rings). (b) Constructions of 2-D models of H1/H8 dimer consisting of 16 particles. Each particle was named 1st through 16-th and assumed to connect with nearby particles by springs where the natural length of the spring connecting the i-th and j-th particles (d Rh ij ; S1 Table) was given to maintain the basic structure of Rh-dimer. The 8th and 16th particles (yellow) indicate the particles corresponding to the palmitoylated areas of H8 and 3rd and 11th particles (purple) indicate the Rh-Rh interface particles. Blue and purple curves indicate examples of connection between the 4th and others particles and the 7th and other particles. The d Rh ij for each pair of basic structures of Rh-dimers are shown in S1 Table. (c) Relationship between the expected H4/H5 dimer structure and particle positions of the 2-D model of the H4/H5 dimer. The expected H4/H5 dimer structure (ribbon model) was artificially constructed from the X-ray crystal structure of the H1/H8 dimer. The positions on the left and right Rh-monomers in H1/H8 dimer were interchanged. (d) Constructions of the 2-D models of the H4/H5 dimer consisting of 16 particles with the natural length of the spring connecting i-th and j-th particles (d Rh ij ; S2 Table). (a) Typical initial state of the simulations to determine the appropriate value of ε raft . One Rh was inserted with the saturated and unsaturated lipids being randomly distributed. A total of 340 saturated lipids and 370 unsaturated lipids were confined to a 24 μm × 24 μm square box with periodic boundary conditions. Since unsaturated particles can overlap with each other more easily than saturated lipids, the area fraction of saturated and unsaturated lipids obtained were approximately the same, even though the absolute number of the two types of lipids differed. (b) Typical snapshot of model simulations from the abovementioned initial conditions. The vicinity of the Rh area was defined by the area surrounded by the white curve, which was closer than 2.4 nm from each particle. (c) Ratio of saturated lipids around the Rh-dimer defined by r ¼ N s N s þN u , predicted by experimental findings [Y. Tanimoto, submitted] and that of various raftophilicity and unsaturated lipid (cyan particles) was set at approximately 8% (80 particles) and 92% (920 particles), respectively. When q saturated = q unsaturated = 1.5, ϕ was estimated to be 0.08 as the ratio of saturated was 8%. In cases of q saturated = 4 and 9, ϕ exhibited significantly larger values than that in the case of q saturated = 1.5. This suggested that due to their rigidity, saturated lipids tended to form domains that may correspond to "raft lipid domain". (b) Snapshots of saturated and unsaturated lipid configurations at q i = 1.5,4,9. The greater the q i value resulted in increased amounts of raft lipid domain that could be observed.