Molecular dynamic simulation of Ca 2+ -ATPase interacting with lipid bilayer membrane

: In biomedical and drug delivery treatments, protein Ca 2+ -ATPase in the lipid bilayer (plasma) membrane plays a key role by reducing multidrug resistance of the cancerous cells. The lipid bilayer membrane and the protein Ca 2+ -ATPase were simulated by utilising the Gromacs software and by applying the all-atom/united atom and coarse-grained models. The initial structure of Ca 2+ -ATPase was derived from X-ray diffraction and electron microscopy patterns and was placed in a simulated bilayer membrane of dipalmitoylphosphatidylcholine. The conformational changes were investigated by evaluating the root mean square deviation, root mean square ﬂ uctuation, order parameter, diffusion coef ﬁ cients, partial density, thickness and area per lipid.


Introduction
Ca 2+ -ATPase, adenosine triphosphate (ATP), which is a serine type lipid, is a transporter protein in the plasma membrane which acts as calcium (Ca 2+ ) remover from the cells and its action is vital for regulating the amount of Ca 2+ within the cells [1]. In fact, Ca 2+ -ATPase is involved in pumping Ca 2+ from all eukaryotic cells [2]. It is thought that the Ca 2+ -ATPase has 10 segments that cross the plasma lipid bilayer membrane, with both termini, namely C and N (shown in Fig. 1), remaining inside the cell. At the C terminus, there is a long 'tail' of between 70 and 200 amino acids in length [2]. This tail is responsible for regulating the pump operation [2]. The pump is powered by hydrolysis of ATP with a stoichiometry of one Ca 2+ ion removed for each molecule of hydrolysed ATP. ATP binds tightly to Ca 2+ ions but does not remove Ca 2+ at a very fast rate [3].
There is a very large transmembrane electrochemical gradient of Ca 2+ driving the entry of the ion into the cells, yet it is very important for the cells to maintain low concentrations of Ca 2+ for proper cell signalling; thus it is necessary for the cells to employ ion pumps to remove the Ca 2+ [4]. Since Ca 2+ -ATPase transports Ca 2+ into the extracellular space; it is also an important regulator of the calcium concentration in the extracellular space [5]. The interaction of membrane proteins with the lipids is essential for their function. Bilayer properties, such as lipid composition, can affect the membrane protein activity [6]. The ability of Ca 2+ to regulate both cell death and proliferation, combined with the potential for pharmacological modulation, offers the opportunity for a set of new drug targets in cancer treatment and Ca 2+ signalling to regulate cells' actions in cancer treatment may be considered as a therapeutic option [7]. Moreover, investigating these biophysical behaviours is important to explain how these molecular motors work.
It is very difficult to determine the structure of membrane proteins experimentally [8] and computational methods [9] have been increasingly used to study their structure and function. Therefore the problem that is dealt with in this study will be to understand the interaction of protein molecules and real cell membranes in the human body.
The Gromacs software [10] with 10 ns timescale was used to simulate four systems of membranes consisting of lipid (dipalmitoylphosphatidylcholine -DPPC) + protein (Ca 2+ -ATPase). In two systems (systems 1 and 2), united atom (UA) was applied to simulate the lipid bilayer membrane (DPPC) and the protein (Ca 2+ -ATPase) was simulated using all atom (AA) model. In two other systems (systems 3 and 4), coarse grained (CG) model was applied to both DPPC and protein. The characteristics of simulated systems are given in Table 1. The simulations started with two initial structures of the protein; proteins 1 and 2 and their structures were obtained, respectively, from X-ray and electron microscopy (EM) patterns. These structures are shown in Fig. 1.
The aim of this simulation study is to investigate the stability and structural changes of proteins 1 and 2 as well as to quantify the perturbations of the bilayer membrane in the presence of the protein. To achieve the aim of this research, the simulations were performed in the following steps: 2. Comparing the results of coarse-grained molecular dynamics (CGMD) and all-atom/united-atom molecular dynamics (AA/UA MD) simulations to assess the reliability of the CGMD method. 3. Comparing the obtained results from the two initial protein structures and examining the effect of protein initial structure on the accuracy of the final results.

Initial structures
In this work, we intend to study the protein Ca +2 ATPase interacting with a model DPPC membrane by molecular dynamics (MD) simulation. The model's components are Fig. 1 Simulations started with two initial structures of the protein; proteins 1 and 2 and their structures were obtained, respectively, from X-ray and electron microscopy patterns a Structure of the Ca 2+ -ATPase in the E1 conformation with bound calcium ions obtained from 1SU4 entry of the PDB [10]. A-domain (residues 1-43 and 124-235), N-domain (residues 360-600), and P-domain (residues 330-359 and 601-739) relative to the transmembrane helices TM1-TM2 (residues 44-123), TM3-TM6 (residues 239-329 and 740-821) and TM7-TM10 (residues 831-994) [11,12] b,c AA (left) and CG (right) structures of the Ca 2+ -ATPase (PMCA). Fig. 1b: X-ray structure and Fig. 1c: EM structure. Fig. 1b is a conformation where Ca 2+ is present (E1 state), whereas Fig. 1c is conformation corresponding to the Ca 2+ -free (E2 state). The protein is far less stable when the calcium ions are removed [13] represented in Table 1. Although the model systems have close similarity with the real membrane, they differ in some aspects from the natural environment of biological systems. However, designing a model membrane with natural environment, both theoretically and experimentally, is a complicated and difficult task and requires extensive, time-consuming computational MD simulation work. Then to keep the computational work at a practical and acceptable level, as is a usual procedure in the MD simulations, makes some simplifications on the applied models inevitable. However, the simplifications made on the membrane models used in this work are less than similar previous MD simulations and the models are much closer to the real natural environment. In what follows, we explain the implementation of the MD simulation to understand and analyse the mechanism of the proteinmembrane interaction for drug delivery applications and the reliability of the membrane models will be justified by comparison of the obtained simulation results with the experimental data.
Four simulations (the four systems in Table 1) were performed using Gromacs 4.5.4 package [10]. Two initial structures, as s h o w ni nF i g s .1b and c, obtained respectively from the X-ray diffraction and EM of the protein Ca 2+ -ATPase (Protein Data Bank -PDB entry of 1SU4 and 1KJU [13]), were utilised in the simulations. The DPPC in the AA MD simulation consisted of 384 lipid molecules and was derived from lipid book [14] and the DPPC in the CGMD simulation consisted of 2048 lipid molecules and was obtained from the Science and Technology Facilities Council (STFC) [15]. The protein-bilayer system was solvated in 85 000 to 105 000 water molecules and 27 positively charged ions Na + were added in order to maintain the electroneutrality of the system (Table 1). Using the steepest descent minimisation method, the final set up of the systems was obtained ( Table 1) which consisted of protein, membrane, ions and water with the z-direction of protein perpendicular to the lipid bilayer membrane. The UA and CG models were employed to reduce the computation time. In the UA membrane; only the interactions of non-polar hydrogen atoms were ignored. The UA membrane and AA protein were modelled using the GROMOS96 53a5 force field [11] and by obtaining the lipid parameters from [16]. Two other systems (systems 3 and 4, Table 1) were modelled by the MARTINI CG [17][18][19], where the atoms of each residue of protein and membrane were mapped into 2-5 beads. The protein CG PDB files were constructed as described in the MARTINI home page [20]. The protein was inserted into the pre-equilibrated, hydrated bilayers by applying inflateGRO script [21]. Briefly, the membrane was expanded four times and then a hole of approximately the same size as the protein was made in the lipid bilayer by removal of lipid molecules. After insertion of protein, the lipid was compressed gradually to its initial structure. This procedure took about 25 steps and at every step the energy minimisation was done using about 6000 steps (Fig. 2).

Simulation conditions
Water was modelled using the SPC/E model [22]. The Columbic interactions were calculated using the particle Table 1 Characteristics of MD simulation systems In systems 1 and 2, UA model for the membrane (DPPC) and AA for the protein were used. In systems 3 and 4, the membrane and protein were simulated by CG model. The structure of proteins 1 and 2 were, respectively, obtained from the X-ray and EM pattern mesh Ewald [23] approach and the van der Waals interactions were truncated at 2 nm. All bonds were constrained using the LINCS algorithm [24] and the equations of motion were integrated using the leap-frog algorithm [25][26][27] with a time step of 2 fs. Periodic boundary condition was applied in all directions (x, y, z) and the centre of mass motion of the system was displaced at every time step. The temperature of all simulated cells was fixed at 323 K by coupling of three Nose-Hoover thermostats [28,29] with the bath coupling constant of 0.5 ps. The box vectors were subjected to semi-isotropic pressure coupling using a Parrinello-Rahman barostat [30] with a reference pressure of 1 bar, a coupling parameter of 2 ps and an isothermal compressibility of 4.5 × 10 −5 bar. In each simulation the protein molecules were located inside the membrane and the whole system was solvated with a sufficient amount of water molecules (see Table 1). To neutralise the charged systems, adequate sodium ions were added to the system (see Table 1). After energy minimisation, the equilibrium step for NVT ensemble (constant number of particles, volume and temperature) and NPT ensemble (constant number of particles, pressure and temperature) were performed for 2 ns. The MD run duration for each simulation was 10 ns.

Root mean square deviation (RMSD)
The structural changes undergone by the proteins were investigated by computing the RMSD of the α-carbon atoms with respect to pre-equilibrated configuration. The RMSD for an N-atom system with the variance σ at time t is defined as where r i (0) is the reference coordinates. Usually r i (0) is taken from the first frame in a MD simulation, though it can also be taken from the PDB [13] or a target structure. RMSD is very useful in monitoring the approach to equilibrium, typically signalled by the appearance of a broad shoulder. It is good practice to keep monitoring RMSD during production runs to ensure that the system stays near equilibrium.
Our simulations showed good stability, for the RMSD reached a well-defined value in the whole simulation time. Having smaller RMSD, as has been observed in a previous simulation work [31], it is seen from Fig. 3 that system 4 with protein 2 is more stable than the systems with protein 1. Fig. 3 indicates that the proteins in systems 1, 2 and 4 are almost stabilised in 2 ns (2000 ps), but a sudden increase in RMSD for system 3 after 2 ns is observed which, as reported previously [32], may be because of the change in the protein structure. However, these anomalous results for system 3, may also depend on the stability condition of the protein initial structure [33]. For example, as mentioned above, shown in Fig. 3 and also reported in a previous work [13], protein 2 (the structure which includes Ca 2+ at E1 state) is more stable than protein 1 at E2 state. System 4 is the best converged of all and therefore CG protein 2 is more stable than AA protein 2. This behaviour is probably because of the presence of an extra energy in the AA systems caused by the application of detailed force fields which requires more time for stabilisation [32]. Except for system 3, which indicates a drift in RMSD in the range 2000-6000 ps (Fig. 3), the proteins exhibit relatively low RMSD values which remain almost constant during the simulation runs.
As expected, and seen in Fig. 3, the RMSD of CG-protein in system 4 approaches a constant value in a short time, which can be considered as an indication of formation of a stable conformational structure, whereas in system 3 with CG protein the approach to a constant value occurs in a much longer time. This behaviour can be attributed to the higher energy of system 3, which causes a continuation of conformational change in the protein initial structure and as a result retards the formation of a stable conformational structure. More simulation results that will be presented in the following sections will substantiate these justifications.

Root mean square fluctuation (RMSF)
The RMSF can be computed by the following equation where (t) is the simulation time and x i (t j ) is the coordinates of atom i at time t j . During the course of simulation, the sum of the squared difference of the x i (t j ) and the mean coordinate x i were calculated and then divided by simulation time t and the root was calculated. Hence, the fluctuation of an atom around its mean in the trajectory files was obtained over 10 ns production run for each simulation. Although the protein structure is globally stable in the course of simulations, this does not imply that no significant conformational fluctuations are taking place. Thus, it is of interest to examine the magnitude of the conformational fluctuations in different regions of the studied systems and compare the results.
In Fig. 4 qualitatively it is seen that the two AA curves and CG curves are well overlapped for each protein, for protein 1 in systems 1 and 3 and for protein 2 in systems 2 and 4, and the highest fluctuations belong to protein 1. It should be noted that the high RMSF values in Fig. 4 are observed for the membrane parts of the proteins which are namely: (TM1-TM2: residues 44-100), (TM3-TM6: residues 270-300) and  Fig. 1a.
Moreover, Fig. 4 shows that the MARTINI CG model is too flexible; the residues in the CG model (system 3) fluctuate much more than the residues represented in the AA model (system 1). This trend for protein 2 is unlike protein 1 and shows that protein 2 in the CG-model is more stable and less flexible. However, although the CG protein 1 indicates more flexibility than the AA protein 1, the same patterns as the AA system's highest peaks have been observed, which corresponded to the large fluctuations of the residues 430-600 and 370-380 (Fig. 4). As stated before, the structure of protein 1 in system 3 is less stable. This is due to conformational change of the protein during the simulation, which may affect the residues' dynamics and as a result causes large differences in the calculated values of RMSF for protein in CG and AA models.

Diffusion coefficient
To evaluate the diffusion coefficient, the mean square displacement (MSD) of the protein was calculated by Einstein's equation in the following form where r (t 0 ) and r (t+t 0 ) are, respectively, the positions of the protein at time t 0 and t+t 0 , and the angle brackets represent their mean square deviation at time t 0 . (Table 2) As seen in Fig. 5 for the CG systems, the larger values of MSD are due to the CG force fields which map the atoms into the sites, and as a result some of the degrees of freedom are reduced. Moreover, the large value of MSD may be due to high-frequency intra-molecular vibrations that incorporate into an averaged effective interaction between the sites.
In the CG models, proteins 1 and 2 have almost the same value of diffusion coefficients and the MSD variations, as shown in Fig. 5, have the same trend of variation up to 8000 ps and then deviate from each other. This deviation comes from the change in the structure of protein, which in turn affects the lipid structure. However, the UA-model membrane (UA model) with protein 1 (system 1) has a better match with the experimental diffusion coefficient. This can be attributed to it having less constraint dynamics compared with protein 2. The decrease in diffusion coefficient in the UA membrane, as has been mentioned by other researchers [36], may be reflected by protein existence in the membrane. These results indicate that by increasing the concentration of protein, a decrease in overall mobility of the membrane occurs, which causes a reduction in the lipid diffusion coefficients. In comparison with the pure membrane, the existence of protein could reduce lipid diffusion coefficients too. However, the reduction in lipid diffusion coefficient in the system with protein 2 is greater than the system with protein 1 and as mentioned earlier the initial structural difference between these two proteins can be the main cause of this difference.

Area per lipid
The average area per lipid can be evaluated by multiplying the two dimensions (x, y) of the simulation box and subtracting the area of protein and then dividing the result by the number of lipid molecules in one leaflet of the lipid bilayer membrane (for the number of lipid molecules used in the systems, see Table 1).
The average saturated area per lipid at T=323 K for a DPPC bilayer is presented in Fig. 6. It is seen that all membrane systems become stable after about 5000 ps for UA models and after about 1000 ps for CG models. This indicates that CG models need less time for stabilisation. However, as there are more force fields acting between the atoms in the UA models, they need longer time for stabilisation. Table 3 represents the Error % of the calculated area per lipid membrane compared with the experimental data. It is seen that the area per lipid for the UA-model membranes are almost matched (Error <10%) with experimental data and in the case of CG-model membranes, the agreements are better (<2%). This means that the CG models are more stable because they need less time for stabilisation.

Partial densities and thickness
To obtain the thickness of membrane in the simulation systems, in Fig. 7 the partial densities of the head groups of the lipid are plotted against the Z-axis of the simulation  www.ietdl.org box. The bilayer thickness can obtained from the distance between the head group peaks in the density profile, as shown in Fig. 7. It is seen in Fig. 7a that for the lipid tail in the density profiles there is only one peak, but a few small peaks which are due to the distance between the tails of the lipid bilayer are seen in the middle area. As shown in Fig. 7a, starting from the outside of the bilayer membrane, an increase in the head group density is observed and after rising to a peak, the density of head group significantly drops down. This trend continues until reaching the peak for the tail groups, where the density of head group is zero and then starts to rise again. Moreover, this figure indicates that two curves for CG models are almost matched and the values of density, as expected, are slightly higher than those of UA models. It should be that, since different box size for UA models was used, these curves are separated from each other in the x-axis, but this separation does not occur in CG models with the same box size.
By considering Table 4, it is seen that the thicknesses in the membranes are close to the experimental data and this reflects the relative stability of the membranes in the studied systems. However, inspecting the reductions in the membrane thickness (as represented in Table 4) and their agreement with previous reports on the presence of protein in the membranes which reduced the membrane thickness [37,38] and by considering the fact that the smallest membrane thickness, in Table 4, belongs to system 4 which is the most stable system, the relation between stability of membrane and its thickness becomes evident. That is, the membrane in system 1 is less stable, because as expected, protein 1 in AA model is more dynamic and this point is confirmed by its higher diffusion coefficient (as reported in Table 2).
It should be noted that TM α-helices are embedded in the lipid bilayer, and as a result the change in the structure of the lipid bilayer leads to change in the function of the protein ATPase and vice versa. Moreover, the release of Ca 2+ ions from the ATPase involves a significant change in the packing of the TM α-helices [40] in the ATPase structure (as shown in Fig. 1) and this can explain why the protein in system 3 has the highest fluctuation as indicated in Fig. 3.I t is evident that because of close packing of ATPase, the lipid bilayer membrane in system 3 has the lowest thickness.
The Z values of the lipid head groups in the XY-plane, which indicate the thickness of the lipid bilayer membrane [37], are shown in Figs. 7b and c. The Z value of the lipid head groups is also a factor that defines the compressibility of the membrane with respect to the location of protein molecule [37]. The main purpose of the Z value analysis is to understand the effect of protein insertion on the lipid membrane accumulation. That is, a compact accumulation of protein is represented by a low Z value of the lipid in the XY-plane and occurs in most of the simulation systems, as shown by the dark regions in Figs. 7b and c.
As shown in Figs. 7b and c, the low Z value occurs for protein location in the middle of the membrane. Fig. 7(c, left) for system 1 indicates the highest thickness and lowest Z value, which can be interpreted by less stability of protein and less compressibility of the membrane.

Order parameter
Movements of lipids, such as rotations around the lipid axis and chemical bonds, fluctuations and other kind of motions in a fluid bilayer, are taking place in a very short time-scale ranging from picoseconds up to milliseconds. The order parameter can be used to analyse the membrane property and compare the simulation results. The order parameter of the lipid chain is defined as [41] S = 3 cos 2 u − 1 2 where θ is the angle between the molecular vector and the vector that is parallel to the bilayer normal (Z ). The angle brackets in the above equation refer to an ensemble time average. For a completely aligned bond, S bond = 1; for a  completely random (isotropic) bond, S bond = 0; and for a bond perfectly perpendicular to the bilayer normal, S bond =½.
Moreover, the membrane entering into a gel phase during the simulation can be verified by order parameter analysis [42].
In Fig. 8, the order parameters for CG and UA simulation models for two chains of DPPC interacting with two proteins are represented. As seen in this figure, the system with UA-DPPC and protein 2 has a higher order compared with the systems with protein 1. As mentioned earlier, the systems with protein 2 are more stable, which means that the chains are more ordered in the membrane. There are differences between the order parameters in UA-models and these differences are more obvious in DPPC chain 2. Although the CG force field lead to a slight  underestimation of the tail order, it represents important structural features corresponding to the lipid bond order and that the lipids have more order compared with the UA-models. In addition, good overall agreement for all lipid bonds is observed between CG and UA simulation results and the (systems 3 and 4) have almost the same order parameters. In Fig. 8a, it is seen that the lipids in system 2 have the lowest order which can be due to more fluctuations of protein in system 2. This appears as a larger dark region in Fig. 7(b, right) which corresponds to the location of protein in the membrane.

Conclusion and final remarks
Large-scale MD computer simulations were performed on the ion pump protein Ca 2+ -ATPase, in a lipid bilayer DPPC membrane. Simulations with 10 ns duration show that the CG-models have almost the same behaviour in this time scale the protein-bilayer in CG models is more stable than AA/UA. Transient of the studied systems including protein 1 and lipid bilayer structures from their initial configurations to the final relaxed equilibrium configurations takes about 4 ns for CG model (system 3) and for AA/UA model (system 1) takes about 6 ns. The timescales for protein 2 are 2 ns for CG and 6 ns for AA/UA models. The timescale will of course depend on how far the initial configuration is from the equilibrium state. However in AA/UA systems, long-term simulation is needed to attain the equilibrium state. The results indicate that the protein 2 with the initial structure obtained by X-ray diffraction is more stable than protein 1 with the initial structure obtained by EM. This difference can be explained by fact that in constructing the crystallographic structure of protein based on X-ray diffraction pattern, the detailed structure of protein which is highly influenced by the presence of calcium ions has been taken into account, whereas in EM pattern this detailed structure cannot be considered. Furthermore, the simulations provide an approach to quantify the fluctuations of protein, and the results indicate a good agreement between application of CG-models and AA-models. The results also show more dynamics in EM-structural protein compared with the X-ray-structural protein and this can explain why the EM-structural protein has a higher diffusion coefficient in the bilayer membrane and consequently the membrane is less stable. This is manifested by smaller area per lipid and larger thickness in comparison to other systems. The activity of Ca 2+ -ATPase because of its structural change affects the membrane thickness and therefore the higher activity of protein can bring about larger thickness for the membrane. By assessing the order parameter of the lipids, it is found that CG-models are more ordered than UA-models. In addition, the systems with protein 2 are more ordered than the systems with protein 1.
By comparing systems 1 and 2 (AA-UA systems), in Table 4, it can be seen that the membrane in system 1 has larger thickness and therefore should have smaller area per lipid, and as a result the lipid chain will be more ordered. This trend appears to follow in CG systems (systems 3 and 4) as well: the thicker membranes have smaller area per lipid and are more ordered. That is, when the lipid is located perpendicular to the membrane surface, it occupies less area and is more ordered.
In system 3, a conformational change occurs in the protein structure, which appears as a shift in the RMSD value (according to Fig. 3) and it causes more fluctuations in the protein structure. As has been reported previously [43], the activity of protein is related to the thickness of the bilayer membrane, where the protein has been embedded, and how the protein and membrane are matched with each other in terms of their hydrophobicity. Therefore, these bring about an optimum thickness for the membrane where the highest activity of the protein will occur. Hence, thick or thin membranes may decrease the activity of protein or make the conformational change in the protein structure. The results obtained indicate that protein in system 3 undergoes a conformational change which is manifested by higher fluctuations compared with other studied systems.