Characterization of Interactions between Curcumin and Different

Feb 2, 2018 - In this study, the interactions of CUR with different types of model lipid bilayers (POPE, POPG, POPC, DOPC, and DPPE), mixtures of mode...
1 downloads 4 Views 3MB Size
Subscriber access provided by UNIVERSITY OF TOLEDO LIBRARIES

Article

Characterization of Interaction Between Curcumin and Different Types of Lipid Bilayers by Molecular Dynamics Simulation Yuan Lyu, Ning Xiang, Jagannath Mondal, Xiao Zhu, and Ganesan Narsimhan J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b10566 • Publication Date (Web): 02 Feb 2018 Downloaded from http://pubs.acs.org on February 6, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Characterization of Interaction between Curcumin and Different Types of Lipid Bilayers by Molecular Dynamics Simulation Yuan Lyu1, Ning Xiang1, Jagannath Mondal2, Xiao Zhu3*, and Ganesan Narsimhan1* 1

Department of Agricultural and Biological Engineering, Purdue University, West Lafayette, IN 47907 USA

2

Tata Institute of Fundamental Research, Centre for Interdisciplinary Sciences, 36/P, Gopanapally Village, Serilingampally Mandal, Ranga Reddy District, Hyderabad-500107, India 3

Research Computing, Rosen Center for Advanced Computing, Purdue University, West Lafayette, IN 47907 USA

* Authors to whom correspondence should be addressed: Ganesan Narsimhan. E-mail: [email protected] Xiao Zhu. E-mail: [email protected]

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Curcumin (CUR) is a natural food ingredient with known ability of targeting microbial cell membrane. In this study, the interaction of CUR with different types of model lipid bilayer (POPE, POPG, POPC, DOPC, and DPPE), mixture of model lipid bilayer (POPE/POPG), and biological membrane mimics (E. coli and yeast) were investigated by all atom explicit solvent molecular dynamics (MD) simulation. CUR readily inserts into different types of model lipid bilayer systems in liquid crystalline state, staying in the lipid tails region near the interface of lipid head and lipid tail. Parallel orientation to the membrane surface is found to be more probable than perpendicular for CUR as indicated by the tilt angle distribution. This orientation preference is less significant as the fraction of POPE is increased in the system, likely due to the better water solvation of perpendicular orientation in POPE bilayer. In E. coli and yeast bilayers, tilt angle distributions were similar to that for POPE/POPG mixed bilayer with water hydration number around CUR for the former being higher. Insertion of CUR resulted in membrane thinning. The results from these simulations can provide insights into the possible differences in membrane disrupting activity of CUR against different types of microorganisms.

2 ACS Paragon Plus Environment

Page 2 of 34

Page 3 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1. Introduction Natural antimicrobial molecules, such as enzymes

1-3

, bacteriocins

4-6

, essential oil

7-9

, organic acid

10-11

,

and phenolic compounds 12-18, have attracted great attention in controlling infectious food borne outbreaks in recent decades. Among these agents, phenolic compounds showed capability of interaction with membrane by altering properties of the cell

19-21

, while the exact mechanism responsible for this altering

abilities remain largely unclear. Therefore, curcumin (CUR [1,7-Bis(4-hydroxy-3-methoxyphenyl)-1,6heptadiene 3,5-dione]), as a typical hydrophobic polyphenol molecule extracted from Curcuma longa, was selected for investigation of its interaction with lipid bilayers in this study. CUR is widely used as a natural food coloring and seasoning from ancient times. It has two aromatic groups with hydroxyl groups and other substituents at each end of the molecule, and these two aromatic rings are linked by a diene chain, which is similar to the structure of stilbenoids. Isomerization of CUR with keto-enol form has been confirmed by NMR studies 22. As the active component of Curcuma longa turmeric, CUR was studied worldwide and reported to possess multiple functions, including: antimicrobial

23-25

, antioxidant

26-28

, anti-Alzheimer’s

29-31

, anticancer

32-33

, anti-inflammatory

34-35

etc.

Although various studies have reported that CUR is able to inhibit the growth of several microorganisms by targeting their cell membranes 36-39, understanding of the detailed interaction mechanism between CUR and cell membranes is still not adequate. Commonly used experimental techniques for directly detecting membrane disruption include: (1) measurement of fluorescence intensity of the dye that leaked from cells or liposomes

36, 38, 40

, and (2) observation of the membrane structure change through scanning electron

microscopy (SEM) 36, 41 and transmission electron microscopy (TEM) 39, 42. All-atom molecular dynamics (MD) simulation provides an efficient way to understand the interactions between molecules and lipid bilayer at an atomic level for timescales not accessible by experiments and therefore complements experimental methods 43. It is well known that cell membrane is very complex with heterogeneous architecture that contains a variety of phospholipids, cholesterol, and membrane proteins. It is also known that the lipid compositions vary with different bacterial strains, as a response to changing environment, or to exposure to other molecules

44

. Model lipid bilayer MD simulation has been widely used for many years due to its

capability of providing feasible insights for understanding the atomic interaction. Simulations of model lipid bilayer not only greatly reduced the complexity of structure construction, but also provides us clearer data interpretation to compare with experimental results. Systematic understanding of the effect of lipid properties on the interaction of molecules with lipid bilayer is important. On the other hand, studying complicated lipid bilayer is necessary to understand the function of drug in a more realistic system. MD studies on more complex lipid bilayers have been shown to describe more accurately the behavior of a 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 34

microorganism membrane (such as E. coli and yeast) than simple model lipid bilayers shown that several specific components, such as ergosterol moieties

45

47

, cardiolipin (CL)

49

45-49

. It has been

, and cyclopropane

affect the properties of membranes.

In this work, we investigated the interaction of CUR with different lipid bilayers through MD simulations: (1)

lipids

with

various

size

and

charge

on

lipid

heads:

phosphatidylcholine

(PC)

vs.

phosphatidylethanolamine (PE) vs. phosphatidylglycerol (PG), (ii) lipids with different degrees of saturation in lipid tails: 1,2-dioleoyl- (DO) vs. 1-palmitoyl-2-oleoyl- (PO) vs. 1,2-dipalmitoyl- (DP), (iii) mixed model lipid bilayer: 1-Palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE) / 1palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol (POPG) with ratio of 3:1 to mimic bacteria cell membrane, and (iv) two realistic microorganism cell membrane systems: E. coli and yeast cell membrane. The results obtained from this study will not only provide us specific insights for interaction of curcumin and different bilayer systems, but also give us a more generalized idea for choosing lipid bilayer based on simulation parameters in the future. 2. Methods 2.1 System description We conducted MD simulations of CUR in six different model membranes (consisting of pure POPE, POPG,

1-Palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine

(POPC),

1,2-Dioleoyl-sn-glycero-3-

phosphocholine (DOPC), 1,2-Dipalmitoyl-sn-glycero-3-phosphoethanolamine (DPPE), and a mixture of POPE/POPG with ratio of 3:1) and two realistic membrane systems (E. coli membrane consisting of 1,2dipalmitoyl-1’-palmytoil-2’-cis-9,10-methylenehexadecanoyl-cardiolipin

(PMCL2),

1,2-1’,2’-tetra-

hexadecenoyl-cardiolipin (TXCL2), 1,2-dipalmitoleic-phosphatidylethanolamine (DXPE/DYPE), 1,2dipalmitoyl-phosphatidylethanolamine (DPPE), and 1,2-dioleoyl-sn-glycero-3-phosphoethanolamine-N(1-pyrenesulfonyl) (PYPE); yeast membrane consisting of DOPC, POPE, 1-palmitoyl-2-oleoyl-snglycero-3-phosphate

(POPA),

1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-L-serine

(POPS),

1,2-

dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) and cholesterol). Detailed information of lipid membrane compositions was listed in Table 1. The E. coli membrane system is almost the same as that employed by Wu et al

49

, with the only difference of 1-palmytoil-2-cis-9,10-methylenehexadecanoyl-

phosphatidylethanolamine (PMPE) replaced by PYPE. Yeast membrane system contains the same components as that employed by Sunhwan et al 50, which has been reported existing similar properties compared to experimental results. The structure of CUR is shown in Figure 1. The initial CUR-membrane systems were constructed using CHARMM-GUI’s Membrane Builder module (http://www.charmmgui.org/?doc=input/membrane ) 48.

4 ACS Paragon Plus Environment

Page 5 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 1. Structure of CUR. Color code: cyan: carbon, red: oxygen, and white: hydrogen. Table 1. System information of lipid bilayers

Lipid Type POPE POPG POPC DOPC DPPE Lipid Type POPE Lipid Type PMCL2 TXCL2 DYPE DPPE PYPE Lipid Type POPA POPS POPE DOPC DPPC CHL

Model lipid bilayers Charge Tail Info. [sn-1/sn-2] Transition Temperature 51 0 16:0 / 18:1 25°C (298.15K) -1 16:0 / 18:1 -2°C (275.15K) 0 16:0/18:1 -2°C (275.15K) 0 18:1/18:1 -17°C (256.15K) 0 16:0 / 16:0 63°C (336.15K) Mixture lipid bilayer Number of lipid Lipid Type Number of lipid 96 POPG 32 49 E. coli membrane Charge Tail Info. [sn-1,sn-2/sn-3,sn-4] -2 16:0,16:0 / 16:0,16:0 -2 16:1,16:1 / 16:1,16:1 0 16:1 / 16:1 0 16:0 / 16:0 0 16:0 / 16:1 yeast membrane 48 Charge Tail Info. [sn-1/sn-2] Transition Temperature -1 16:0 / 18:1 28°C (301.15K) -1 16:0 / 18:1 14°C (287.15K) 0 16:0 / 18:1 25°C (298.15K) 0 18:1/18:1 -17°C (256.15K) 0 16:0 / 16:0 41°C (314.15K) 0 Not applicable Not applicable

Number of lipid 128 128 128 128 128 Total number of lipid 128 Number of lipid 48 18 12 6 6 Number of lipid 20 10 60 100 20 60

CUR was initially placed both parallel and perpendicular to the lipid bilayer respectively. The closest distance between CUR atom and one of the bilayer leaflets was about 1 nm (Figure 2a and 2b). For E. coli and yeast membrane system, CUR was additionally placed in transmembrane orientation (Figure 2c). Initial structure of CUR corresponding to lipid bilayer was shown in Figure 2. 3 nm minimum height of water was added on both top and bottom of the system to simulate a fully hydrated bilayer system. A 0.15

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 34

M salt (both potassium and chloride were added) concentration was set by mimicking the physiological conditions, and the final system was neutral in charge. All simulations used CHARMM36 force field for lipids

52

and ions, and CHARMM modified TIP3P water model

53-55

. The CUR parameter was generated

from CHARMM general force field (CGenFF) 56.

Figure 2. Initial position of CUR at different location of lipid bilayer. (a) Parallel, (b) perpendicular, and (c) transmembrane. Color code: cyan: carbon, red: oxygen, white: hydrogen, brown: phosphate, water molecules are not shown for clarity. 2.2 Curcumin force field validation To validate molecular model of CUR generated from CGenFF, careful benchmark calculations has been performed, which included a single CUR simulation in both water and 1-octanol and solvation free energy in both solvent, partition coefficient in water/1-octanol system (logPoctanol/water), and multiple CUR molecules in aqueous solution. For a single CUR simulation, the solute molecule was solvated separately in a box (5 × 5 × 5 nm3) containing 4409 water molecules and the same size box with 471 octanol molecules. CHARMM modified TIP3P water model 53-55 was used for water and 1-octanol from CGenFF used in octanol simulation. After the energy minimization by the method of steepest descent the system was equilibrated for 1 ns. Both systems were extended for 200 ns under constant temperature and pressure. In simulations of multiple copies of CUR, 8 or 12 or 16 molecules were randomly placed in the simulation box, containing 20448, 31488 and 40084 water molecules, respectively. Similar simulation procedures described for single CUR system were applied except that the length of production run was 100ns. Periodic boundary conditions and PME method were applied in all simulations. The force based switching functions with range of 1.01.2 nm for the LJ interactions were used. The reference temperature was set at 303.15 K using NoseHoover

57-58

extended ensemble thermostat, and the reference pressure was set at 1 atm, coupled 6 ACS Paragon Plus Environment

Page 7 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

isotropically with the Parinello-Rahman barostat 59-60. The bonds with H-atoms were constrained using the LINCS algorithm 61 The final structure of CUR in water and octanol were used as the starting structures for the calculation of solvation free energy and partition coefficient. The solvation free energy of CUR in water (∆Gwater) and octanol (∆Goctanol) were calculated using the thermodynamic integration (TI) method

62

. The coupling

parameters, λ, were scaled to 21 points (0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, and 1.00). For each value of λ, a complete workflow was conducted by steepest descents minimization, NVT equilibration (100ps), NPT equilibration (100ps), and production run (2ns). The initial 200ps production simulation was removed when calculated the solvation free energy. Partition coefficient (logPoctanol/water) was calculated from the equation 63: log ܲ௢௖௧௔௡௢௟/௪௔௧௘௥ =

∆‫ܩ‬௪௔௧௘௥ − ∆‫ܩ‬௢௖௧௔௡௢௟ 2.303ܴܶ

where ∆Gwater and ∆Goctanol were obtained from the previous calculation, R is the molar Boltzmann constant, and T is the temperature (300K in this case). All the molecular dynamics simulations were performed using GROMACS (version 5.1.1) software package 64. In addition, solvation free energy of CUR in water was also obtained from Density Functional Theory (DFT) in conjunction with Truhlar and coworker’s SMD solvation model calculations were performed using the Gaussian 09 package

66

65

. The quantum mechanical

. The equilibrium molecular geometry

under vacuum is obtained by geometry optimization, using the both M05-2X

67

and B3LYP

68-69

density

functionals with both 6-31G* and 6-31G** quality basis sets. 2.3 Molecular dynamics simulations All simulations were performed with both GROMACS (version 5.1.1)

70

and NAMD (version 2.10) 71.

Simulations performed with NAMD were used to monitor the insertion of CUR (upto 200 ns), whereas GROMACS simulation results were used for analyzing behavior of CUR inside lipid bilayer analysis (upto 1 µs). In NAMD simulation, periodic boundary conditions were applied in all directions. In all simulations, the cut-off distance for the Coulomb and van der Waals interactions was set to 1.2 nm. The force-based switching function was used for the van der Waals (Lennard-Jones, LJ) interactions with switching range of 1.0-1.2 nm, as this range was considered to be a standard to compare the lipid properties with other studies

72

. Coulomb interactions were calculated by the Particle Mesh Ewald (PME) method

73

. The

temperature was maintained at 303.15K, which is above the gel-liquid transition temperature of 298.15K for POPE, 275.15K for POPG, 275.15K for POPC, 256.15K for DOPC, and below 336.15K for DPPE. 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 34

The pressure was maintained at 1atm. The constant temperature was controlled by Langevin dynamics, and constant pressure was controlled by Nose-Hoover Langevin piston method covalent bonds were fixed using the RATTLE algorithm

76

74-75

. The hydrogen

.The simulation was performed with

minimization of 50000 steps, equilibrium of 5ns, followed by production of 200 ns with a time step of 2fs. First 100 ns of simulation were taken for analysis of insertion of CUR into the bilayer. GROMACS simulation setup was similar to NAMD simulation. Periodic boundary conditions and PME method were also applied in the simulations. The force based switching functions with range of 1.0-1.2 nm for the LJ interactions were used. The reference temperature was set at 303.15 K using Nose-Hoover 57-58

extended ensemble thermostat, and the reference pressure was set at 1 atm, coupled with semi-

isotropically using the Parinello-Rahman barostat

59-60

and a time constant of 5 ps. The bonds with H-

61

atoms were constrained using the LINCS algorithm . The simulation was performed with minimization of 50000 steps, equilibrium of 5ns, followed by production of 1000 ns with a time step of 2 fs. First 100 ns of simulation were taken for analysis of insertion of CUR into the bilayer, and the last 300 ns were used for analysis of CUR inside the bilayer. The summary of all simulation times was shown in Table 2. Table 2. Details of simulation time GROMACS CUR Orientation

POPE

POPG

POPC

DOPC

POPE/POPG

E. coli

Yeast

Parallel

1000 ns

1000 ns

1000 ns

1000 ns

1000 ns

1000 ns

1000 ns

Perpendicular

1000 ns

1000 ns

1000 ns

1000 ns

1000 ns

1000 ns

1000 ns

1000 ns

1000 ns

300 ns

300 ns

300 ns

300 ns

300 ns

300 ns

300 ns

Transmembrane Control

NAMD CUR Orientation

POPE

POPG

POPC

DOPC

POPE/POPG

DPPE

E. coli

Yeast

Parallel × 3*

200 ns

200 ns

200 ns

200 ns

200 ns

200 ns

300 ns

300 ns

Perpendicular × 2**

200 ns

200 ns

200 ns

200 ns

200 ns

200 ns

300 ns

300 ns

200 ns

300 ns 300 ns

300 ns 300 ns

Transmembrane Control

200 ns

200 ns

200 ns

200 ns

200 ns

*Three independent MD runs with different random seed; **two independent MD runs with different random seed. 2.4 Data analysis Order parameter (SCD), area per lipid (APL), and 2D lateral membrane thickness were calculated using the membrane plugin 77 implemented in VMD. The analysis accounted for lateral area occupied by curcumin, and APL was calculated through a Voronoi diagram using the qvoronoi program from the Qhull package 78

. Bilayer thickness was defined as the distance between peaks of the phosphate atoms calculated from

density profile. RDF, number of hydrogen bond, density profile, and potential energy were calculated using tools implemented in VMD79 and GROMACS. The hydrogen bond was determined based on the 8 ACS Paragon Plus Environment

Page 9 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

following criteria: (i) the cutoff for the distance Donor - Acceptor within 0.35 nm, and (ii) the cutoff for the angle Hydrogen - Donor - Acceptor within 30°. The standard error of the mean (SEM) were obtained from six independent MD simulations. 3. Results and Discussion 3.1 Simulation of curcumin in water and octanol We performed the simulation of CUR in water and octanol respectively. Several general properties of CUR (including dipole moment, radial distribution function (RDF), and dihedral angles as shown in Figure S1) obtained from water simulation were compared with the work of Ilnytskyi et al, who used OPLS united atom (OPLS-UA) force field 80. For example, the most probable dipole moment obtained from the histogram of CUR in water appeared around 3.5 D (Figure S1b), agreed well with the estimation of CUR in vacuum and water

80

. The histogram of the probability of dihedral angle PhD2 (denoted in

Figure S1 caption) showed a peak at 0° and small non-zero wings close to ±180° (Figure S1c), indicating H17 was able to point outwards to O5. The RDF for the hydrogen (H5, H17) and oxygen (O1, O5, and O2, O6) of side parts of CUR with respect to the oxygen of water (OW) were presented in Figure S1d. The peak values and positions of [H5, H17]-OW and [O2, O6]-OW were consistent with Ilnytskyi et al’s work; whereas the RDF of [O1, O5]-OW showed a local maximum (0.75) at the position of 0.28 nm, but this value was still much lower than the first peak of [O2, O6]-OW (1.41) at the position of 0.29 nm. This qualitatively comparison indicated that, O1 and O5 in our case were has higher possibility to interact with water molecules than CUR with OPLS-UA force field. The structure and solvation properties of CUR (RDF) in octanol were compared with Samanta et al’s work

81

using GROMOS96 parameters. The

distribution of distance between H9 and O4 was plotted in Figure S1e. However, in our case, there was a single peak for both water (0.18 nm) and octanol (0.17 nm) cases, indicating the internal hydrogen bond between H9 and O4 was very stable in our case. The RDF for the oxygen of side parts (O1, O2, and O5, O6) and central part (O3, O4) of CUR with respect to the oxygen atoms of the solvent molecules were shown in Figure S1f. It was obvious that the RDF values of CUR in octanol were higher than that in water indicating octanol has higher attraction to CUR. The solvation free energy of CUR was obtained with the values of -75.53 ± 2.25 kJ/mol in water and 100.24 ± 3.15 kJ/mol in octanol. The resulting partition coefficient logPoctanol/water is 4.3, which was higher than the one obtained by Samanta and coworkers(1.17) 81 using GROMOS96 force field, but much closer to other theoretically calculated value (3.07 ± 0.4)

82

or 2.517 (with CONFLEX/PM3 method)

83

. Our

result agrees with the available values, at least qualitatively, and is consistent with the known solvation property of CUR, i.e. it is more soluble in 1-octanol than water. On the other hand, CGenFF tends to

9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 34

slightly overestimate favorable interaction of CUR with the organic solvent. It is interesting to note that solvation free energies of CUR based on CGenFF model are quite negative in water and 1-octanol. Previous TI simulations reported positive solvation free energy of CUR in aqueous media using different force fields, showing opposite affinity toward water. While there are no available chemical data to verify the results from our simulations, DFT/SMD gives similar solvation free energy in water, for example, 86.5 kJ/mol with MX05-2X/6-31G* and -72.9 kJ/mol with B3LYP/6-31G*, respectively. In SMD calculations, the non-electrostatic energy indeed is positive (~12 kcal/mol in MX05-2X/6-31G*), however the electrostatics still dominates (-31 kcal/mol) and the dipole moment of the molecule is 5.3 Debye in water with the MX05-2X/6-31G* method. Other functional and basis set combinations give consistent results. Our calculations suggest that parameters for CUR from CGenFF better describe a single CUR molecule interacting with the aqueous media, compared to other force field models that display qualitatively opposite effect. CUR is also known to aggregate in water. Harza et al. reported hydrophobic hydration driven selfassembly of curcumin in water from 40 ns simulation using GROMOS 53a6 force filed

84

. Bonab and

coworkers also observed curcumin-curcumin aggregation in aqueous media by MD simulations

85

. Our

simulations of multiple CUR molecules in water qualitatively concur with these previous studies. Curcumin monomers form clusters with hydrophilic groups pointing outwards, and representative structures are shown in the Fig S1g-S1i. RDF of middle carbon atom in CUR shows a large peak at 0.55 nm for all cases (see Fig S1j-k), but we did not observe the different layers of cluster formation as Harza and his coworker. We speculate that the well-known insolubility of CUR in water is caused by its selfaggregation. CUR does have a favorable interaction with water (by those polar hydroxyl and carbonyl) but less compared to organic solvent and also itself. 3.2 Insertion of curcumin into model lipid bilayer The snapshots of simulation systems with CUR and different lipid bilayers after 1000 ns are shown in Figure 3 and Figure S2. The replicate snapshots (Figure S2) are found to be similar to the ones presented in Figure 3. We therefore discuss the behavior of CUR based on Figure 3 in the following. CUR is found to insert into most of the lipid bilayers except DPPE bilayer. Simulation results show that CUR could insert into model lipid bilayer within a few hundred nanoseconds (Table S1). DPPE bilayer (transition temperature: 336.15K) is found to be in gel phase at the simulation temperature (303.15K), as a result, the lipid tails are not flexible to accommodate CUR. Consequently, CUR cannot be inserted into DPPE bilayer up to 200 ns simulation. For other lipid bilayers, CUR stays in the lipid tails region, close to the interface of lipid head and lipid tail, which is the glycerol region. Karewicz A el al

10 ACS Paragon Plus Environment

86

observed similar

Page 11 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

phenomenon through the fluorescence quenching experiment, which is consistent with our simulation results. It is also evident that CUR could interrupt membrane and may change some properties of lipid bilayer as shown in Figure 3. For example, bending of lipid heads occurs in POPE and POPC bilayers at the local region near CUR (Figure 3a, 3c). Time evolution (only the first hundred nanoseconds are shown for clarity) of center of mass of CUR interacting with different lipid bilayers are shown in Figure S3, where the red arrows indicate insertion of CUR into lipid bilayer. The permeation of CUR into the lipid bilayer is favored and occurs within hundred nanoseconds timescale. The initial orientation has no impact on both insertion time and insertion terminal of CUR into model lipid bilayer (Table S1, S2). It is also found that DOPC and POPE/POPG bilayers are relatively easier for CUR insertion as shown in the insertion time comparison (Table S1). The number of hydrogen bonds between CUR and water, and between CUR and lipid bilayer is monitored during the insertion of CUR (Figure S4). As CUR inserts into the bilayer, the number of hydrogen bonds of CUR and water decreases, while that between CUR and lipid bilayer increases. After insertion, CUR stays inside lipid bilayer, below the phosphate lipid head, suggesting hydrophobic interaction also plays an important role during this process. We also calculated the time evolution of the positions of O1 and O5 (Figure 1) with respect to the center of lipid bilayer. The end (either O1 or O5) that first inserts into lipid bilayer in different trajectories are listed in Table S2. Not surprisingly, results show that there’s no clear preference for insertion of either end considering the symmetric structure of CUR.

Figure 3. Snapshots of CUR in different model lipid bilayers after 1000 ns simulation. DPPE is 200 ns simulation result. Color codes are the same as Figure 2. Water molecules are not shown for clarity. 3.3 Position of curcumin inside model lipid bilayer In this section, we characterize the binding position of CUR within lipid bilayer for the last 300 ns of 1000 ns simulation trajectory. The number density profiles of CUR, water, lipid head, glycerol, and terminal methyl groups as a function of distance to bilayer center are shown in Figure 4. CUR distribution is in the range of 0.5-2 nm to bilayer center, and is closer to glycerol than lipid head. This result is consistent with both earlier experimental and simulation results

86-89

. Employing high-resolution X-ray

diffraction and MD simulation, Alsop et al recently found that CUR could be inserted into 1,211 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

dimyristoyl-sn-glycero-3-phosphocholine (DMPC) bilayers and stay in the region of glycerol group region 89, consistent with our results. CUR also stays in one of the leaflets without moving through to the other leaflet, which implies existence of an energy barrier for penetration through the lipid bilayer. CUR distributions in POPE and POPE/POPG are found to be wider than the other three lipid bilayers, which indicates the movement of CUR is more flexible in POPE and POPE/POPG bilayers. The lipid head peak of POPG (located at ~1.88 nm) is the widest compared to the other lipid bilayers. The breadth of the peaks increases in the order of POPE (located at ~2.16 nm) < POPE/POPG mixed membrane (located at ~2.11 nm) < DOPC (located at ~1.93 nm) < POPC (located at ~1.84 nm). The peaks of glycerol for these lipid bilayers are following the same trend (POPE < POPE/POPG < DOPC < POPC < POPG). The widths of lipid peaks are related to rigidity of lipid as discussed below. The peak of CUR distribution with respect to the distance to bilayer center is the largest (~1.41 nm) in POPE, with this distance being similar in other bilayers (POPE/POPG: 1.19 nm, POPC: 1.14 nm, POPG: 1.23 nm and DOPC: 1.14 nm). Penetration of water molecules into lipid tail region also strongly depends on the lipid types. Figure 4 displays that water distributions extend closest to bilayer center in POPG and POPC bilayers, followed by DOPC, POPE/POPG, and POPE bilayer.

12 ACS Paragon Plus Environment

Page 12 of 34

Page 13 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 4. Number density profile of CUR and different model lipid bilayers. 3.4 Orientation of curcumin in model lipid bilayer In order to better describe the orientation of CUR in lipid bilayer, we calculated the tilt angle of CUR from the normal to the lipid bilayer for the last 200 ns out of 1000 ns simulation. The tilt angle refers to the angle between the vector (from O1 to O5 on CUR) and the normal vector to lipid bilayer. The orientation distribution of the CUR with respect to the normal vector of lipid bilayer was calculated as tilt angle as shown in Figure 5. The mean values for this angle are 80.79° (POPE), 99.25° (POPE/POPG), 85.63° (POPC), 87.54° (POPG), and 98.23° (DOPC). It is interesting to note that both the maximum and minimum tilt angles of CUR occur in lipid bilayer containing POPE.

Figure 5. Orientation distribution of curcumin in different model lipid bilayers.

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

To further compare the orientation distribution, the tilt angle was further separated into two groups: if the calculated value is 0°-60° and 120°-180°, it refers to CUR perpendicular to bilayer normal; if the tilt angle is 60°-120°, CUR takes parallel orientation to bilayer normal as shown in Figure 6a. Probability of CUR for parallel and perpendicular orientation inside different model lipid bilayers is shown in Figure 6b. The results show that CUR tends to be parallel to the lipid bilayer predominantly, which means the tilt angle of CUR is mostly from 60° to 120°. The angle distribution of CUR in POPC, POPG and DOPC are similar whereas the probabilities of parallel orientation are all above 0.8. Compared to CUR tilt angle in POPC, POPG, and DOPC bilayer systems (0.19, 0.20 and 0.17, respectively), the probabilities of perpendicular orientation in POPE and POPE/POPG lipid bilayer increase to 0.37 and 0.27, respectively.

Figure 6. Definition and comparison of CUR in different lipid bilayer. (a) Definition of tilt angle orientation of CUR according to bilayer norm. (b) Tilt angle probability of CUR inside different types of model lipid bilayer. To explain this interesting behavior, we analyzed the solvation properties of CUR in lipid bilayer as well as its interactions with different components of lipid bilayers. The RDFs of water (oxygen) and lipid molecules (oxygen of lipid head and glycerol, carbon of terminal methyl) with respect to the oxygen of CUR molecule in different lipid bilayers were shown in Figure 7. From the RDFs, it can be seen that, the probabilities of observing water in the first solvation shell around the oxygen of CUR in POPG and DOPC are relatively higher than that in the other three bilayers (Figure 7a). The first peaks of RDFs in Figure 7b and 7c are at the same position for all lipid bilayers with the peak value being in the opposite order (CUR-lipid head: POPE > POPE/POPG > POPC > POPG > DOPC; CUR-glycerol: DOPC > POPG > POPC > POPE/POPG > POPE) suggesting that this behavior may be attributed to competition between lipid head and glycerol. The first peak of CUR-glycerol are stronger than CUR-lipid head 14 ACS Paragon Plus Environment

Page 14 of 34

Page 15 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

indicating that CUR interacts more strongly with glycerol. The RDFs of CUR-terminal methyl (Figure 7d) are also different for different lipid heads even though the structures of fatty acid chains are exactly the same (except DOPC). CUR in POPE bilayer has lowest probability followed by POPE/POPG, DOPC, POPC, and POPG, which is due to the higher order parameter of POPE compared to other lipid bilayers as discussed below. The number of hydrogen bonds between CUR and water, CUR and membrane, CUR and lipid head, CUR and glycerol group were also calculated for different model lipid bilayers (Figure S5). The results show that the number of hydrogen bonds between CUR and water (~1.76-2.20) is more than that between CUR and membrane (~0.48-0.70). CUR did not form more hydrogen bonds with PE lipids than other lipid bilayers (Figure S5a), but the number of hydrogen bonds between CUR and glycerol decreased in the order of POPE, POPE/POPG, POPC, POPG, and DOPC (Figure S5b). The results suggest that the interaction between CUR and water is more pronounced than CUR with lipids, therefore, the solvation property of CUR was further investigated.

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 7. Radial distribution function g(r) of CUR with respect to different parts of solvent molecules. (a) CUR-WAT: oxygen of CUR to the oxygen of water. (b) CUR-lipid head: oxygen of CUR to the oxygen of lipid head. (c) CUR-glycerol: oxygen of CUR to the oxygen of glycerol group. (d) CUR-terminal methyl: oxygen of CUR to the carbon of terminal methyl. Since the first solvation shell peak appears at the distance of 0.36 nm from the distance of oxygen atom of CUR with respect to the oxygen atom of water molecules in all the lipid bilayers (Figure 7a), number of water molecules within 0.36 nm of CUR is calculated and compared for different types of lipid bilayers in Figure 8a. In POPE and POPE/POPG lipid bilayer, number of water molecules is higher in perpendicular orientation (3.36 and 2.94 respectively) than parallel orientation (2.69 and 2.74 respectively). While in other three bilayers, number of water molecules around CUR is higher in parallel orientation than perpendicular orientation. This implies that the solvation of CUR is dependent on the structure and 16 ACS Paragon Plus Environment

Page 16 of 34

Page 17 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

property of lipid bilayer. With the presence of POPE, solvation of perpendicular orientation is higher than parallel orientation. Different components of interaction energy between CUR and membrane, CUR and water are shown in Figure S6 and their differences between parallel and perpendicular orientations are given in Figure 8b. Positive value indicates a specific energy component is stronger for CUR in perpendicular orientation than in parallel orientation, while negative difference suggests an opposite effect. Interestingly, only POPE system shows positive differences in both short-range electrostatic and LJ interactions between CUR and water. This is consistent with the observation of higher water hydration number for CUR in perpendicular orientation in POPE bilayer. In contrast, parallel orientation is more favorable in DOPC bilayer, presumably due to the stronger CUR and water interactions when CUR aligns parallel to the membrane surface. In addition, LJ interaction of CUR and lipid in POPE system is negative (-4.44 kJ/mol), while this value becomes positive in other cases, particularly in DOPC bilayer. Such large positive LJ energy is likely related to the insertion depth of CUR in different systems. Indeed, CUR in perpendicular orientation inserts deeper than that in parallel in DOPC, while the distances between the center of mass and lipid center are similar in both orientations in POPE (Figure S7).

Figure 8. (a) Number of water within 0.36 nm of CUR at different types of model lipid bilayer. (b) Interaction energy difference between parallel and perpendicular orientation of CUR in different lipid bilayer systems. Coul-CUR-MEMB and Coul-CUR-WAT refers to the Coulombic potential energy between CUR and membrane, CUR and water respectively; LJ-CUR-MEMB and LJ-CUR-WAT refers to the Lennard-Jones potential energy between CUR and membrane, CUR and water respectively. 3.5 Effect of curcumin insertion on the bilayer thickness and area per lipid for model lipid bilayers 17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 34

Effect of CUR insertion on the bilayer thickness and area per lipid (APL) for different types of model lipid bilayers are evaluated in Figure 9. Bilayer thickness is defined as the distance between two phosphate peaks. In the absence of CUR, bilayer thickness decreases in the order: POPE/POPG (4.24 ± 0.09 nm) to POPE (4.21 ± 0.12 nm), POPC (4.01 ± 0.07 nm), DOPC (3.96 ± 0.06 nm) and POPG (3.94 ± 0.05 nm). These values are in good agreement with reported simulation results

92-93

72, 90-91

and the experimental

. Although not significant, insertion of CUR causes membrane thickness to increase in POPE

bilayer (4.23 ± 0.09 nm), whereas thinning of membrane occurred for all other bilayers (POPE/POPG: 4.11 ± 0.08 nm, POPC: 3.91 ± 0.10 nm, POPG: 3.67 ± 0.05 nm, and DOPC: 3.84 ± 0.09 nm). We also evaluated the lateral 2D membrane thickness and the corresponding density profile of CUR (Figure S8). It shows that CUR insertion could induce local membrane thinning where CUR is located, which is more obvious from the one-dimensional membrane thickness and corresponding CUR distribution at fixed location (Figure S9). Membrane thinning induced by CUR insertion has been reported previously. Huang et al 87 measured the thickness change of DOPC bilayer as a function of the CUR/lipid ratio. Their results showed that there is a non-linear membrane thinning effect by CUR. Similar results of membrane thickening were reported by Ng et al for POPE when assembled with an adenosine A2a receptor protein 94. In the absence of CUR, APL increase from POPE (0.5517 ± 0.0054 nm2) to POPE/POPG (0.5748 ± 0.0072 nm2), POPC (0.6331 ± 0.0054 nm2), POPG (0.6625 ± 0.0024 nm2) and DOPC (0.6760 ± 0.0068 nm2) (Figure 9b). Our results are consistent with reported both experimental results results

98-99

95-97

and simulation

2

. Insertion of CUR causes APL to increase (POPE: 0.5651± 0.0081 nm , POPE/POPG: 0.5848

± 0.0077 nm2, POPC: 0.6482 ± 0.0099 nm2, POPG: 0.6825 ± 0.0120 nm2, and DOPC: 0.6818 ± 0.0134 nm2), and this increase is significant in POPG lipid bilayer. Sun Y et al

100

reported the responses of

individual giant unilamellar vesicles (GUVs) to the binding of CUR from solution. Their results indicated an increase in the fractional area of DOPC in GUVs due to CUR binding, which is also consistent with our simulation results.

18 ACS Paragon Plus Environment

Page 19 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 9. (a) Bilayer thickness of different types of lipid bilayer within and without the absence of CUR. (b) Area per lipid (APL) of different types of lipid bilayer within and without the absence of CUR. ** Pvalue