Vibrational spectroscopy of water in hydrated lipid multi-bilayers. III. Water clustering and vibrational energy transfer S. M. Gruenbaum and J. L. Skinner Citation: The Journal of Chemical Physics 139, 175103 (2013); doi: 10.1063/1.4827018 View online: http://dx.doi.org/10.1063/1.4827018 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/139/17?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Predicting solute partitioning in lipid bilayers: Free energies and partition coefficients from molecular dynamics simulations and COSMOmic J. Chem. Phys. 141, 045102 (2014); 10.1063/1.4890877 Vibrational spectroscopy of water in hydrated lipid multi-bilayers. II. Two-dimensional infrared and peak shift observables within different theoretical approximations J. Chem. Phys. 135, 164506 (2011); 10.1063/1.3655671 Vibrational spectroscopy of water in hydrated lipid multi-bilayers. I. Infrared spectra and ultrafast pump-probe observables J. Chem. Phys. 135, 075101 (2011); 10.1063/1.3615717 Hydration dependent studies of highly aligned multilayer lipid membranes by neutron scattering J. Chem. Phys. 133, 164505 (2010); 10.1063/1.3495973 Effect of hydration on the structure of oriented lipid membranes investigated by in situ time-resolved energy dispersive x-ray diffraction Appl. Phys. Lett. 86, 253902 (2005); 10.1063/1.1952583

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19

THE JOURNAL OF CHEMICAL PHYSICS 139, 175103 (2013)

Vibrational spectroscopy of water in hydrated lipid multi-bilayers. III. Water clustering and vibrational energy transfer S. M. Gruenbaum and J. L. Skinner Theoretical Chemistry Institute and Department of Chemistry, 1101 University Ave., University of Wisconsin-Madison, Madison, Wisconsin 53706, USA

(Received 24 June 2013; accepted 23 September 2013; published online 4 November 2013) Water clustering and connectivity around lipid bilayers strongly influences the properties of membranes and is important for functions such as proton and ion transport. Vibrational anisotropic pumpprobe spectroscopy is a powerful tool for understanding such clustering, as the measured anisotropy depends upon the time-scale and degree of intra- and intermolecular vibrational energy transfer. In this article, we use molecular dynamics simulations and theoretical vibrational spectroscopy to help interpret recent experimental measurements of the anisotropy of water in lipid multi-bilayers as a function of both lipid hydration level and isotopic substitution. Our calculations are in satisfactory agreement with the experiments of Piatkowski, Heij, and Bakker, and from our simulations we can directly probe water clustering and connectivity. We find that at low hydration levels, many water molecules are in fact isolated, although up to 70% of hydration water forms small water clusters or chains. At intermediate hydration levels, water forms a wide range of cluster sizes, while at higher hydration levels, the majority of water molecules are part of a large, percolating water cluster. Therefore, the size, number, and nature of water clusters are strongly dependent on lipid hydration level, and the measured anisotropy reflects this through its dependence on intermolecular energy transfer. © 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4827018] I. INTRODUCTION

The hydration of lipid molecules is vital to the proper functioning of cell membranes, as water greatly influences the fluidity, stability, phase, and dynamics of lipid bilayers.1, 2 The local structure and dynamics of water surrounding biological molecules including lipids,3–5 DNA,6–9 and proteins10–12 have been extensively studied through a variety of experimental13–27 and computational28–40 techniques, and in general the behavior of such “biological” water is quite distinct from bulk water.41, 42 For example, in phosphatidylcholine (PC) lipid bilayers, water preferentially forms strong hydrogen bonds to phosphate oxygens,43 and molecular dynamics simulations have suggested a clathrate-like water cage around the PC choline.44–46 Water molecules deep within the lipid bilayer also display significantly slowed diffusive as well as reorientational motion,47–51 the latter likely due to a paucity of hydrogen-bonding acceptors.52–54 Water molecules also are preferentially oriented in space near lipid headgroups, with both anionic and zwitterionic (cationic) lipids resulting in a preference for water hydrogens (oxygens) to point towards the lipid interface.55 These changes to water dynamics and hydrogen bonding will in turn affect the properties and reactivity of biological systems. In addition to the dynamics of hydration water and lipid molecules, another important issue is the degree to which water is clustered and connected at the lipid bilayer interface. It has been suggested that the onset of biological activity in proteins is related to the hydration level at which water forms a connected, or percolating,56 network around a protein surface.57–60 Oleinikova and Brovchenko, among others, have

0021-9606/2013/139(17)/175103/10/$30.00

in addition studied the effects of water connectivity around model solutes,61 as well as the relation between water percolation to polymorphism and ion conductivity in DNA.62 Water clustering is likely also of importance in artificial membranes for use in, e.g., fuel cells, batteries, and desalination.63–71 These systems require a connected water network for proton (and ion) transport in order to function. Connected quasi-one dimensional water “wires” are also important for ion and proton transport in a number of trans-membrane proteins,72–74 as well as in water channels created in bicontinuous lipid gyroid phases.75 Thus, for lipid bilayers, our goal is to better understand the variables that affect water clustering and connectivity. As the frequency of the water OH (or OD) stretching mode is sensitive to its local environment,76 vibrational (infrared, IR) spectroscopy of this mode is a useful experimental probe of water dynamics, hydrogen bonding, and clustering.77 For neat H2 O (or D2 O), intra- and intermolecular coupling between OH (OD) chromophores results in rapid energy transfer and thus tends to complicate the interpretation of experimental spectra. However, for isotopically dilute HOD in H2 O (or D2 O), vibrational techniques such as linear and two-dimensional infrared spectroscopy (2DIR)78–87 provide information on the hydrogen-bonding environment and picosecond dynamics (spectral diffusion) of water, while pump-probe spectroscopy43, 88–91 can investigate vibrational relaxation and reorientational motion. Surfaceselective vibrational sum-frequency generation spectroscopy (SFG) has also been used to study the orientation and dynamics of water interfaces, including the interaction of water with lipid monolayers.92–100 Previously, our group has

139, 175103-1

© 2013 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19

175103-2

S. M. Gruenbaum and J. L. Skinner

utilized theoretical vibrational spectroscopy to study water in both the liquid101, 102 and solid103 phases, in small clusters,104 at interfaces,105 in reverse micelles106 and salt solutions,107 and in lipid multi-bilayers.108–110 By comparing calculations to experimental results, these studies were able to help interpret the microscopic water structure and dynamics that underlie the experimental spectra. While studies on isotopically dilute HOD are useful for investigating water dynamics, vibrational energy transfer111–113 between OH (or OD) chromophores in H2 O/D2 O mixtures can be used to investigate water clustering and connectivity. Recently, Piatkowski et al. performed anisotropic pump-probe spectroscopic experiments on the OD stretch of water in lipid multi-bilayers as a function of both the lipid hydration level and deuterium fraction.114 The decay of the anisotropy signal with time results from both the reorientation of the water OD bonds as well as energy transfer to other nearby OD chromophores. For neat D2 O, such energy transfer dominates the signal. The timescale and degree to which energy flows out of an initially excited OD stretch depend strongly on both intraand intermolecular couplings. As intermolecular couplings are sensitively distance and orientation dependent, analysis of the anisotropy decay signal provides information on the size of water clusters. In order to extract the relative positions of water molecules in lipid bilayers, Piatkowski et al. utilized a Förster energy transfer model89, 115–117 in which both intramolecular energy transfer and molecular reorientation were included separately.114 Based on fits of their model to the experimental anisotropy decays, they suggested that for all lipid hydration levels studied, most water molecules are found in small clusters with an average intermolecular distance of 3.4 Å. As the lipid hydration level increases, the density of these water clusters increases, but the average intermolecular distance does not appear to change. In this paper, we will use molecular dynamics simulations and our previously developed spectroscopic models118, 119 to calculate the anisotropic pump-probe signal for water in dioleoylphosphatidylcholine (DOPC) lipid multi-bilayers as a function of hydration level and deuterium fraction, as in the study by Piatkowski et al. Previous studies on liquid water suggest that a Förster model for intermolecular energy transfer may not be the most appropriate model for describing energy transfer in such disordered systems as hydration water around lipids.119 By comparing our calculations to the recent experiments, and by directly analyzing water clustering in our dynamics simulations, we therefore hope to provide a microscopic interpretation to the measured pump-probe anisotropy data. The organization of this paper is as follows. In Sec. II we review our methods for calculating vibrational spectra, while additional details for calculating the frequency-resolved anisotropic pump-probe signal for coupled chromophores are discussed in the Appendix. In Sec. III we compare our calculations to the experimental results of Piatkowski et al.,114 while in Sec. IV we further investigate water clustering and connectivity in our simulations. Finally, we conclude in Sec. V.

J. Chem. Phys. 139, 175103 (2013)

II. COMPUTATIONAL METHODS

In a vibrational anisotropy experiment, a chromophore interacts first with a pump laser pulse and then a probe pulse, separated by a delay time t. These pulses are polarizationresolved, and the frequency-resolved anisotropic pump-probe signal r(t; ω) is given by r(t; ω) =

S (ω, t) − S⊥ (ω, t) , S (ω, t) + 2S⊥ (ω, t)

(2.1)

where S (S⊥ ) is the transient absorption change when the pump and probe pulses are polarized parallel (perpendicular) to one another.114, 118 These absorption changes S are in turn calculated from the third-order response functions120   ∞     S ppp p (ω3 , t2 ) = Re R ppp p (t3 , t2 ) exp(iω3 t3 )dt3 , 0

(2.2) where  indicates p = p and ⊥ indicates p = p . The response   function R ppp p can be written in terms of a sum of 4-point dipole correlation functions and in general contains contributions due to ground-state bleaching (GB), stimulated emission (SE), and excited-state absorption (EA) terms.121 In the   Appendix we derive our working expression for R ppp p (t3 , t2 ), Eq. (A10), for the case of multiple, coupled chromophores and for frequency resolution near the peak frequency of the linear infrared spectrum. In order to compute Eq. (2.1) for the case of H2 O/ HOD/D2 O mixtures in lipid bilayers, we will utilize the mixed quantum/classical method employed in our previous studies of water.101, 102, 122, 123 In this model, the water OD stretching modes are treated quantum mechanically, all low-frequency modes are treated classically with molecular dynamics simulations, the water bend is ignored, and vibrational relaxation is included phenomenologically. The fluctuating frequencies, transition dipoles, and both intra- and intermolecular couplings for the OD chromophores are determined using ab initio-based electrostatic maps119, 124 at each time-step during a molecular dynamics simulation. In our calculations, we also include heterogeneous vibrational lifetimes.125–128 Each OD stretch in a D2 O molecule is assigned a lifetime of 0.4 ps, while lipid- and water-associated OD stretches in HOD molecules are assigned hydration-dependent lifetimes as discussed in Ref. 108. This procedure was recently utilized in our calculations of the 1D and 2DIR spectroscopy of isotopically dilute hydration water in lipid multi-bilayers,108, 109 and was found to yield good agreement with experimental results.43 Classical molecular dynamics simulations for water in DOPC multi-bilayers were performed as discussed previously108 for lipid hydration levels of X = 2, 4, 6, 8, and 16, where X indicates the number of water molecules per DOPC molecule. Briefly, our simulations consist of 128 lipid molecules arranged in a bilayer, with an average area per lipid of 66 Å2 for X = 16.31, 129 These bilayers are stacked in the z direction through periodic boundary conditions, and our dynamics simulations are performed in the NpT ensemble using GROMACS 3.3.1.130 Water is modeled with the SPC/E force field,131 while the united atom force field and partial

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19

175103-3

S. M. Gruenbaum and J. L. Skinner

J. Chem. Phys. 139, 175103 (2013)

0.4

a

c

b

d

rt

0.3

0.2

0.1

0.0 0.4

rt

0.3

0.2

0.1

0.0 0

1

2

3

4

t ps

5 0

1

2

3

4

5

t ps

FIG. 1. The frequency-resolved anisotropic pump-probe decay is shown for deuterium fractions of 10%, 25%, 50%, and 100% in panels (a)–(d), respectively. The solid lines indicate our calculated results, while the circles are the experimental results of Piatkowski et al.114 In each panel, the results for lipid hydration levels of X = 2, 4, 6, 12, and bulk water are indicated by the red, orange, yellow, turquoise, and purple symbols, respectively.

charges of Berger et al.132, 133 are used for the DOPC parameters. Spectra for each hydration level were calculated from multiple 100 ps simulations, with data collected every 2 fs. For each hydration level, all water molecules were simulated as H2 O rather than an H2 O/HOD/D2 O mixture. III. ANISOTROPIC PUMP-PROBE DECAY AND ENERGY TRANSFER

Using molecular dynamics simulations and the mixed quantum/classical spectroscopic model discussed in Sec. II, we have calculated the anisotropic pump-probe signal r(t), Eq. (2.1), for water in DOPC multi-bilayers as a function of the hydration level and deuterium fraction. In Figure 1, we compare our calculated results (solid lines) to the experimental data of Piatkowski et al.114 (circles). As in the experiments, we have investigated lipid hydration levels of X = 2, 4, 6, 12, and bulk water (red, orange, yellow, turquoise, and purple, respectively) for deuterium fractions of 10%, 25%, 50%, and 100% (panels (a)–(d), respectively). Our calculations are frequency resolved at the peak frequency of the linear infrared spectrum for each hydration level.108 Note that our calculations were actually performed for hydration levels of X = 8 and 16; the X = 12 result shown in Figure 1 is simply an average of these two calculations. As can be seen in Figure 1, for each deuterium concentration, the rate of decay of the anisotropy r(t) increases as the hydration level X is increased towards bulk water. At low deuterium concentrations (panel (a)), this is primarily a result of the slowed hydrogen-bond rearrangement dynamics at low lipid hydration levels.43, 51, 108 At higher deuterium fractions, however, the decay of r(t) is dominated by vibrational energy transfer, both within a water molecule as well as to nearby wa-

ters. At low hydration levels, there are typically fewer neighboring water molecules, and therefore intermolecular energy transfer can be significantly slowed compared to high hydration levels. In panels (a)–(c) of Figure 1, our calculated results are in qualitative agreement with the experiment. In each case, the theoretical r(t) for X = 12 (turquoise lines) decays somewhat more quickly than the experimental data. Interestingly, as shown in Figure 9 of Ref. 108, previous calculations of the anisotropy decay for dilauroylphosphatidylcholine (DLPC) multi-bilayers were found to decay more slowly than the experimental results of Zhao et al.43 In each case, this discrepancy for low deuterium fractions is most likely due to slightly different lipid bilayer structures and/or water dynamics in the molecular dynamics simulations, as compared to the experimental lipid multi-bilayers. Another potential source of uncertainty in our calculations involves the environmental dependence of the OD stretch vibrational lifetime. As discussed by Piatkowski et al. in Figure 5 of Ref. 114, the vibrational lifetime T1 for water OD stretches bound to lipid phosphate oxygens is typically shorter than T1 for water bound to other water molecules. In addition, the OD stretches in D2 O molecules have a significantly shorter vibrational lifetime than OD stretches in HOD molecules (0.4 versus 1.8 ps in the bulk liquid). These heterogeneous vibrational lifetimes significantly affect the calculated r(t), especially for intermediate deuterium fractions where both HOD and D2 O molecules contribute to the observed signal. For example, if vibrational lifetime effects are ignored, the calculated anisotropy decay for X = 2 (red line) in panel (c) of Figure 1 decays much more quickly than if heterogeneous lifetime effects are included. This discrepancy primarily results from the increased

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19

175103-4

S. M. Gruenbaum and J. L. Skinner

J. Chem. Phys. 139, 175103 (2013)

contribution to the signal of D2 O molecules, in which intramolecular coupling results in a rapid decay of the anisotropy r(t). Unlike in panels (a)–(c) of Figure 1, for the case of 100% D2 O shown in panel (d), our calculated results for the anisotropy decay significantly more quickly than the experimental data. This is particularly evident for the low hydration case X = 2 (red), as after a waiting time of 5 ps, the calculated r(t) has decayed to roughly 10% of its initial value, compared to 30% for the experimental data. As the anisotropy decay for a deuterium fraction of 100% primarily results from vibrational energy transfer, we wish to understand better the influence of both intra- and intermolecular energy transfer on our computations. In panel (a) of Figure 2, we thus plot the anisotropy r(t) for X = 2 and 100% D2 O, but with only certain types of coupling “turned on.” The solid line indicates the case where both intra- and intermolecular couplings are included in the calculation, while the dotted-dashed line shows the case where both couplings are turned off. This latter case corresponds to the isotopically dilute HOD in H2 O limit, and the initial drop in r(t) followed by a quite slow decay simply results from rapid librational motion of water molecules and then very slow reorientational motion due to hydrogen-bond making and breaking dynamics.51, 108 The dashed curve in panel (a) of Figure 2 corresponds to our calculated r(t) including only intramolecular couplings

0.4

rt

0.3

0.2

0.1

a

0.0 1.0 0.8

St

0.6 0.4 0.2 0.0 0.0

b 0.2

0.4

0.6

0.8

1.0

t ps FIG. 2. The (frequency-integrated) anisotropic pump-probe decay for 100% D2 O and X = 2 is plotted in panel (a) for the cases of no intra- or intermolecular couplings (dotted-dashed line), inter-molecular coupling only (dotted line), intra-molecular coupling only (dashed line), and both intra- and inter-molecular couplings (solid line). Panel (b) plots the survival probability S(t), Eq. (A14), for the same coupling cases, hydration level, and deuterium fraction as in panel (a).

within D2 O molecules. In this case, r(t) quickly drops to roughly half of its initial value and then levels off. While part of this initial drop is due to librational water dynamics, the greater part of the decay is due to rapid intramolecular energy redistribution on the time-scale of a few hundred femtoseconds. Using our previously described coupling maps,119 we have found that D2 O molecules in the X = 2 lipid multi-bilayer have an average intramolecular coupling of −35 cm−1 with a standard deviation of 7 cm−1 , compared to −38 cm−1 (7 cm−1 ) for the bulk liquid. As such, most D2 O molecules in the lipid bilayer undergo this rapid intramolecular energy transfer. The only exceptions are likely those molecules buried deep in the lipid layer that have one hydrogen-bonded and one free OD group. In this case, these frequencies of these OD chromophores are significantly off resonance, and intramolecular energy transfer is slowed. Finally, in panel (a) of Figure 2 we have also plotted r(t) including only intermolecular couplings (dotted curve). Although this curve does not initially decay as quickly as r(t) including only intramolecular couplings (dashed line), we see that vibrational energy does flow to nearby OD chromophores on multiple time-scales. An initial drop on the time-scale of several hundred femtoseconds is followed by a gradual decay on the time-scale of a few picoseconds, as energy flows out into water molecules beyond the first solvation shell. Although this intermolecular energy flow has previously been interpreted in terms of vibrational Förster energy transfer,89, 114 this theory makes several assumptions (incoherent hopping of energy between donor and acceptor, random orientation of acceptors, irreversible energy transfer, etc.) that are likely to be invalid for liquid water,119 as well as water around biological molecules. As such, while Förster theory can be used to provide a good fit to anisotropy data, care should be taken in attaching physical meaning to the associated parameters, including the Förster radius r0 .119 As the anisotropy r(t) also includes contributions from water reorientational motion, we have tried to isolate the effects of vibrational energy transfer by calculating the survival probability S(t), Eq. (A14), for X = 2 and 100% D. These results are shown in panel (b) of Figure 2 for the same coupling cases as in panel (a). Note that after 1 ps, vibrational energy has equilibrated between both OD stretches on a D2 O molecule (dashed curve), and approximately half of the initial vibrational excitation has flowed out into the first solvation shell of other water molecules (dotted curve). These time scales are similar for each hydration level X studied, though in the case of high hydration levels (or bulk water), intermolecular coupling is of even greater importance due to a complete first solvation shell of water molecules.119 Energy flow into the second solvation shell and beyond is relatively unimportant on the time-scale of the OD stretch vibrational lifetime, though this process determines the decay of the survival probability on the time scale of tens of picoseconds. IV. WATER CLUSTERING AND CONNECTIVITY

As discussed in Sec. III, the decay of the anisotropic pump-probe signal r(t) is sensitive to the size and arrangement of water clusters, as intermolecular coupling and energy

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19

175103-5

S. M. Gruenbaum and J. L. Skinner

J. Chem. Phys. 139, 175103 (2013)

100

Ρ g r nm

3

80 60 40 20 0 3

4

5

6

r FIG. 3. The water oxygen-oxygen radial distribution function g(r) is plotted for lipid hydration levels of X = 2, 4, 6, 8, 16, and bulk water (red to purple, respectively). Each curve is weighted by the water number density ρ.

transfer are strongly distance and orientation dependent. However, rather than attempt to extract water-water distances and cluster sizes from our calculated spectra, we will instead directly analyze our molecular dynamics simulations in order to understand the arrangement of water molecules in hydrated lipid multi-bilayers. In our previous studies of DLPC lipids,108, 109 we investigated the local hydrogen-bonding environment of hydration water. Briefly, as the lipid hydration level is decreased, a relatively larger fraction of water molecules becomes strongly hydrogen-bonded to lipid phosphate oxygens, as opposed to other water molecules. Thus, intermolecular energy transfer to other water OD stretching modes would likely be slowed, although the extent of this slowing should depend on whether water molecules are isolated or instead tend to form small water pools. In Figure 3 we begin to investigate this by plotting the (density-weighted) water oxygen-oxygen radial distribution function g(r) for lipid hydration levels of X = 2, 4, 6, 8, 16, and the bulk liquid (red to purple, respectively). As is evident, for each hydration level, most water molecules have a welldefined first solvation shell, and the first maximum in g(r) at r = 2.75 Å is essentially unchanged for all X. The area under this first peak (and thus the average number of water molecules in the first solvation shell) decreases dramatically, however, from bulk water to low hydration lipids. In fact, for X = 2, we find that each water molecule only has on average about 1.2 other waters within an O–O distance of 3.3 Å. At this hydration level, a small fraction of waters are indeed isolated, typically located deep within the lipid bilayer bound to carbonyl oxygens.108 However, most waters have one or two other nearby water molecules, and thus we expect some degree of intermolecular energy transfer to take place. Note that if most water molecules were located in water pools at X = 2, we would expect that g(r) for small r would more closely resemble that of the bulk liquid than is observed. In order to analyze more fully the nature of water clustering in our lipid systems, we will define two water molecules as connected (and thus belonging to the same water cluster) if either (a) they are hydrogen-bonded, or (b) if a chain of hydrogen-bonded waters connects the two molecules. For the purposes of this article, we have utilized a previously de-

FIG. 4. Snapshots from molecular dynamics simulations of a DOPC multibilayer system for lipid hydration levels of X = 2, 4, and 6 are shown. The lipids are depicted in gray, while water molecules are colored based on cluster size. Red, orange, yellow, and green indicate water clusters with N = 1, 2, 3, and 4, respectively, while shades of blue indicate larger water clusters with N ≥ 5, and purple indicates a percolating water cluster.

scribed distance-angle definition of a hydrogen bond,134 although other reasonable definitions yield similar results for water cluster sizes. From our molecular dynamics simulations at each lipid hydration level, we have then sorted each water molecule into clusters of size N. In Figure 4 we show snapshots from our simulations at X = 2, 4, and 6 (top to bottom panels, respectively). While the DOPC lipid molecules are depicted in gray, the water molecules are colored according to their cluster size. Isolated waters (N = 1) are shown in red, while orange, yellow, and green indicate cluster sizes of N = 2, 3, and 4, respectively. Larger clusters with N ≥ 5 are shown as different shades of blue, with each shade indicating a separate larger cluster. Finally, water clusters that span the simulation box (percolating clusters) are shown in purple. Note that the higher hydration levels X = 8 and 16 appear

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19

175103-6

S. M. Gruenbaum and J. L. Skinner

J. Chem. Phys. 139, 175103 (2013)

0.30

Fraction in cluster size N

0.1

0.25 0.01

0.20 0.001

0.15 2

4

6 N

0.10

8

10

0.05 0.00 2

4

6

8

10

N FIG. 5. The fraction of water molecules in clusters of size N is plotted for the same hydration levels as in Figure 2. The inset is the same data plotted on a log scale.

similar to the X = 6 snapshot, with most of the water molecules contained in a percolating (purple) water cluster. The distribution of cluster sizes is quantified in Figure 5, where we plot the fraction of water molecules in a cluster of size N for the same hydration levels as in Figure 3. The inset displays this same data on a log scale, and the case of bulk liquid water is not shown as essentially all water molecules in the liquid are part of the same percolating water network. In addition, in Figure 6 we plot the total fraction of waters that are in small (N ≤ 5) water clusters for hydration level X as the “×” symbols (the dotted line is a guide to the eye). As seen in Figures 5 and 6, and as can be visualized from the snapshots in Figure 4, the fraction of water molecules in very small clusters decreases rapidly as the lipid hydration level is increased. While nearly 30% of water molecules are isolated (N = 1) for X = 2, only about 2% of waters are isolated at X = 16. For X = 4, we observe a wide distribution of cluster sizes, including many isolated water molecules and water pairs, as well as clusters containing up to around 20 water molecules. At higher hydration levels, the minority of waters that are not part of a percolating cluster are typically those deep within the lipid bilayer, while waters surrounding the DOPC choline or hydrogen-bonded to the phosphate group

1.0

Fraction

0.8 0.6 0.4 0.2 0.0 0

4

8

12

16

X FIG. 6. The fraction of water molecules in a percolating cluster is plotted against the lipid hydration level X as the filled circles (and dashed line). The “×” symbols (and dotted line) indicate the fraction of water molecules contained in clusters of five or fewer molecules.

are usually part of the large connected network. Finally, in Figure 6 we also plot the fraction of water molecules that are part of this percolating water cluster as the colored circles (and the dashed line). This fraction rises from zero to essentially unity over the range X = 4–8, indicating a fundamental change in the nature of water clustering as the lipid hydration level is increased. For each hydration level shown in Figure 5, a fraction of water molecules are either isolated or in very small water clusters (N = 2 or 3). It is only for these very small water clusters that we would expect slow intermolecular energy transfer to take place. As soon as several water molecules are contained in a cluster, intermolecular coupling between OD chromophores results in rapid energy transfer to first solvation shell waters. As the lipid hydration level is increased, the average size of this first solvation shell increases and thus the anisotropy decays to a smaller value, but the time-scale of such energy flow is relatively unchanged. Energy flow to the second solvation shell (or out into a percolated water network) certainly occurs, but this process is on the time-scale of several picoseconds, and by this time the anisotropy signal has already mostly decayed. Therefore, our calculations suggest that the anisotropy decay is sensitive primarily to small water clusters, and measurements of r(t) cannot easily distinguish between medium, large, and percolated clusters. V. CONCLUSIONS

In this article, we have combined molecular dynamics simulations with a mixed quantum/classical model for vibrational spectroscopy in order to interpret more fully recent anisotropic pump-probe experiments114 of hydration water interacting with lipid multi-bilayers. The decay of our calculated vibrational signals is due to both reorientational motion of water molecules as well as intra- and intermolecular vibrational energy transfer. At low deuterium fractions (primarily HOD in H2 O), reorientational effects predominate, and our calculations are in good agreement with the experiments. For 100% D2 O, both intra- and intermolecular couplings between OD chromophores result in a rapid decay of the anisotropy. In this case, we find that our calculated results tend to decay more quickly and have a lower plateau value than the experiments, especially for low lipid hydration levels. This rapid decay is a result of the fact that (a) almost all water molecules undergo fast intramolecular energy redistribution, resulting in a roughly 50% drop in r(t), and (b) many (but not all) water molecules have at least one neighboring water with reasonably strong intermolecular couplings between OD chromophores. The discrepancy between theory and experiment could be the result of issues with the underlying molecular dynamics simulations or with our theoretical model, though recent testing on this model suggests that it is surprisingly robust for water in a variety of environments.124 Heterogeneous, environmentally dependent vibrational lifetimes can also significantly affect the anisotropy decay, and part of the discrepancy could be a result of an inadequate lifetime model.108 As the decay of the anisotropy is strongly dependent on water clustering (due to the distance and orientational dependence of intermolecular couplings between OD groups), we

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19

175103-7

S. M. Gruenbaum and J. L. Skinner

J. Chem. Phys. 139, 175103 (2013)

have analyzed the local environment of water molecules in our lipid multi-bilayer simulations. Unsurprisingly, we find that the average number of neighboring waters dramatically decreases as the lipid hydration level decreases. However, even for low lipid hydration levels, most water molecules have at least one other neighboring water in the first solvation shell, and the location of the first peak in the water oxygen-oxygen radial distribution function is hardly changed from bulk liquid water. At the low hydration level X = 2, up to 30% of waters are in fact isolated (mostly waters bound deep in the lipid headgroup region to carbonyl oxygens), but the remaining 70% form small clusters or connected chains and thus do exhibit intermolecular energy transfer. At higher hydration levels, only a small fraction of waters are isolated, and for X ≥ 6, most water molecules are part of a single percolating water cluster. Therefore, our simulations clearly suggest that the number and size distribution of water clusters are strongly dependent on lipid hydration level. The nature of this water clustering in turn will influence the properties of the lipid bilayer, and water connectivity will likely be important for understanding properties such as proton or ion transport in membrane materials. Overall, vibrational anisotropic pump-probe spectroscopy is a powerful tool for understanding energy transfer and clustering in complicated systems, and simulations and theoretical spectroscopy are useful for the interpretation of the experimental results.

where g, e, and f represent the ground state, set of singleexcitation states, and set of double-excitation states, respecp tively, μij is the transition dipole matrix element with polarization p = x, y, or z, and Uee and Uff are the single- and double-excitation blocks of the time propagator. Now consider the case of multiple coupled chromophores. The single-excitation states will be denoted |i , while the double-excitation states are denoted |ij , where i, j label the chromophores. Double excitation on a single chro¯ The dipole matrix elements are mophore will be denoted |i . 119 then given by  μeg = μj |j 0|, (A4) j

μf e =



μ¯ j |j¯ j | +

j



μk |j k j |.

Inserting Eqs. (A4) and (A5) into Eqs. (A1)–(A3) yields pppp

RGB (t3 , t2 ) =



 p 2 p p μj (0) μk (t2 + t3 )μl (t2 )Ukl (t2 + t3 , t2 ) ,

pppp RSE (t3 , t2 )

=



p

p

p

† × Uik (t2 , 0)Ulj (t2

pppp REA (t3 , t2 )

 + t3 , 0) ,

i,j,k,l p × k|μef (t2

+ t3 )Uff (t2 +

 .

p t3 , t2 )μf e (t2 )|l

(A8)

The rephasing third-order response function with wavevector k I = −k 1 + k 2 + k 3 is given by the sum of three terms: GB, SE, and EA. To calculate the anisotropic pumpprobe signal, we will consider the case of t1 = 0. Ignoring vibrational lifetime effects, the terms of the response function are then given by121 pppp

RGB (t3 , t2 )     = μpge (0)μpeg (0)μpge (t2 + t3 )Uee (t2 + t3 , t2 )μpeg (t2 ) , (A1)    † (t3 , t2 ) = μpge (0)Uee (t2 , 0)μpeg (t2 )μpge (t2 + t3 )  × Uee (t2 + t3 , 0)μpeg (0) , (A2)

 pppp p † REA (t3 , t2 ) = − μpge (0)Uee (t2 + t3 , 0)μef (t2 + t3 ) × Uff (t2 +

(A7)

 p p † =− μi (0)μj (0)Uik (t2 + t3 , 0)Ulj (t2 , 0)

APPENDIX: CALCULATION OF RESPONSE FUNCTIONS

pppp

p

μi (0)μj (0)μk (t2 )μl (t2 + t3 )

i,j,k,l

RSE

(A6)

j,k,l

ACKNOWLEDGMENTS

This work was primarily supported by the Department of Energy through Grant No. DE-FG02-09ER16110. We also wish to acknowledge Professor H. J. Bakker and Dr. L. Piatkowski for helpful discussions involving their experimental measurements and for providing experimental data in advance of publication.

(A5)

k=j

j

 p t3 , t2 )μf e (t2 )Uee (t2 , 0)μpeg (0) , (A3)

The excited-state absorption response function, Eq. (A8), contains matrix elements of the two-quantum propagator. These terms can be numerically computed, as described by Jansen et al.135 However, we wish to calculate the frequencyresolved anisotropic pump-probe signal near the IR peak frequency ω10 , and so we can safely ignore terms in the response function which oscillate with a frequency ω21 (that is, we will make a two-level approximation for each chromophore and ¯ We can thereignore all terms in Eq. (A5) involving |k ). fore also make a harmonic approximation to the two-quantum propagator matrix elements,135 Uab,cd (t + τ, t) 1 =√ (Uac (t + τ, t)Ubd (t + τ, t) (1 + δab )(1 + δcd ) + Uad (t + τ, t)Ubc (t + τ, t)).

(A9)

Plugging Eq. (A9) into Eqs. (A6)–(A8) and canceling terms between the GB, SE, and EA response functions119 result in

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19

175103-8

S. M. Gruenbaum and J. L. Skinner

J. Chem. Phys. 139, 175103 (2013)

In the case of uncoupled chromophores, |Uki (t, 0)|2 = 0 for k = i, and we obtain the expressions for the heterogeneous lifetime discussed in Refs. 108 and 109.

our final expression for the total response function,    p p p p μi (0)μj (0)μk (t2 + t3 )μl (t2 ) R ppp p (t3 , t2 ) = 2 i,j,k,l



†

× −Uik (t2 + t3 , 0)Ulj (t2 , 0)Ukl2 (t2 + t3 , t2 ) †

+ Uik (t2 + t3 , 0)Ukj (t2 + t3 , 0)Ukl (t2 + t3 , t2 ) 

† + Uil (t2 , 0)Ulj (t2 , 0)Ukl (t2 + t3 , t2 ) . (A10) In the limit of frequency-integrated detection (t3 = 0), Eq. (A10) reduces to Eq. (9) of Yang, Li, and Skinner.119 The one-quantum propagator matrix elements Ujk (t, t ) are propagated as119 i U˙ (t, t  ) = − H (t)U (t, t  ), (A11) ¯ where the Hamiltonian is given by   H (t) = ¯ωj (t)|j j | + ¯ωj k (t)|j k|, (A12) j

j,k

and ωj (t) and ωjk (t) are described by frequency and coupling electrostatic maps.124 The probability that an excitation initially on chromophore k is on chromophore j at time t is described by Pj k (t) = |Uj k (t, 0)|2 ,

(A13)

where the brackets indicate an average over the ensemble of trajectories. The survival probability S(t) is finally given as N 1  S(t) = Pjj (t). N j =1

(A14)

Within Fo¨rster theory, this survival probability is then related to the measured anisotropy r(t).119 In the above expressions for the third-order response functions, we have neglected vibrational lifetime effects. One reason for this is because if the lifetime is constant and independent of environment, then the vibrational lifetime terms exactly cancel in the expression for the anisotropy, Eq. (2.1). However, if the lifetime T1 depends on the environment of the chromophore, then the lifetime terms do not necessarily cancel from the anisotropy but can in fact significantly change the calculated signal. To treat this heterogeneous lifetime case, let the lifetime for each chromophore j vary with time, T1, j (t) (as the chromophore’s environment changes). If we only consider lifetime effects during the waiting time t2 (reasonable, as t3 will typically be short), then the summand in Eq. (A10) should be multiplied by a factor   t  −1 (A15) fi (t) = exp − dτ ϒ1,i (τ ) , 0

−1 where we define the time-dependent effective lifetime ϒ1,j (τ ) as  |Uki (τ, 0)|2 −1 . (A16) ϒ1,i (τ ) ≡ T1,k (τ ) k

1 A.

G. Lee, Biochim. Biophys. Acta 1666, 62 (2004). van Meer, D. R. Voelker, and G. W. Feigenson, Nat. Rev. Mol. Cell. Biol. 9, 112 (2008). 3 J. F. Nagle and S. Tristram-Nagle, Biochim. Biophys. Acta 1469, 159 (2000). 4 J. Milhaud, Biochim. Biophys. Acta 1663, 19 (2004). 5 M. L. Berkowitz and R. Vacha, Acc. Chem. Res. 45, 74 (2012). 6 S. K. Pal, L. Zhao, and A. H. Zewail, Proc. Natl. Acad. Sci. U.S.A. 100, 8113 (2003). 7 K. E. Furse and S. A. Corcelli, J. Am. Chem. Soc. 130, 13103 (2008). 8 S. Sen, D. Andreatta, S. Y. Ponomarev, D. L. Beveridge, and M. A. Berg, J. Am. Chem. Soc. 131, 1724 (2009). 9 A. Paciaroni, A. Orecchini, G. Goracci, E. Cornicchi, C. Petrillo, and F. Sacchetti, J. Phys. Chem. B 117, 2026 (2013). 10 P. Ball, Chem. Rev. 108, 74 (2008). 11 F. Sterpone, G. Stirnemann, and D. Laage, J. Am. Chem. Soc. 134, 4116 (2012). 12 D. Russo, M. A. Gonzalez, E. Pellegrini, J. Combet, J. Ollivier, and J. Teixeria, J. Phys. Chem. B 117, 2829 (2013). 13 S. R. Wassall, Biophys. J. 71, 2724 (1996). 14 H. Binder and W. Pohle, J. Phys. Chem. B 104, 12039 (2000). 15 V. Kurze, B. Steinbauer, T. Huber, and K. Beyer, Biophys. J. 78, 2441 (2000). 16 J. Sykora, P. Kapusta, V. Fidler, and M. Hof, Langmuir 18, 571 (2002). 17 S. Tristram-Nagle and J. F. Nagle, Chem. Phys. Lipids 127, 3 (2004). 18 J. Sykora, P. Jurkiewicz, R. M. Epand, R. Kraayenhof, M. Langner, and M. Hof, Chem. Phys. Lipids 135, 213 (2005). 19 K. Gawrisch, H. C. Gaede, M. Mihailescu, and S. H. White, Eur. Biophys. J. 36, 281 (2007). 20 H. Binder, Eur. Biophys. J. 36, 265 (2007). 21 V. V. Volkov, D. J. Palmer, and R. Righini, J. Phys. Chem. B 111, 1377 (2007). 22 V. V. Volkov, D. J. Palmer, and R. Righini, Phys. Rev. Lett. 99, 078302 (2007). 23 K. J. Tielrooij, D. Paparo, L. Piatkowski, H. J. Bakker, and M. Bonn, Biophys. J. 97, 2484 (2009). 24 F. Foglia, M. J. Lawrence, C. D. Lorenz, and S. E. McLain, J. Chem. Phys. 133, 145103 (2010). 25 P. Berntsen, C. Svanberg, and J. Swanson, J. Phys. Chem. B 115, 1825 (2011). 26 J. G. Davis, K. P. Gierszal, P. Wang, and D. Ben-Amotz, Nature (London) 491, 582 (2012). 27 L. Comez, L. Lupi, A. Morresi, M. Paolantoni, P. Sassi, and D. Fioretto, J. Phys. Chem. Lett. 4, 1188 (2013). 28 C. F. Lopez, S. O. Nielsen, M. L. Klein, and P. B. Moore, J. Phys. Chem. B 108, 6603 (2004). 29 A. H. de Vries, S. Yefimov, A. E. Mark, and S. J. Marrink, Proc. Natl. Acad. Sci. U.S.A. 102, 5392 (2005). 30 S. Y. Bhide and M. L. Berkowitz, J. Chem. Phys. 125, 094713 (2006). 31 N. Sapay, W. F. D. Bennett, and D. P. Tieleman, Soft Matter 5, 3295 (2009). 32 R. Vacha, S. W. I. Siu, M. Petrov, R. A. Bockmann, J. BaruchaKraszewska, P. Jurkiewicz, M. Hof, M. L. Berkowitz, and P. Jungwirth, J. Phys. Chem. A 113, 7235 (2009). 33 M. L. Berkowitz, Biochim. Biophys. Acta 1788, 86 (2009). 34 R. Vacha, M. L. Berkowitz, and P. Jungwirth, Biophys. J. 96, 4493 (2009). 35 M. Stepniewski, A. Bunker, M. Pasenkiewicz-Gierula, M. Karttunen, and T. Rog, J. Phys. Chem. B 114, 11784 (2010). 36 F. J.-M. de Meyer, A. Benjamini, J. M. Rodgers, Y. Misteli, and B. Smit, J. Phys. Chem. B 114, 10451 (2010). 37 T. Broemstrup and N. Reuter, Biophys. J. 99, 825 (2010). 38 A. Godec and F. Merzel, J. Am. Chem. Soc. 134, 17574 (2012). 39 C. Ibarguen, M. Manrique-Moreno, C. Z. Hadad, J. David, and A. Restrepo, Phys. Chem. Chem. Phys. 15, 3203 (2013). 40 M.-L. Tan, J. R. Cendagorta, and T. Ichiye, J. Am. Chem. Soc. 135, 4918 (2013). 2 G.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19

175103-9 41 D.

S. M. Gruenbaum and J. L. Skinner

Zhong, S. K. Pal, and A. H. Zewail, Chem. Phys. Lett. 503, 1 (2011). 42 M. D. Fayer, Acc. Chem. Res. 45, 3 (2012). 43 W. Zhao, D. E. Moilanen, E. E. Fenn, and M. D. Fayer, J. Am. Chem. Soc. 130, 13927 (2008). 44 H. E. Alper, D. Bassolino-Kilmas, and T. R. Stouch, J. Chem. Phys. 99, 5547 (1993). 45 K. V. Damodaran and K. M. Merz, Langmuir 9, 1179 (1993). 46 K. V. Damodaran and K. M. Merz, Biophys. J. 66, 1076 (1994). 47 T. Rog, K. Murzyn, and M. Pasenkiewicz-Gierula, Chem. Phys. Lett. 352, 323 (2002). 48 S. Y. Bhide and M. L. Berkowitz, J. Chem. Phys. 123, 224702 (2005). 49 K. Murzyn, W. Zhao, M. Karttunen, M. Kurdziel, and T. Rog, BioInterphases 1, 98 (2006). 50 H. Khandelia and Y. N. Kaznessis, Biochim. Biophys. Acta 1768, 509 (2007). 51 Z. Zhang and M. L. Berkowitz, J. Phys. Chem. B 113, 7676 (2009). 52 D. Laage and J. T. Hynes, Science 311, 832 (2006). 53 M. Ji, M. Odelius, and K. J. Gaffney, Science 328, 1003 (2010). 54 D. Laage and W. H. Thompson, J. Chem. Phys. 136, 044513 (2012). 55 S. Nihonyanagi, J. A. Mondal, S. Yamaguchi, and T. Tahara, Annu. Rev. Phys. Chem. 64, 579 (2013). 56 D. Stauffer, Introduction to Percolation Theory (Taylor and Francis, London, 1985). 57 A. Oleinikova, N. Smolin, I. Brovchenko, A. Geiger, and R. Winter, J. Phys. Chem. B 109, 1988 (2005). 58 A. Oleinikova and I. Brovchenko, Mol. Phys. 104, 3841 (2006). 59 A. Oleinikova, N. Smolin, and I. Brovchenko, Biophys. J. 93, 2986 (2007). 60 A. Oleinikova and I. Brovchenko, J. Phys. Chem. Lett. 2, 765 (2011). 61 A. Oleinikova and I. Brovchenko, Phys. Chem. Chem. Phys. 14, 5686 (2012). 62 I. Brovchenko, A. Krukau, and A. Oleinikova, Phys. Rev. Lett. 97, 137801 (2006). 63 D. E. Moilanen, I. R. Piletic, and M. D. Fayer, J. Phys. Chem. C 111, 8884 (2007). 64 M. Zhou, P. R. Nemade, X. Lu, X. Zeng, E. S. Hatakeyama, R. D. Noble, and D. L. Gin, J. Am. Chem. Soc. 129, 9574 (2007). 65 R. L. Kerr, S. A. Miller, R. K. Shoemaker, B. J. Elliot, and D. L. Gin, J. Am. Chem. Soc. 131, 15972 (2009). 66 D. B. Spry and M. D. Fayer, J. Phys. Chem. B 113, 10210 (2009). 67 M. K. Petersen, A. J. Hatt, and G. A. Voth, J. Phys. Chem. B 112, 7754 (2008). 68 R. Devanathan, A. Venkatnathan, R. Rousseau, M. Dupuis, T. Frigato, W. Gu, and V. Helms, J. Phys. Chem. B 114, 13681 (2010). 69 N. Idupulapati, R. Devanathan, and M. Dupuis, J. Phys. Chem. A 114, 6904 (2010). 70 M. J. Cheah, I. G. Kekrekidis, and J. Benziger, J. Phys. Chem. B 115, 10239 (2011). 71 R. Jorn, J. Savage, and G. A. Voth, Acc. Chem. Res. 45, 2002 (2012). 72 N. Agmon, Chem. Phys. Lett. 244, 456 (1995). 73 S. Wolf, E. Freier, M. Potschies, E. Hofmann, and K. Gerwert, Angew. Chem., Int. Ed. Engl. 49, 6889 (2010). 74 E. Freier, S. Wolf, and K. Gerwert, Proc. Natl. Acad. Sci. U.S.A. 108, 11435 (2011). 75 G. P. Sorenson, K. L. Coppage, and M. K. Mahanthappa, J. Am. Chem. Soc. 133, 14928 (2011). 76 C. J. Tainter, Y. Ni, L. Shi, and J. L. Skinner, J. Phys. Chem. Lett. 4, 12 (2013). 77 H. Bakker and J. L. Skinner, Chem. Rev. 110, 1498 (2010). 78 J. B. Asbury, T. Steinel, C. Stromberg, S. A. Corcelli, C. P. Lawrence, J. L. Skinner, and M. D. Fayer, J. Phys. Chem. A 108, 1107 (2004). 79 J. B. Asbury, T. Steinel, K. Kwak, S. A. Corcelli, C. P. Lawrence, J. L. Skinner, and M. D. Fayer, J. Chem. Phys. 121, 12431 (2004). 80 C. J. Fecko, J. J. Loparo, S. T. Roberts, and A. Tokmakoff, J. Chem. Phys. 122, 054506 (2005). 81 J. J. Loparo, S. T. Roberts, and A. Tokmakoff, J. Chem. Phys. 125, 194521 (2006). 82 J. J. Loparo, S. T. Roberts, and A. Tokmakoff, J. Chem. Phys. 125, 194522 (2006). 83 S. T. Roberts, K. Ramasesha, and A. Tokmakoff, Acc. Chem. Res. 42, 1239 (2009). 84 R. A. Nicodemus, K. Ramasesha, S. T. Roberts, and A. Tokmakoff, J. Phys. Chem. Lett. 1, 1068 (2010).

J. Chem. Phys. 139, 175103 (2013) 85 R. A. Nicodemus, S. A. Corcelli, J. L. Skinner, and A. Tokmakoff, J. Phys.

Chem. B 115, 5604 (2011). Ramasesha, S. T. Roberts, R. A. Nicodemus, A. Mandal, and A. Tokmakoff, J. Chem. Phys. 135, 054509 (2011). 87 D. E. Rosenfeld, J. Nishida, C. Yan, S. K. K. Kumar, A. Tamimi, and M. D. Fayer, J. Phys. Chem. C 117, 1409 (2013). 88 D. E. Moilanen, E. E. Fenn, Y.-S. Lin, J. L. Skinner, B. Bagchi, and M. D. Fayer, Proc. Natl. Acad. Sci. U.S.A. 105, 5295 (2008). 89 L. Piatkowski, K. B. Eisenthal, and H. J. Bakker, Phys. Chem. Chem. Phys. 11, 9033 (2009). 90 R. Costard, N. E. Levinger, E. T. J. Nibbering, and T. Elsaesser, J. Phys. Chem. B 116, 5752 (2012). 91 R. Costard, C. Greve, I. A. Heisler, and T. Elsaesser, J. Phys. Chem. Lett. 3, 3646 (2012). 92 M. Sovago, E. Vartiainen, and M. Bonn, J. Chem. Phys. 131, 161107 (2009). 93 J. A. Mondal, S. Nihonyanagi, S. Yamaguchi, and T. Tahara, J. Am. Chem. Soc. 132, 10656 (2010). 94 X. Chen, W. Hua, Z. Huang, and H. C. Allen, J. Am. Chem. Soc. 132, 11336 (2010). 95 M. Bonn, H. J. Bakker, A. Ghosh, S. Yamamoto, M. Sovago, and R. K. Campen, J. Am. Chem. Soc. 132, 14971 (2010). 96 Y. Nagata and S. Mukamel, J. Am. Chem. Soc. 132, 6434 (2010). 97 Y. Nagata and S. Mukamel, J. Am. Chem. Soc. 133, 3276 (2011). 98 Z. Zhang, L. Piatkowski, H. J. Bakker, and M. Bonn, J. Chem. Phys. 135, 021101 (2011). 99 J. A. Mondal, S. Nihonyanagi, S. Yamaguchi, and T. Tahara, J. Am. Chem. Soc. 134, 7842 (2012). 100 P. C. Singh, S. Nihonyanagi, S. Yamaguchi, and T. Tahara, J. Chem. Phys. 137, 094706 (2012). 101 B. M. Auer, R. Kumar, J. R. Schmidt, and J. L. Skinner, Proc. Natl. Acad. Sci. U.S.A. 104, 14215 (2007). 102 B. M. Auer and J. L. Skinner, J. Chem. Phys. 128, 224511 (2008). 103 L. Shi, S. M. Gruenbaum, and J. L. Skinner, J. Phys. Chem. B 116, 13821 (2012). 104 C. J. Tainter and J. L. Skinner, J. Chem. Phys. 137, 104304 (2012). 105 Y. Ni, S. M. Gruenbaum, and J. L. Skinner, Proc. Natl. Acad. Sci. U.S.A. 110, 1992 (2013). 106 P. A. Pieniazek, Y.-S. Lin, J. Chowdhary, B. M. Ladanyi, and J. L. Skinner, J. Phys. Chem. B 113, 15017 (2009). 107 Y.-S. Lin, B. M. Auer, and J. L. Skinner, J. Chem. Phys. 131, 144511 (2009). 108 S. M. Gruenbaum and J. L. Skinner, J. Chem. Phys. 135, 075101 (2011). 109 S. M. Gruenbaum, P. A. Pieniazek, and J. L. Skinner, J. Chem. Phys. 135, 164506 (2011). 110 J. L. Skinner, P. A. Pieniazek, and S. M. Gruenbaum, Acc. Chem. Res. 45, 93 (2012). 111 S. Woutersen and H. J. Bakker, Nature (London) 402, 507 (1999). 112 D. Kraemer, M. L. Cowan, A. Paarmann, N. Huse, E. T. J. Nibbering, T. Elsaesser, and R. J. D. Miller, Proc. Natl. Acad. Sci. U.S.A. 105, 437 (2008). 113 Y. Zhang, H. Fujisaki, and J. E. Straub, J. Phys. Chem. A 113, 3051 (2009). 114 L. Piatkowski, J. de Heij, and H. J. Bakker, J. Phys. Chem. B 117, 1367 (2013). 115 G. R. Fleming and G. D. Scholes, Nature (London) 431, 256 (2004). 116 M. L. Cowan, B. D. Bruner, N. Huse, J. R. Dwyer, B. Chugh, E. T. J. Nibbering, T. Elsaesser, and R. D. J. Miller, Nature (London) 434, 199 (2005). 117 A. Paarmann, T. Hayashi, S. Mukamel, and R. J. D. Miller, J. Chem. Phys. 128, 191103 (2008). 118 Y.-S. Lin, P. A. Pieniazek, M. Yang, and J. L. Skinner, J. Chem. Phys. 132, 174505 (2010). 119 M. Yang, F. Li, and J. L. Skinner, J. Chem. Phys. 135, 164505 (2011). 120 S. Mukamel, Principles of Nonlinear Optical Spectroscopy (Oxford, New York, 1995). 121 T. L. C. Jansen and J. Knoester, J. Phys. Chem. B 110, 22910 (2006). 122 J. R. Schmidt, S. A. Corcelli, and J. L. Skinner, J. Chem. Phys. 123, 044513 (2005). 123 J. L. Skinner, B. M. Auer, and Y.-S. Lin, Adv. Chem. Phys. 142, 59 (2009). 124 S. M. Gruenbaum, C. J. Tainter, L. Shi, Y. Ni, and J. L. Skinner, J. Chem. Theory Comput. 9, 3109 (2013). 86 K.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19

175103-10 125 C.

S. M. Gruenbaum and J. L. Skinner

P. Lawrence and J. L. Skinner, J. Chem. Phys. 117, 5827 (2002). Chieffo, J. Shattuck, J. J. Amsden, S. Erramilli, and L. D. Ziegler, Chem. Phys. 341, 71 (2007). 127 L. R. Chieffo, J. T. Shattuck, E. Pinnick, J. J. Amsden, M. K. Hong, F. Wang, S. Erramilli, and L. D. Ziegler, J. Phys. Chem. B 112, 12776 (2008). 128 J. L. Skinner, Theor. Chem. Acc. 128, 147 (2011). 129 D. P. Tieleman and H. J. C. Berendsen, J. Chem. Phys. 105, 4871 (1996). 130 E. Lindahl, B. Hess, and D. van der Spoel, J. Mol. Model. 7, 306 (2001). 126 L.

J. Chem. Phys. 139, 175103 (2013) 131 H.

J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987). 132 S.-W. Chiu, M. Clark, V. Balaji, S. Subramaniam, H. L. Scott, and E. Jakobsson, Biophys. J. 69, 1230 (1995). 133 O. Berger, O. Edholm, and F. Jahnig, Biophys. J. 72, 2002 (1997). 134 R. Kumar, J. R. Schmidt, and J. L. Skinner, J. Chem. Phys. 126, 204107 (2007). 135 T. L. C. Jansen, B. M. Auer, M. Yang, and J. L. Skinner, J. Chem. Phys. 132, 224503 (2010).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.117.10.200 On: Sun, 21 Jun 2015 15:06:19