Molecular Dynamics Simulations of the Permeation of Bisphenol A and Pore Formation in a Lipid Membrane

Bisphenol A (BPA) is particularly considered as one of the most suspicious endocrine disruptors. Exposure to BPA may bring about possible human toxicities, such as cancerous tumors, birth defects and neoteny. One of the key issues to understand its toxicities is how BPA enters cells. In this paper, we perform molecular dynamics simulations to explore the interactions between BPA and a phospholipid membrane (dipalmitoylphosphatidylcholine, DPPC bilayer). The simulation results show that BPA can easily enter the membrane from the aqueous phase. With the increasing concentrations of BPA in the membrane, BPA tends to aggregate and form into cluster. Meanwhile, several DPPC lipids are pulled out from each leaflet and adsorbed on the cluster surface, leading to pore formation. Detailed observations indicate that the lipid extraction results mainly from the dispersion interactions between BPA cluster and lipid tails, as well as weak electrostatic attractions between lipid headgroups and the two hydroxyl groups on BPA. The lipid extraction and pore formation may cause cell membrane damage and are of great importance to uncover BPA’s cytotoxicity.

All MD simulations were carried out under the isothermal-isobaric (NPT) ensemble using the Gromacs package 4.5.5 26,27 and GROMOS53a6 force field 28 . Periodic boundary conditions were employed in all directions. The vdW interactions were treated with smooth cutoff at a distance of 12 Å, whereas the particle-mesh Ewald method was adopted to calculate the long-range electrostatic interactions 29,30 . Water was represented by the SPC model 31 . The temperature was kept stable at 323 K using the V-rescale thermostat 32 and the pressure was controlled semi-isotropically by a Berendsen barostat 33 . Bond lengths within BPA/DPPC and water molecules were constrained by the LINCS and the SETTLE algorithms 34,35 , respectively.
The free energy of BPA across the bilayer was computed from the potential of mean force (PMF) using umbrella sampling 36 . First, we conducted steered MD simulation to pull the molecule from the aqueous phase to the middle of the membrane. Then, 30 configurations were generated along the z-axis direction (reaction coordinate). The z coordinates of COM distance between BPA and membrane in each configuration differed by about 0.1 nm. Each window was equilibrated for 5 ns and a production run of 5 ns was continued for sampling. Eventually, the PMF profile was obtained by the Weighted Histogram Analysis Method (WHAM) 37 , implemented in the GROMACS package as 'g_wham' 38 . The convergence and sufficiency and of each sampling were presented in the Supplementary Information Figures S1 and S3.

Results and Discussion
The permeation of BPA in the membrane. Because of the extremely low water solubility, we only put one BPA in the aqueous phase and one in the bilayer center (see Fig. 1A) to look for the preferable location of the molecules. Figure 1 showed the initial and final structures and the time evolution of the z-coordinates of their COMs. During the first 25 ns, the outer BPA just moved randomly in the water. However, at t~25.32 ns, the outer BPA rapidly entered the bilayer and then stayed at the interface between water and the lipid headgroups. The detailed dynamic process of the outer BPA entering the lipid bilayer was presented in Fig. 1D-G. First, BPA was adsorbed on the surface of bilayer because of the electrostatic interactions between the two hydroxyl groups of BPA and the lipid headgroups (see Fig. 1D, t = 23.76 ns). Due to the hydrophobic interactions between BPA and lipid tails, the main body of BPA rotated and then was pulled into the lipid tails (see Fig. 1E-G), while the hydroxyl groups were still near the headgroups. Such location and orientation had no essential changes but slight fluctuations for the rest of the simulation. Meanwhile, the inner BPA went straight and quickly to another interface and remained there with the similar location and orientation till the end of simulation. The two symmetric locations were in line with NMR determination 39 . We have rerun the system with more accurate partial charges, obtained by DFT calculation at wb97xd/6 − 31 + g(d, p) level. As shown in Figure S4 in the Supplementary Information, the results are in good agreements with those in Figure 1.
The translocation of BPA can be more quantitatively illustrated by the PMF of BPA at different transmembrane positions, as shown in Fig. 2B. The PMF of BPA in the aqueous phase far away from the interface was defined as zero. We observed that BPA penetrating into the bilayer should overcome an energy barrier, though it was not high (~0.6 kcal/mol). However, the free energy difference between BPA at the center and at the interface reached 5.1 kcal/mol (Δ G pen ), resulting in its fast moving towards the lipid headgroups. The minimum of the PMF was − 4.5 kcal/mol at the symmetric positions ± 1.1 nm from the bilayer center. Mass density profile (see Fig. 2A) showed clearly that the preferable locations of BPA were symmetric and close to the lipid headgroups (the densest region of the membrane). These two locations were in good agreement with the two minima of the PMF, implying that they were energetically favorable. Figure 3 showed the detailed interactions among DPPC, BPA and water. The orientation of BPA in the membrane was that its hydroxyl groups pointed to the hydrophilic interface and the rest of the molecule was immersed in the lipid tails (see Fig. 3C), due to the following two facts: 1) the electrostatic interactions between the two hydroxyl groups on BPA and the lipid headgroups; 2) the hydrophobic interactions between the two benzene rings and methyl groups on BPA and the lipid tails. Radial distribution functions (RDF) were calculated to elucidate the former. The first peak between polar hydrogen atoms on BPA and oxygen atoms on phosphate groups (Fig. 3A, black line, at ~0.18 nm) corresponded to a hydrogen bond, while the first peak between oxygen atoms on BPA and phosphorus atoms on phosphate groups (Fig. 3A, red line) was at ~0.39 nm, close to their collision radii. These peaks indicated that the hydroxyl groups on BPA strongly interacted with the polar lipid headgroups.
Water molecules could easily penetrate into the interspaces between lipid headgroups, and act as bridges between lipid headgroups as well as between BPA and headgroups. Figure 3B showed the RDF of water molecules around BPA. The first (at ~0.18 nm) and second peak (at ~0.35 nm) represented the hydrogen bonds, meaning that water in the membrane preferred to surround BPA. After this, the RDF rose gradually with the increasing distance, since more and more water molecules were included out of the membrane.
The permeation of BPA cluster in the membrane. We then increased the concentration of BPA in the water to 0.43 mol%, in which mol% was defined as the number of BPA divided by the sum of BPA and water molecules. As shown in Fig. 4A, BPA formed into cluster quickly at t = 12 ns, and then moved randomly in the water. At about t = 150 ns, BPA cluster began to slowly enter the membrane. Figure 4B showed that the difference between the COMs of BPA cluster and lipid bilayer was close to zero after entering the membrane, since the COM of BPA cluster was near the bilayer center (see Fig. 4A, snapshot at t = 300 ns). The BPA cluster was kept almost intact except that 3 BPA molecules went away from the cluster. Mass density profile (Fig. 4C) also showed that different from single BPA in the membrane, BPA cluster preferred to stay in the middle of the bilayer. That is to say, BPA presents low lipid solubility.
To confirm the low solubility of BPA in the lipid membrane, we randomly embedded BPA molecules in the bilayer center with different concentrations, namely 3.8 mol%, 7.2 mol%, 10.5 mol%, 13.5 mol%, 16.3 mol% and 23.3 mol%, Correspondingly, the number (abbreviated to 'n') of BPA was n = 5, 10, 15, 20, 25 and 45. Here, mol% was defined as the number of BPA divided by the sum of BPA and DPPC molecules. The results were given in Fig. 5. At extremely low concentration, BPA could disperse in the membrane (see Fig. 5A) and BPA was close to the lipid headgroups. Then, BPA began to aggregate into small clusters, when the concentration was about 7.2~13.5 mol% (see Fig. 5B-D). And BPA formed into big clusters, when the concentration was more than 16.3 mol% (Fig. 5E,F). In general, BPA has weak lipid solubility according to our simulation results.
Interestingly, it was found that pore formed at the entrance (see Fig. 6A) where BPA cluster entered the membrane. And water molecules flew in and filled the hole. Detailed observation showed that these holes were located on both sides of BPA cluster, and several lipids broke away from the bilayer and were adsorbed on the surface of the cluster. The pore formation was also found in the systems, in which the concentration of BPA in the membrane was more than 16.3 mol%. It has been reported that the dispersion interactions between sp 2 carbon atoms and lipid tails are stronger than the self-attraction among the lipids in the bilayer 40 . Clearly, this strong attraction between BPA cluster and lipids originates from such dispersion interactions between the benzene rings on BPA and lipid tails. Figure 6C showed the number of water molecules in the bilayer as a function of simulation time. These water molecules were mainly distributed on both sides of the BPA clusters (see Fig. 6B), where the holes existed.

Conclusions
The translocation and aggregation of BPA in a DPPC bilayer have been investigated by molecular dynamics simulations. Under low concentration, BPA can readily penetrate into the membrane and prefer to stay at the interface between water and lipids, with its hydroxyl groups towards the lipid headgroups and its carbon rings in the lipid tails. This location and orientation enable electrostatic interactions between the hydroxyl groups of BPA and the polar headgroups of DPPC as well as water, which can be confirmed by hydrogen bonds found in the RDFs. With high concentration, BPA clusters first and then enters the membrane, leaving a hole at the entrance. The aggregation of BPA in the membrane can also cause pore formation, since several lipids are extracted from each leaflet and adsorbed on the cluster surface. The dispersion interactions between lipid tails and the benzene rings on BPA as well as electrostatic attractions among their polar groups contribute to the lipid extraction. The results provide detailed atomic mechanism of the interactions between BPA and a model cell membrane. And BPA may possess potential cytotoxicity when clustered in the membrane.