A computational model to unravel the function of amyloid-β peptide in contact with a phospholipid membrane

The normal synapse activity involves the release of copper and other divalent cations in the synaptic region. These ions have a strong impact on the membrane properties, especially when the membrane has charged groups, like it is the case of synapse. In this work we use an atomistic computational model of dimyristoyl-phosphatidylcholine (DMPC) membrane bilayer. We perturb this model with a simple model of divalent cation (Mg2+), and with a single amyloid-β (Aβ) peptide of 42 residues, both with and without a single Cu2+ ion bound to the N-terminus. In agreement with experimental results reported in the literature, the model confirms that divalent cations locally destabilize the DMPC membrane bilayer, and, for the first time, that the monomeric form of Aβ helps in avoiding the interactions between divalent cations and DMPC, preventing significant effects on the DMPC bilayer properties. These results are discussed in the frame of a protective role of diluted Aβ peptide floating in the synaptic region. Author summary We modelled the behavior of a Mg2+ divalent cation, with the size of Zn2+ and Cu2+, in contact with a phosphatidyl lipid bilayer. We also modelled the monomeric amyloid-β peptide 1-42, both free and Cu-loaded, the latter mimicking the final step of the binding between the peptide and the divalent cation. On the basis of the simulation results, we propose that the peptide hinders the strong interactions between the divalent cation and the membrane.

Set-up of molecular dynamics simulations 79 The amyloid-β peptide of 42 residues (Aβ(1-42)), with and without a single bound 80 copper ion in 2+ oxidation state (Cu 2+ ), was simulated with constant temperature 81 molecular dynamics (conventional MD, CMD hereafter) and with the replica exchange 82 molecular dynamics (REMD) method, in order to sample the configurational statistics 83 at the physiologically relevant temperature of 311 K (38 • C). The peptide and the ions 84 were put in contact with a bilayer composed of 85 1,2-dimyristoyl-sn-glycero-3-phosphocholine (also abbreviated as 86 dimyristoyl-phosphatidylcholine, DMPC hereafter) lipid molecules. 87 The sequence of Aβ(1-42) is: 88 DAEFRHDSGY 10 EVHHQKLVFF 20 AEDVGSDKGA 30 IIGLMVGGVV 40 IA 89 with aminoacids indicated with the one-letter code. We used the Amber16 package [56], 90 with the FF14SB [57] force-field for the peptide and monovalent ions (KCl), TIP3P 91 water model [58] for the explicit water solvent, and LIPID14 [59] for the DMPC 92 molecules. AMBER FF14SB force-field is an improved version of FF99SB [60] used in 93 our previous simulations [52,61]. Although older CHARMM force-fields tend to provide 94 better results for Aβ peptide than old AMBER force-fields [62], new AMBER versions, 95 especially AMBER FF14SB and CHARMM36m provide good agreement with 96 experimental data for Aβ [63,64]. Moreover, AMBER FF14SB is fully compatible with 97 LIPID14 force-field [59], which is expected to provide optimal accuracy for both lipids 98 and peptide in the simulations. The use of more recent force-fields for IDPs will be 99 pursued in the future, after a detailed comparison between experiments and simulations 100 in generalized ensembles will be reported in the context of amyloid peptides. 101 We assumed the physiological (pH∼7) protonation state for aminoacid sidechains 102 and free termini. Thus, the charge of Aβ(1-42) is -3 (the N-terminus is protonated and 103 the C-terminus deprotonated). The parameters for copper and copper-bound 104 aminoacids were the same used in our previous MD and REMD simulations [52,61]. Cu 105 is bound to N and O of Asp 1, Nδ of His 6 and N of His 13, the latter protonated at 106 Nδ. His 14 is neutral and protonated in N , like His 6.

107
Bond distances and angles involving Cu contribute to harmonic energy terms, with 108 stretching constants, bending constants, and equilibrium values set as fitting parameters 109 of quantum-mechanics calculations at the density-functional level of approximation for 110 truncated models (see Methods in [52]). All the dihedral angles where Cu has index 2 or 111 3, do not contribute to the potential energy, while those with Cu with index 1 or 4 are 112 obtained by the AMBER99SB force-field where heavy atoms have the same dihedral 113 indices of Cu. Point charges are derived from the restrained electrostatic potential 114 (RESP) method [65,66], where the electrostatic potential mapped onto the  Lennard-Jones parameters for Cu are reported in the literature [67]. The Cu 2+ 120 coordination geometry in this empirical force-field is approximately square-planar, with 121 a fifth axial coordination always available to electrostatic interactions, as shown in 122 previous simulations performed with the same force-field [36]. The distance between 123 configurations obtained with this empirical force-field and minimal-energy 124 configurations obtained including explicit electrons (like in density-functional theory 125 applied to truncated models) is small. 126 As for the free divalent cation, we used the so called "dummy" cation model for 127 Mg 2+ [47]. This model has been used together with AMBER99SB phosphate 128 groups [48], where it showed reasonable electrostatic properties. Even though this 129 model is a very crude approximation of divalent cations, it is far more reliable than a 130 single site with point charge 2+. A comparison between the affinity of divalent and 131 monovalent cations for the DMPC membrane has been performed by umbrella sampling 132 estimates of free energy differences (see Supporting Information, file ).

133
An initial lattice model of DMPC bilayer was built, using 77 DMPC molecules per 134 layer, with an approximate area per molecule of 62Å 2 . An orthorhombic simulation cell 135 was built, with the cell side along zeta, the latter the direction normal to the DMPC 136 layer, initially set to 70Å. The space between the periodic images of the bilayer was 137 filled with 13511 water molecules, initially at the density of 1 g/cm 3 , according to the 138 TIP3P model of bulk water at room conditions. KCl was added in the same space, 139 according to an approximate bulk concentration of 0.1 M. Ions were added randomly 140 replacing water molecules in the initial configuration. The number of Clanions was 141 adapted to the change of net charge due to addition of the peptide (see below). The net 142 charge of the simulation cell was always zero.

143
Initial configurations of amyloid−β monomer, without copper (charge 3-) and with 144 copper (charge 2-, because of N-terminus deprotonation), were inserted in the space 145 filled by the water molecules. The same was done for the single divalent cation. The  To remove eventual bad contacts produced by each initial configuration set-up, we 152 performed 25000 steps of steepest decent energy minimization, followed by other 25000 153 steps of conjugate gradient energy minimization.

154
The initial coordinates for the CMD and REMD simulations are included as 155 Supporting Information in the protein data-bank (PDB) file format (the first 156 configuration) and as compressed (Bzip2) XYZ format. See , , , , , and .

157
Molecular dynamics simulation protocol 158 We simulated MD trajectories in the isobaric-isothermal (NPT) statistical ensemble, at 159 the constant temperature of 311 K and at the pressure of 1 Atm. Temperature was 160 controlled by a Langevin thermostat [68] with collision frequency of 2 ps −1 . Pressure 161 was controlled by a stochastic barostat, with relaxation time of 100 fs. The SHAKE 162 algorithm [69] was applied to constrain bonds involving hydrogen atoms. A cut-off of 163 10Å was applied for non-bonded interactions and the particle mesh Ewald 164 algorithm [70] was used to compute long-range Coulomb and van der Waals interactions. 165 The simulation time-step was 2 fs. 166 In order to increase the sampling, we collected several trajectories for each system, 167 starting from different initial conditions, that is the initial velocities (DMPC) and the 168 position of the solute within the water layer (other systems). Composition of each 169 system and some parameters related to sampling are reported in Table 1. The REMD simulation was carried out with 56 replicas (or trajectories) corresponding 172 to 56 temperatures ranging from 273K to 500K. The minimized structure was 173 distributed among 56 replicas, and each replica was equilibrated in 200000 steps at the 174 temperature chosen in the temperature distribution. After equilibration, the replica 175 exchange molecular dynamics simulation started, lasting a total time of 400 ns. The 176 exchange of temperature between pair of replicas was attempted every 500 steps of 177 simulation. The REMD simulation is used here mainly to capture the statistical 178 contribution of extended peptide configurations and partially disordered layers, 179 configurations that are rarely sampled at temperatures in the range where the 180 force-field is accurate. The acceptance rate of REMD simulations was, on average, 20% 181 and 21% for, respectively, Aβ(1-42) and Cu-Aβ(1-42).

182
The behavior of lipid order parameters as a function of temperature (data not shown 183 here) shows that the DMPC bilayer is, at the temperature closest to that of the human 184 body (37 • C, 310 K), in the liquid crystalline phase. The configurations sampling the 185 temperature of 311 K are, therefore, analyzed in detail in the following.

186
To avoid possible bias due to initial configuration construction, we used 187 approximately the second half of each simulation for analysis (see Table 1).

199
Radial distribution function (RDF) measures the probability to have the distance 200 between two sites within a given distance range, N (r). As usual for liquids and 201 polymers, this quantity is then divided for the same probability for the ideal gas with 202 the same uniform density of sites, N id (r): g(r) = N (r)/N id (r). The function g(r) 203 approaches the limit g(r) =1 when r → ∞, i.e. when the two sites in the pair become 204 not correlated.

205
The bilayer thickness is defined as the distance between the two planes formed by 206 phosphor atoms belonging to each layer. The roughness of a layer is defined as the 207 standard deviation of z coordinates of phosphor atoms within each layer.

208
The number of contacts is defined as the count of the usual distance-based step-like 209 variable: (1) with i and j running over different sets of atom pairs, each term of the pair contained in 211 a different portion of the system. When the two sets of atoms identify, respectively, 212 atoms belonging to positively charged groups (Nζ in Lys and Nη in Arg) and negatively 213 charged groups (Cγ in Asp and Cδ in Glu), we address the contact as an intramolecular 214 salt-bridge. The number of such contacts is indicated as SB and the d 0 parameter is 215 chosen as 4Å. As for generic inter-residue contacts, we measured the distance between 216 the centers of mass of sidechains in the two involved residues. In this case, d 0 is chosen 217 as 6.5Å. When the contact between aminoacids and lipid molecules is addressed, the 218 center of mass of DMPC molecules is used and the d 0 distance is 4.5Å.

219
Elastic moduli 220 Elastic moduli of lipid bilayer were calculated by fitting suitable ensemble averages with 221 the following equations [72]: (2) January 14, 2020 6/29 where K c , K Θ , K tw are bending, tilt, and twist elastic moduli, respectively, k B is 223 Boltzmann constant, T is temperature, andn q is the reciprocal space vector determined 224 as summarized below (see also supplementary information of Ref. 72 and Ref. 73).

225
The xy plane of the membrane is discretized to a square 8×8 grid. The orientation 226 vector of lipid molecule j is n (α) j (x, y, z) with α 1 or 2 for upper and lower layers, 227 respectively. Each vector points from the midpoint between P and C2(glycerol) atoms 228 to the midpoint between the terminal C atoms of the lipid tails. The orientation vectors 229 are projected onto the xy plane and are mapped onto the 8×8 grid, providing n (α) (x, y). 230 Fast Fourier transform is used to obtain n (α) q , where q is the reciprocal space index.

231
From n (α) q we obtain the quantity that is decomposed into longitudinal (n q ) and transverse (n ⊥ q ) components: Finally Eq. 2 is used to average according to the collected sampling of lipid molecules. 234

235
Addition of a divalent cation to the DMPC bilayer 236 The affinity of Mg 2+ for the DMPC bilayer was measured by the umbrella sampling 237 method (see Supporting Information, Fig.S1 in ). The free energy minimum was found 238 at 17Å from the bilayer center, thus corresponding to the average distance of P atoms 239 (see below). The flatter shape of free energy around the minimum in the case of Na + is 240 due to the equivalent interactions of Na with phosphate and carbonyl groups of DMPC. 241 These interactions allow a deeper penetration of Na into the bilayer than Mg. should favour Na compared to Mg, being the hydration free energy at 300 K about five 246 times more negative for Mg compared to Na [74]. This effect is due to the strong  Methods), the average over the 3 trajectories is analyzed in the following. We indicate 258 the cation-bound layer as layer 1 (L1) and the layer not affected by the binding as layer 259 2 (L2). The difference between g calculated for L1 and L2 is displayed in Fig 1. The 260 divalent cation (black) is bound to the phosphate oxygen atoms, thus displaying the 261 coordination distance of 2.9Å with respect to P atoms. Including the second-shell P 262 atoms (the peak at 3.5Å), the number of P atoms around the cation is 4. This 263 coordination affects the average distance between charged groups within L1, as it is 264 displayed by the P-P distances (red line), respectively, within each layer L1 and L2. 265 Conversely, atoms farther than P from the perturbing cation are less affected, as shown 266 by the difference in N-P distance distribution among the two layers (blue line).

275
The low symmetry of K-P distribution in the absence of Mg (red curves) is due to 276 sampling limitations. Indeed, the presence of Mg on the L1 layer displays a "hole" in K 277 distribution where there is a little excess in the absence of Mg. Because of the change in 278 interactions between K + and P at short distance (the peaks at the left), there is also a 279 decrease of bulk concentration within a distance of 1 nm from the P atoms. This change 280 of the electrostatic properties between the two sides of the bilayer is equivalent to a weak 281 polarization of the membrane. This asymmetry is caused by the asymmetry in the P-P 282 radial distribution (Fig 2B), that is due to the formation of the Mg-O(P) coordination. 283 The asymmetry of the interactions between divalent cations added from one side of 284 the bilayer is consistent with experimental data reported for exogenous addition of Cu 2+ 285 and Zn 2+ to bilayer models (POPC/POPS mixtures) [53]. The comparison between 2 H 286 and 31 P ss-NMR spectra of POPC/POPS molecules shows that P atoms are strongly  The effect of the divalent cation on the elastic property of DMPC is also significant. 291 In Table 2 we report the elastic constants determined by the different simulations, with 292 averages of Eq. 2 (see Methods) computed over all the acquired trajectories (see 293   Table 1).

294
The values are in the range of those found in DPPC atomistic simulations [72],    The effect of Mg addition to L1 does not significantly alter other structural 313 parameters of the bilayer at the same temperature (see Table 3 [75]); average of 10 CMD simulations for Aβ/DMPC at 303 K and DMPC at 311 K (blue squares).
In Fig 5 we display the radial distribution function g for selected pairs, to show the 338 extent of penetration of N-and C-termini (respectively Nt and Ct) through the 339 membrane surface (using P atoms in the pair) or towards the membrane center (using 340 the terminal C atom in the two acyl chains of DMPC, Cf hereafter). The g function is 341 measured at T =311 K, i.e. the physiological temperature of the synaptic membrane.

342
The REMD trajectory at 311 K shows that the propensity for Aβ and Cu-Aβ N-termini 343 to interact with the membrane surface is limited to the head groups of the DMPC   The representation of this change in penetration is better understood examining the few 375 snapshots contributing to the contribution to g at short distances in, respectively, Cf-Nt 376 (Aβ/DMPC, Fig 5C) and Cf-Ct (Cu-Aβ/DMPC, Fig 5D). In Fig 6 we display, left and 377 right panels, one of such configurations for, respectively, each of the two systems. It can 378 be observed that a common feature of the peptide structure in these configurations is 379 the breaking of cross-talk between the N-and C-termini. This cross-talk is always  The number of intramolecular salt-bridges (SB) within the peptide (Table 4)   The bilayer structure (Table 3) shows only a moderate propensity for larger thermal 398 fluctuations, induced by the perturbation due to weak interactions with the peptide, 399 and a small increase in thickness.  Table 3 and 427 discussion below), thus the average distance between P atoms and the central plane of 428 the bilayers is never below 17Å. The approach of Mg towards the bilayer central plane 429 does not significantly drift, on average, the P atoms towards the center of the bilayer, 430 because the density of P atoms projected along the z axis does not change (data not 431 shown here). However, as described above, the perturbation makes a little hollow over 432 the bilayer surface affected by Mg binding (see Fig 3 and   The area per lipid as a function of temperature measured by REMD simulation (see 441 above) and consistent with experimental data [75] shows that the area per lipid 442 increases with temperature. Therefore, most of the changes displayed in Table 3   Cu-Aβ/DMPC, we notice that the more symmetric is the interaction between the 472 peptide among the two layers (left panels), the more symmetric is the distribution of K + 473 (right panels). It is also interesting to notice that the strong interaction of trajectory 1 474 for Aβ/DMPC (see above), produces a polarization of K + that is opposite to that 475 produced by Mg 2+ (Fig 2A, black curve).

512
The number of intramolecular salt-bridges is, on average over the 10 trajectories, not 513 altered in the presence of DMPC with respect to the case of water solution (Table 4).

514
The SB quantity increases when the association of the peptide with DMPC is more  Trajectory 5 (middle) ends with configurations significantly penetrating the membrane 533 bilayer, but with interactions almost confined to the surface. Finally, in trajectory 1 the 534 peptide rapidly achieves the penetration of the bilayer from the side of its C-terminus 535 (bottom). In the latter conditions, it can be noticed that the region of Aβ crossing the 536 layer surface is small, separating the N-terminus (above the surface) and the C-terminus 537 (below the surface). This configuration, again, represents the requirement of removing 538 the cross-talk between the N-terminus and the C-terminus (exerted by the bending of 539 N-terminus towards the C-terminus) before a deeper penetration of the peptide into the 540 membrane from the side of the C-terminus. This configuration is similar to that 541 obtained by REMD of Cu-Aβ/DMPC displaying the deepest penetration into the 542 bilayer ( Fig 6B), with the main difference that the N-terminus is not partially 543 neutralized by Cu binding.  freezing, displayed by an increase, with respect to the free peptide, of highly populated 561 contacts between residues far in the sequence. Some of them involve Glu 22,Asp 23,562 and Lys 28, with these charged sidechains interacting mostly with the N-terminus and 563 not between themselves. In the case of a peptide that is more significantly embedded 564 into the bilayer (trajectory 1, bottom panel), one notices the disappearing of contacts 565 within residues in the C-terminus and the extension of the N-terminal domain up to Lys 566 28, with the void observed for trajectory 2 (top) almost filled. This change in cross-talk 567 is induced by the formation of contacts between the C-terminus and DMPC. In the 568 right panels of the same figure, we display the mass density for different atomic sets in 569 Aβ(1-42). S1 is the N-terminus, S4 the C-terminus, while S2 is the hydrophobic segment 570 and S3 contains the charged residues involved in one of the intramolecular SBs. When 571 the peptide is out from the bilayer (trajectory 2, top-right panel) only the N-terminus 572 (S1) is approaching the bilayer surface. The analysis of the trajectories not displaying 573 penetration into the bilayer (all trajectories except 1 and 5, data not shown here) shows 574 that there is no preference among the different segments for weak interactions with the 575 bilayer surface. When a more significant interaction with the bilayer surface occurs 576 (trajectory 5, middle-right panel) the hydrophobic segment S2 is projected towards the 577 bilayer because of the stronger interactions among S3 and S1 (as shown in the 578 middle-left panel). When the penetration is deeper (trajectory 1, bottom-right panel), 579 the S4 segment overtakes the layer of P atoms, with the latter interacting with S3. 580 Interestingly, in these conditions the S2 segment is projected towards the water layer,  and diffraction studies [81]. Consistently, dramatic changes of peptide/membrane 619 interactions are observed at conditions where the peptide is truncated to be more 620 hydrophobic (Aβ(25-35)) or forms fibril assemblies [81,82].

621
The picture of Aβ monomers floating over the membrane surface is consistent with 622 other observations reported in the literature. A recent FRET experimental work [4] 623 describes the strong interactions among growing fibrils and the DOPC membrane, 624 modelled as a lipid vesicle. The same study confirms that monomers do not directly 625 bind the lipid bilayer, as already observed in previous studies.

626
As for the impact on oligomer formation, our results point out the possible role of 627 charged groups of the bilayer in organizing monomers into oligomers. Indeed, several 628 simulations showed that a strong association between Aβ and zwitterionic and charged 629 membranes occurs starting from tetrameric Aβ assemblies [15]. Since it is known that 630 the lag-time of monomers associated to Cu is larger than that of Cu-free Aβ when in 631 the water solvent [80], it is not surprising that the DMPC association with Aβ(1-42) in 632 the absence of divalent cations does not decrease the chance of inter-monomer contacts 633 compared to the water solution. The bilayer-water interface, when bilayer has charged 634 groups on the surface, exerts a mild attraction for Aβ(1-42) thus decreasing the freedom 635 of monomers by reducing the space dimensionality. Conversely, at oligomeric level the 636 bilayer surface can assist the formation of larger oligomers and protofibrils. This type of 637 association has been observed in models of preformed protofibrils interacting with lipid 638 bilayers [14]. 639 We notice here that in ss-NMR experiments, the effect of the addition of free Cu 2+ 640 and Zn 2+ ions on the membrane properties is more dramatic than in the presence of Aβ 641 peptide [53]. Similar strong effects have been observed both experimentally and 642 computationally for free Ca 2+ ions [41][42][43], and Mg and Cu divalent cations are even 643 smaller than Ca 2+ in size. For the first time, we show in this study that the Aβ-bound 644 Cu 2+ ion does not exert the strong perturbation on the membrane exerted by a free 645 divalent ion. Indeed, the effect of Cu-Aβ monomer on the membrane is weaker than 646 that of the more charged Aβ peptide.

647
Therefore, the formation of the Cu-Aβ complex before an eventual incorporation into 648 the membrane and before an increase in peptide concentration, appears as a protection 649 against membrane destabilization and oxidation. This hypothesis is confirmed by the Cu [54,81]. Since the N-truncated peptide does not bind Cu, the addition of Cu to the 652 system has an effect on the bilayer that is similar to that of free Cu.

654
We perturbed an atomistic model of the DMPC bilayer, representing a very crude 655 approximation of a portion of the cellular membrane facing the synaptic region between 656 neurons, with a single divalent cation (Mg 2+ ) and with Cu-free and Cu-loaded 657 amyloid-β peptides of 42 aminoacid residues in the monomeric form.  The model has many limitations. Beyond the limitations in size and number of 666 components, that are common to application of atomistic models, there is the lack of 667 working approximations to interactions between a ion like Cu 2+ , with available 3d 668 orbitals, and molecules providing a pletora of possible ligand atoms, like phosphate, 669 carboxylate, imidazole, and carbonyl groups, not to mention deprotonated amide 670 backbone nitrogen that are known to bind Cu 2+ at physiological pH. These limitations 671 will be removed by polarizable and reactive force-fields, not yet available.