Membrane and Membrane Protein Dynamics Studied ...

Membrane and Membrane Protein Dynamics Studied with Time-Resolved Infrared Spectroscopy by

Paul Stevenson M.Chem. Chemistry University of Oxford, 2011 Submitted to the Department of Chemistry in partial fulfillment of the requirements for the degree of

Doctor of Philosophy in Physical Chemistry at the

Massachusetts Institute of Technology June 2017 c 2017 Massachusetts Institute of Technology. All rights reserved.

Signature of Author: Paul Stevenson April 28, 2017 Certified by: Andrei Tokmakoff Professor of Chemistry, The University of Chicago Thesis Supervisor Accepted by: Robert W. Field Haslam and Dewey Professor of Chemistry Chairman, Department Committee on Graduate Theses

This doctoral thesis has been examined by a committee of the Department of Chemistry that included:

Professor Keith A. Nelson Thesis Committee Chair

Professor Robert G. Griffin

Professor Andrei Tokmakoff Thesis Supervisor

2

3

Membrane and Membrane Protein Dynamics Studied with Time-Resolved Infrared Spectroscopy by

Paul Stevenson Submitted to the Department of Chemistry on April 28, 2017 in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Physical Chemistry Abstract Proteins are the machinery of the cell, performing functions essential for life. Proteins do not operate in isolation, however. Their function is intimately coupled to their environment; changes in this environment modulate the behavior of the protein. One of the most striking examples of protein-environment coupling is the interaction between membrane proteins and membranes. These interactions govern some of the most fundamental processes in biology, yet the origins of protein-membrane coupling are not well understood. Infrared (IR) spectroscopy offers a route to non-invasively probing these interactions. However, despite sustained interest in the problem over many decades, only limited progress has been made using IR spectroscopy to study protein-membrane interactions. One of the main reasons for this is the density of information encoded into a small frequency range – many hundreds of oscillators may contribute to a signal which spans a 160 residues in length [28]. KcsA is a bacterial potassium channel, but shares such a strong sequence and structural similarity to many more complex potassium channels (such as the KV , KIR channels) that it is a widely used model system [29]. The structure of KcsA can be coarsely divided into three regions - the selectivity filter, a region of ≈20 residues which line the conducting pore of the channel; the pore helices, a structural scaffold which supports the selectivity filter; and the remainder of the protein, which can change conformation in response to environmental changes, such as pH [30]. The selectivity filter and pore helices are shown in Figure 1.3. KcsA selectively conducts K+ across the lipid membrane, at rates approaching 108 s−1 with a selectivity of >1000:1 between the similarly-sized K+ and Na+ ions [31] - discrimination between ions is believed to be performed by the selectivity filter. Though KcsA is not as sensitive to voltage as the potassium channels associated with voltage gating, changes in pH can be used to induce conformational changes similar to those suspected of being responsible for changes in conduction in voltage-sensing potassium channels. Though KcsA is a “model” potassium channel, few aspects of conduction through it are well-understood. Even fundamental questions, such as how many ions are present in the channel at any one time, are still actively debated [32, 33]. The potential for IR spectroscopy to provide insight into the structural changes associated with different states in KcsA is addressed in §5. Functional studies of KcsA complicate this picture further. KcsA will spontaneously enter a non-conducting state, even when experimental conditions should result in conduction [34]. This phenomenon (C-type inactivation) has been linked to regulation of excitability in cells [35]. Simulations of this process suggests that it arises from a complicated interplay between water, protein and lipid interactions [36], though detailed experimental insight into the mechanism is still lacking. Direct coupling between KcsA and its membrane environment has been suggested by several studies. Notably, the absence of anionic lipids in the membrane inhibits conduction

30

Chapter 1. Introduction

Figure 1.4: Timescales of some common fluctuation modes in lipid membranes. in KcsA [19] - binding of phosphatidylglycerol (PG) lipids appears necessary for proper function. Changes in the conduction of KcsA have also been observed when lipid membranes undergo a phase transition [37] (see 1.2.2). The similarity between KcsA and many other potassium channels suggests these membrane-dependent effects may be present in other potassium channels. Experimental insight into the origins of this interaction between KcsA and the membrane is limited, however. §5 describes efforts to develop strategies for using IR spectroscopy as a tool to gain molecular-level insight into the dynamics of KcsA, with the aim of eventually using IR spectroscopy to probe these protein-membrane interactions.

1.2.2

Membranes

Understanding the coupling between membrane proteins and their membrane environment is not possible without first understanding at least some aspects of the membrane itself. Membranes are challenging systems to study, however. One of the main reasons for this is that lipids are rather unique among the common biomolecules - they are not extended polymers (such as proteins or DNA), but they do aggregate to form large structures. These lipid assemblies are held together by very many, very weak interactions - dynamics in membranes are highly collective. The non-covalent nature of the interactions in lipid structures is one of the defining features that makes them suitable for their roles in the cell. The cell membrane must

1.2 Membranes and Membrane Proteins

31

be able to deform in response to stress (particularly important in the case of red blood cells [38]), modify their structure (such as in endocytosis) and enable embedded proteins to laterally diffuse - these processes are facilitate by having a large number of weak interactions - some of these interactions may be broken in the processes described, but this does not result in catastrophic failure of the cell. If only a few, strong interactions were present, then disrupting one of these interactions would significantly impair the integrity of the membrane. From an energy-landscape perspective, the membrane is characterized by a very rough landscape, with no clear single minimum energy configuration. Since many configurations, separated by barriers < kB T , are accessible, it is reasonable to expect that membranes are highly dynamic systems. Indeed, this is the case - fluctuations over many length- and timescales have been reported [39,40]. Some of these fluctuations are summarized in Figure 1.4 - the range of timescales spanned by these motions suggests the hierarchical dynamics discussed in §1.1.2 will be of relevance here. Thus, any properties of the membrane reported (such as density, viscosity, effective dielectric constant) represent an average value - at any point in time and space there may be a large variance in the values of these parameters. One of the fundamental concepts in this thesis is that proteins can respond not only to changes in the average value of some parameter, but are sensitive to the time-dependence of that parameter.

Membrane Structures In biological settings, lipids will preferentially form bilayer-based structures. This is driven by the amphipathic nature of the lipid - in aqueous systems, a bilayer enables the hydrophobic acyl chain to be shielded from the water, while retaining the favorable interactions between the hydrophilic headgroups and water. This bilayer structure can adopt many different morphologies, some of which are shown in Figure 1.5. One of the most common is vesicles, where the bilayer forms a hollow sphere. Vesicles come in many sizes (from < 50 nm to > 100 µm), and may be either unilamellar (a single bilayer forms a sphere,

32


Figure 1.5: Common structures of lipid membranes. Large Unilamellar Vesicles and Multilamellar Vesicles are shown on the top. A lipid monolayer at the air-water interface and stacked lipid bilayers are shown on the bottom.


33

encapsulating a water core), or multilamellar (where many vesicles are “nested” inside each other, separated by a layer of water). Note that the leaflets of a bilayer need not share the same lipid composition! The experiments presented in this thesis use symmetric leaflets to simplify analysis, but in cell membranes the inner and outer leaflets have very different compositions, effectively doubling the number of unique interactions an integral membrane protein may have with the membrane. In addition to the more cell-like vesicle structures, other experimentally-advantageous membrane structures may be formed with careful preparation. Two of these are shown in Figure 1.5 - monolayers and stacked bilayers. In particular, stacked bilayers are useful for optical spectroscopy experiments which require a high density of lipids, such as the experiments reported in §7. A lipid bilayer has three distinct regions - the hydrophobic core, the polar headgroup region, and the interfacial linking region between them. Integral membrane proteins typically span the bilayer - different regions of the protein will experience different interactions based on their position in the membrane. In this way, membranes provide more than just a structural scaffold for proteins; the heterogeneity of the membrane provides multiple protein-membrane interactions, which may be used to tune the function of the protein. The observation noted earlier (§1.2.1) that KcsA requires a membrane containing PG lipids to function is an example of this; lipid identical in every aspect other than the headgroup identity do not result in the same effect. KcsA responds specifically to the protein-headgroup interaction. One important structural aspect of lipid membranes is the packing of the lipids. This packing is extremely sensitive to temperature; lipid membranes undergo a phase transition between an ordered (gel) state and a disordered (fluid) state. A comparison of the two phases is shown in Figure 1.6. Though the changes between the two phases may seem rather minor - primarily a change in the disorder of the acyl chain - the impact on the properties of the membrane are rather significant. The lateral diffusion constant is an order of magnitude greater in the fluid phase [41], the thickness of the membrane changes by

34


Figure 1.6: The low-temperature (gel) phase and high-temperature (fluid) phase of a single-component lipid membrane.


35

≈ 20% [42, 43], and the rate of water permeation across the bilayer increases by an order of magnitude between the gel and fluid phase [44]. Biological membranes are typically in the fluid phase - some bacteria will even change the composition of their cellular membranes in response to changes in ambient temperature to ensure the membrane is fluid [45]. Though the phase transition of the membrane is not a common biological process, it does provide a convenient way to monitor how proteins respond to changes in their membrane environment. §8 utilizes this in a time-resolved manner to probe the dynamics of how gD responds to its environment.

Heterogeneity in Cell Membranes The discussion so far has focused on membranes of uniform composition. These are useful model systems; protein-lipid interactions are already challenging systems to study without the complication of multiple species of lipid. However, it is important to recognize that these model systems are not necessarily representative of the cellular environment. Real membranes are extremely heterogeneous; the plasma membrane is composed of many different types of lipid, and many proteins and small molecules, shown schematically in Figure 1.7. This picture is a representation of the fluid-mosaic model of membranes [46], one of the first models to recognize the heterogeneity in cellular membranes. There are two consequences of this heterogeneity which are of particular relevance to the experiments in this thesis. First, the spatial heterogeneity of the membrane may introduce a spatial dependence to the function of proteins; the same protein, present in different parts of the membrane, may perform different functions. This potential for spatial heterogeneity motivates experiments which aim to understand specific protein-lipid interactions; with a sufficiently large library of these interactions, it may be possible to better understand the behavior of proteins in real cellular systems. The second consequence of the heterogeneous nature of cellular membranes is the potential to fine-tune the properties (mechanical and chemical) of the membrane. A mixture of lipids and small molecules may produce a membrane with properties (e.g. spontaneous

36


Figure 1.7: Schematic diagram of the plasma membrane, containing different lipids, membrane proteins, and small molecules. curvature, fluidity, thickness) distinct from those of its components. A particularly striking example of this is cholesterol. Cholesterol is a vital part of eukaryotic membranes, comprising up to 30% of the plasma membrane. In gel-phase membranes, cholesterol decreases the order of the membrane while in fluid-phase membranes, cholesterol increases the order parameter of the membrane. The presence of cholesterol lowers the rate of lateral diffusion in fluid-phase membranes, and is believed to form small domains in the membrane. This has significant functional consequences; cholesterol content modulates the behavior of a number of ion channels. This potential for small molecules to modulate the properties of membranes is what motivates the studies in this thesis which attempt to describe protein-membrane interactions in terms of the mechanical properties of the membrane.

1.3

Experimental Strategies

Developing a picture of how membrane proteins and membrane interact, explicitly considering the many tiers of dynamics present, is no small task. Challenges, both conceptual and experimental, abound - no single experimental technique is able to fully describe the

1.3 Experimental Strategies

37

system. X-ray diffraction can provide exquisite structural information, and with the advent of serial femtosecond crystallography is no longer a purely structural technique [47]. This technique is not well-suited to characterizing fluctuations, however. At the other extreme, fluctuation correlation spectroscopy provides a unique window into the equilibrium fluctuations of a fluorescent dye (which, with careful experimental design, may be related to the dynamics of a biomolecule of interest) [48]. In between these extremes are magnetic resonance techniques which can give some insight into both structure and dynamics [49]. A deep and meaningful understanding of the interactions between proteins and membranes can only come from leveraging the strengths of many techniques. The experiments presented in this thesis focus on the development of time-resolved vibrational spectroscopy as a compliment to the techniques discussed above. Though the bulk of the following chapters is devoted to discussion of vibrational spectroscopy, it is important to acknowledge that without the structural foundation provided by X-ray crystallography and NMR (which were used to determine the structures of KcsA and Gramicidin D, respectively), it would be extremely challenging to interpret the vibrational spectra.

1.3.1

IR Spectroscopy

IR spectroscopy provides a direct route to studying the vibrations of a biomolecule. Vibrations are inherently sensitive to bonding within a molecule. Carboxylic acid groups (deprotonated) have an intense absorption at ≈1580 cm−1 . On protonation, this band disappears completely, and a new absorption appears at ≈1730 cm−1 . The formation of the O-H bond fundamentally alters the nature of the vibration; the deprotonated absorption is assigned to an asymmetric stretch of the COO− group, while the protonated absorption is assigned to the C=O stretching mode. This is shown schematically in Figure 1.8. Though the making and breaking of bonds is fundamental to biology (e.g. the synthesis and subsequent use of ATP), many important processes in biology do not alter the chemical structure of the biomolecule. An example of this is KcsA - the protein can switch between a conductive and non-conductive state without altering its chemical composition. Vibrational

38


Figure 1.8: Protonation of the carboxylic acid group fundamentally changes the normal modes of the system. This results in a a new peak appearing at a different frequency. spectroscopy can still provide insight in this case. The vibrational frequency of a (wellchosen) probe is sensitive to the environment around it - for example, the frequency of the carbonyl stretch is sensitive to the local electrostatic environment via the vibrational stark effect. In the case where there is a high-density of oscillators in the system (such as amide carbonyls in proteins [50], or O-H stretches in H2 O [51]), these vibrations will couple to one another. The resulting eigenstates will be sensitive to the structure, and structural fluctuations, of the system. These factors need not be considered in isolation taken together, they provide insight into the environmental and structural fluctuations of the biomolecule and functional group specificity. Only vibrational modes which result in a change in the dipole moment of the system can be directly observed with IR spectroscopy. Typically the change in dipole moment of the system between the ground state and first vibrationally excited state is much smaller than that of the ground state and first electronically excited state. This has the unfortunate effect (from a practical point of view) of small vibrational extinction coefficients (εvib ≈


39

5×102 M−1 cm−1 for carbonyl stretches [52], εelec ≈ 1×105 M−1 cm−1 for the dye Rhodamine 6G). The comparatively small transition dipole moments of vibrational transitions (≈0.3 D for the amide carbonyl stretch vs > 8 D for Rhodamide 6G) makes vibrations a useful reporter, however, since there is minimal perturbation to the surrounding environment on excitation of the vibration. Put another way, the vibrational Stokes shift is negligible, so the vibrational probe reports on the native dynamics of the system. Electronic excitations, though they may be detected more readily than vibrational excitations, are also more perturbative. This is a particular concern in membranes; while Stokes shift measurements can be used to study solvation dynamics, they are challenging to relate directly to the native dynamics of the solvent (in this case, the lipids).

Carbonyl Stretch Vibrations Carbonyl stretches are attractive vibrational probes of biomolecules for several reasons. They are native to a wide range of biomolecules (proteins, DNA and many lipids contain carbonyl groups), and so do not require the addition of exogenous vibrational labels. The carbonyl stretch has, for a vibration, a large transition dipole moment. As a consequence of this, the detection threshold for carbonyl stretches is lower than many other vibrations (the nonlinear IR experiments in this thesis can detect carbonyls at concentrations down to 0.1-1 mM). Finally, the carbonyl stretch has a large Stark tuning rate (≈1 cm−1 /(MV/cm)), which makes the carbonyl stretch extremely sensitive to the surrounding electrostatic environment [53]. The two carbonyl groups studied in this thesis are the backbone amide carbonyl found in proteins, and the ester carbonyl in phospholipids. These absorb at ≈1650 cm−1 and ≈1730 cm−1 respectively, though this frequency can change significantly depending on the local environment and structure. The amide carbonyl stretch involves minor contributions from nuclear motions in the amide group, and so is usually referred to as the amide I mode (other vibrational modes arising from the amide group are labeled II, III, etc.) The high density of amide carbonyls in proteins (nearest neighbors are separated by ≈3˚ A) results in

40


significant coupling between amide carbonyls. The consequence of this is that the amide I band arises from delocalized vibrations, and so is sensitive to the secondary structure of the protein (discussed below). Ester groups in phospholipids are less structured than their amide protein counterparts, and so do not display the same degree of delocalization. The dominant effect in the ester carbonyl spectrum is the local electric field. The large local electric field from the formation of a hydrogen bond makes the ester carbonyl stretch an excellent reporter on hydration levels in lipid membranes. (The amide carbonyl stretch is also very sensitive to the formation of hydrogen bonds, but interpretation of this effect is complicated by the coupling between amide I oscillators.)

Amide I Structural Sensitivity One of the richest sources of information in the amide I band is the dependence of the lineshape on the secondary structure content of the sample [54]. As noted earlier, this arises because the amide I oscillators couple to each other to form delocalized vibrational modes (calculation of these couplings is discussed further in §2). Two of the most common structural motifs, α-helices and β-sheets give distinctly different spectra. This effect arises from two factors. First, the spatial arrangement of amide groups is different in the two structures, giving rise to different non-zero coupling elements in the Hamiltonian. This in turn gives rise to different eigenstates for the systems. Of these eigenstates, only some will be observable with IR spectroscopy. This is determined by the arrangement of the dipoles of the amide I oscillators. A β-sheet is an essentially planar array of amide groups, arranged into strands (as shown in Figure 1.9). An ideal β-sheet has two IR active modes - an intense, low-frequency mode labelled ν⊥ , and a weaker high-frequency mode labelled νk [55]. The ν⊥ mode has amide carbonyls oscillate in-phase with their neighbors on adjacent strands, and out-ofphase with their neighbors on the same strand. Hence, the net transition dipole for this mode is perpendicular to the direction the strands run in (hence the name, ν⊥ ). The νk


41

Figure 1.9: IR-active vibrational modes of an ideal β-sheet. Arrows indicate the direction of the β-sheet strand (N- to C-terminus direction). Red/blue colors indicate relative phase of the oscillators. mode shows the opposite behavior - neighbors on adjacent strands are out-of-phase and neighbors on the same strand are in phase. If the amide I transition dipole moments were perfectly parallel to the carbonyl bond axis, this mode would be virtually undetectable by IR spectroscopy. However, since the transition dipole moment is not perfectly parallel to the carbonyl bond axis [56], this mode is readily observable as a high-frequency peak in the IR spectra of β-sheets. The delocalized nature of the amide I vibrations forms the core of the analysis in §5, §6 and §8. Though α-helices and β-sheets are the most common secondary structure types, they are not the only spatial arrangement of amide groups which give distinct spectral signatures. The selectivity filter of KcsA and the unusual β-helix structure of Gramicidin D give rise to unique spectral signatures. Nonlinear IR spectroscopy (see §2) can be used to determine not only the frequencies and intensities of IR active modes, but also the angle between their transition dipole moments (for the β-sheet example, θ = 90◦ ).

Time-Resolved IR Spectroscopy The linear absorption spectrum of the carbonyl stretch region of a protein or membrane already contains a considerable amount of information. As noted earlier, β-sheets and α-

42


helices have sufficiently different amide I spectra that it is possible to estimate secondary structure content simply from the shape of the amide I band [54]. Even in absorptions which do not arise from delocalized vibrational modes, the lineshape of the peak contains information on both the environmental dynamics of the vibration (specifically, the dynamics which give rise to dephasing, and thus contribute to the homogeneous linewidth) and the range of environments experienced by the vibrational probe. However, this information is densely packed in the linear spectrum - these effects cannot easily be separated. The dynamics of membranes and membrane proteins, though they may make understanding these systems more challenging, provide a route to alleviating this congestion. Two-dimensional IR spectroscopy (2D IR) measures the sub-ns vibrational correlations in the system, which enables the fastest fluctuations in the system to be directly measured. Temperature-jump (T-jump) experiments, where the system is perturbed by a rapid temperature change, provides insight into the ns-to-ms dynamics in the system. Separating the spectral responses based on their dynamics is a potentially powerful tool for interrogating the origins of these responses, and thus studying the protein-membrane interaction. The experiments in this thesis combine linear FTIR spectra with 2D IR and T-jump experiments to develop a picture of the dynamics of membranes and membrane proteins ranging from femtoseconds to milliseconds. Utilizing this vast range of experimentallyaccessible timescales, it is possible to develop the description of hierarchical dynamics discussed in §1.1.2 beginning with the very fastest tier of dynamics.

1.4

Thesis Outline

The ultimate goal of this thesis is to experimentally probe not just the protein and membrane dynamics, but the dynamics of the interactions between the protein and membrane. To achieve this aim, new approaches to IR spectroscopy (both experimental and conceptual) are required. Chapters 2 and 3 describe the practical aspects of time-resolved IR spectroscopy, and how combining different experimental implementations can provide access to timescales

1.4 Thesis Outline

43

from ps to ms. Chapter 4 describes efforts to develop new analysis tools, which are used in later chapters. The remaining chapters are devoted to experimental studies of membranes and membrane proteins. Chapter 5 describes efforts to understand the conformational behavior of the ion channel KcsA in response to the presence of different ions, and how this may be related to the conduction mechanism. Chapter 6 outlines a vibrational spectral model to describe the nonlinear spectra of different conformations of Gramicidin D. Chapter 7 describes how 2D IR can be used to determine the fluctuations of electric fields in lipid bilayers. Chapter 8 utilizes the spectral model developed for Gramicidin D to examine the structural changes in this protein in response to changes in the membrane environment. Finally, Chapter 9 discusses how a combination of 2D IR and T-jump experiments can be used to develop an understanding of the effect of small molecules on the dynamics of lipid bilayers.

44


Chapter 2

Two-Dimensional Infrared Spectroscopy 2.1

Introduction

Two-dimensional Infrared Spectroscopy (2D IR) is a spectroscopic technique which tracks the fate of a vibrational excitation in a system. Mapping these vibrations to structural coordinates in the system enables 2D IR to make “molecular movies” of the system with a time resolution fundamentally limited by the period of the vibration. (In practice, the time-resolution of the measurement is limited by the length of the mid-IR pulse used in the experiment.) The principle behind 2D IR is deceptively simple. An initial vibrational population (or coherence) is created, and the evolution of this excitation with time reveals details of the underlying Hamiltonian of the system. In practice, this is achieved by using a series of short (> τc is not sufficient to recover a pure homogeneous lineshape. The timescale τc can be recovered from analysis of the time-dependence of the lineshape, as shown in §2.4.1. Note that the subscripts p and q are missing from the above expressions. The expressions for the Kubo lineshape were originally developed with reference to stochastic fluctuations experienced by a single oscillator. This can be extended to two oscillators by reintroducing the correlation coefficient, ρpq , yielding

Cpq (t) = ρpq ∆p ∆q exp (−|t|/τpq )

(2.21)

Vibrational Relaxation Vibrational relaxation is an extremely challenging phenomenon to describe in anything more than a heuristic manner for condensed matter. One of the most challenging obstables to overcome is that, by describing the system as an isolated group of oscillators, there is nowhere for the vibrational excitation to go! A complete description of vibrational relaxation in condensed matter, however, would involve an intractably large number of states. The simplest approach to incorporate this into the description of the 2D IR spectrum is to include an exponential decay term exp(−t/T1 ) in the lineshape function. Though this

2.3 Response Functions

55

Figure 2.2: Simulated 2D IR for a single oscillator. Upper left shows the inhomogeneous broadening case, upper right shows the homogeneous limit. The bottom row shows the time-dependence of the lineshape using the Kubo lineshape function, Eq. 2.20, with a correlation time τc = 10 ps.

56

Chapter 2. Two-Dimensional Infrared Spectroscopy

might not provide any immediate insight into the nature of the relaxation, this term at least incorporates the effect of vibrational lifetime on the lineshape (which may be nonnegligible). The relaxation rate, however, is much more than a phenomenological line-broadening parameter. The vibrational relaxation rate of the system is sensitive to the dynamics of the surrounding environment, though it reports on the environment in a less direct manner than the lineshape. One widely-used model for vibrational energy transfer which demonstrates this divides the problem into “system” states of interest, coupled to a general “bath” of other modes, which are treated in some average sense. Relating this model back to the 2D IR spectrum, states considered explicitly might be the vibrational mode of interest, and some mode just outside the frequency window of the experiment. (An example of this is relaxation from the amide I carbonyl stretch to the amide II’ mode [65].) Models of this nature typically result in an expression of the form [66] Z

∞

kab ∝

dthVab (t)Vab (0)ie−iωab t

(2.22)

−∞

where kab is the energy transfer rate and Vab is the time-dependent coupling between a and b. The time-dependence of this coupling arises from the fluctuations of the bath. Thus the energy transfer rate is sensitive to the fluctuations of the environment - specifically, it is sensitive to the magnitude of the fluctuations occuring with a frequency of ωab [67].

2.4

Six-Level System

A six-level system is a useful illustration of the principles outlined in the previous section. This system describes two eigenstates a and b, with one-quantum energies ωa and ωb . The anharmonicity of each eigenstate is ∆a , ∆b . A combination state is also considered, where ωab = ωa + ωb − ∆ab . Note that this does not necessarily imply any interaction between the eigenstates - this case is described by ∆ab = 0. This six-level system is summarized in Figure 2.3. There are 28 pathways in total to consider in a six-level system (neglecting pathways

2.4 Six-Level System

57

Figure 2.3: Energy-level diagram for two coupled oscillators, a and b. ωi is the transition frequency, ∆i is the anharmonicity of a state. which involve two-quantum transitions). These are summarized in Figure 2.4. Additional pathways are possible if energy transfer or chemical exchange occurs in the system. Crosspeaks are a feature unique to multidimensional spectroscopies. This is one of the most celebrated elements of multidimensional spectroscopy; crosspeaks provide a unique insight into the relationship between different modes of the system. Crosspeaks arise because excitation of one vibrational mode in the system changes the frequency of another mode in the system - this requires that the vibrational modes be in some way anharmonic. For example, if one of the stretching modes of an isolated water molecule is excited, this will change the average O-H bond length; the bending mode of the excited water molecule will be changed. This anharmonicity is reflected in the term ∆ab - if this is zero, pathways which would give rise to oppositely-signed doublets will overlap and cancel exactly. Peaks on the diagonal also appear as a pair of oppositely-signed doublets (referred to as the GroundState Bleach, GSB, and Excited State Absorption, ESA, respectively), separated by the anhamonicity of that mode (∆a , ∆b ). This raises a rather interesting point: though for many purposes it may be convenient to use the relations for a harmonic oscillator, a truly harmonic system would give no 2D IR signal! As a final general comment, the pathways shown in Figure 2.4 do not take into account experimental considerations. The derivations here assume pulses of the form δ(t), and thus have an infinite bandwidth. There are (currently) three distinct implementations of 2D IR

58


Figure 2.4: Rephasing and Non-Rephasing pathways for a six-level system.

spectroscopy: a frequency domain implementation which uses a narrowband (< 10 cm−1 bandwidth) pump with a broadband (> 100 cm−1 bandwidth) probe [68], a time-domain implementation with all broadband pulses [59], and a time-domain implementation using a mix of broadband and continuum (> 1000 cm−1 bandwidth) pulses [69]. Pathways which involve the state |aihb| (such as the R20b0a NR pathway) are only possible if the excitation is sufficiently broad to cover the transitions of both eigenstates. An experimental comparison of the narrowband-excitation, broadband probe and uniformly-broadband excitation probe found that these differences were relatively minor [68]. However, the effect on the crosspeaks is expected to be significantly larger, particularly if the crosspeaks are being used to determine the angle between transition dipoles (discussed further in §2.4.2 and 2.4.3).

The information content of the 2D IR spectrum of even a model six-level system is vast, and has been discussed extensively elsewhere [59]. The remainder of this section is therefore devoted to highlighting the particular aspects of the 2D IR spectrum which are used in later chapters.


2.4.1

59

Center-Line Slope

Center-Line Slope (CLS) analysis is an extension of the qualitative discussion in §2.3.1 relating the time-evolution of the lineshape to the timescale of the FFCF. The CLS of a peak is obtained by determining the value of ω3 value which corresponds to the maximum of the peak for a given ω1 . This is then fit to a line, the slope of which is the CLS. (In some literature, the CLS is defined as the inverse of this slope [70].) For a single, isolated peak experiencing stochastic fluctuations, the timescale of the decay of the CLS is exactly the correlation time of the FFCF [70]. A simulated example is shown in Figure 2.5. For a purely inhomogeneous distribution the slope would be exactly unity; deviations from a slope of unity at very short waiting times are indicative of very fast fluctuations in the system (sometimes referred to as “motional narrowing”). A more challenging, and unfortunately more common, situation arises when multiple peaks are present in the system. This is the case highlighted in Figure 2.6, where each peak experiences fluctuations with different timescales. In the case where the two peaks are clearly separated, extracting each CLS decay is straightforward. However, if there is overlap between the peaks, the situation is markedly less clear [71]. For the case shown in Figure 2.6, it is not clear from visual inspection that spectral diffusion is even occurring. This case is found experimentally with the carbonyl stretch in lipids, discussed further in §7 and §9. When the number of overlapping peaks is known, it is still possible to extract information on the FFCF decays for each peak by modeling the third-order response function and optimizing the parameters against the experimentally-determined spectra.

2.4.2

Crosspeaks

Crosspeaks are a unique source of information regarding how two vibrational modes are related. There are two distinct types of crosspeaks - those which are evident even at τ2 = 0 (“zero-time” crosspeaks) and those which “grow in” with increasing τ2 . Zero-time crosspeaks indicate coupling between different vibrational modes (as discussed earlier). This by itself is often useful information; crosspeaks are one of the primary signa-

60


Figure 2.5: Simulated 2D IR spectra at various waiting times for a single oscillator experiencing stochastic fluctuations. The CLS at each time is drawn in white. The CLS decay yields the input correlation time of 2 ps.


61

tures of β-sheet secondary structure content in proteins [55] and can be used to disentangle the spectrum arising from a mixture of structures [72]. As discussed below (§2.4.3), the intensity of crosspeaks can be used to determine the spatial relationship of two coupled vibrational modes. Time-dependent crosspeaks can arise from one of two sources. Either the excitation created moves to a new site (Vibrational Energy Transfer, discussed in §2.3.1), or the excitation stays where it is, and the environment around the excitation changes (Chemical Exchange) [73]. VET and CE are spectrally identical - they cannot be distinguished by 2D IR alone. In this case, (such as in §7), additional experiments (such as partial isotopic substitution) must be performed. Partial isotopic substitution can confirm the presence of absence of VET; if crosspeaks grow in between the natural isotope and the substituted isotope, this must be from VET. However, this does not rule out CE as the effects are not mutually exclusive.

2.4.3

Polarization Control

One of the most fundamental differences between 2D IR and FTIR measurements is the multiple interactions with different pulses in a 2D IR experiment. After the first interaction, the sample is no longer isotropic; transition dipoles lying parallel to the polarization vector of the incident field are preferentially excited. Carefully controlling the polarization of the different fields can therefore be used to determine spatial relationships and dynamics in the system [60]. This information is contained within the orientational response function, Ynabcd . Each pathway in the response function has an orientational term associated with it, which depends on both the polarization of the electric fields and the angle between the dipoles of different vibrational modes involved in that pathway. These terms are tabulated extensively elsewhere [60], so only a few key results are summarized here. The two most common polarization schemes used in 2D IR are ZZZZ, where all fields have the same polarization, and ZZYY, where the first two fields to interact with the sample

62


Figure 2.6: Simulated 2D IR spectra at various waiting times for two oscillators experiencing stochastic fluctuations. The CLS at each time is drawn in white. The top row illustrates the case where the two peaks are clearly separated. The bottom row shows the case where both peaks overlap significantly.

2.5 Carbonyl Vibrations

63

have polarization vectors perpendicular to the third interaction and detected signal field. The relative intensities of peaks in the 2D IR spectrum will change with the polarization scheme. For a zero-time crosspeak, the depolarization ratio ρ = IZZY Y /IZZZZ can be used to determine the angle between the two vibrational modes:

ρ=

7 − cos2 θ 12 cos2 θ + 6

(2.23)

where θ is the angle between the transition dipoles of the two vibrational modes involved in the crosspeak. This relationship is used in §6 to determine the angle between two vibrational modes of a protein. Note, however, that this expression only holds for the time-domain allbroadband implementation of 2D IR discussed in §2.4. The time-dependent behavior of the depolarization ratio (or alternatively, the anisotropy r(t) = (IZZZZ − IZZY Y )/(IZZZZ + 2IZZY Y ) ) provides insight into the dynamics of the system. The simplest case is physical reorientation of the molecule, which will scramble the initial anisotropy created by the interaction with the first pulse. Though this has been observed in small molecules, this effect is negligible for the larger systems studied in this thesis. The rotational correlation time of protein-sized molecules is ≈ 10 ns [74], much longer than the timescales typically accessible with 2D IR. Decay of the depolarization ratio faster than physical reorientation can occur is a sign that there is some energy exchange between modes with different transition dipole vectors.

2.5

Carbonyl Vibrations

The description of 2D IR developed in the previous section outlines how to interpret, or even reconstruct, the spectrum if the eigenstates are known. This sidesteps one of the central questions of vibrational spectroscopy - what are the eigenstates of the system? There is no simple, general answer to this question. However, for the systems considered in this thesis, it is possible to describe a method for determining the eigenstates (in principle, even if it is not always computationally tractable).

64


Figure 2.7: Left, schematic of an amide unit from a protein. The carbonyl stretch transition dipole is shown in green; this is rotated away from the C=O bond axis towards the N by 20◦ . Right, schematic of a site basis amide I Hamiltonian. Site energies are shown in red on the diagonal, coupling terms in blue.

The vibrations of interest here are carbonyl stretches. These vibrations are both commonly found in biology (as amide or ester groups) and have a large transition dipole moment (compared to other vibrations found in biomolecules). Though typically described as a carbonyl stretch, the vibrational mode does in fact include minor contributions from surrounding atoms - hence amide and ester groups are observed at different frequencies. The transition dipole moment of these carbonyl stretches lies approximately along the carbonyl bond axis; in the case of the amide vibration (usually referred to as amide I), the transition dipole is rotated ≈ 20◦ towards the nitrogen in the amide group [56]. This is shown schematically in Figure 2.7.

The first simplifying approximation is that the carbonyl groups are isolated from other vibrations - they form an isolated manifold. With this assumption, it is possible to construct a Hamiltonian H0 (this is the same Hamiltonian as discussed in §2.3) in the basis of carbonyl stretches, and diagonalize to determine the eigenstates. This site-basis Hamiltonian is composed of two elements: site energies and coupling between carbonyl groups.

2.5 Carbonyl Vibrations

2.5.1

65

Site Energies

The primary factor influencing the site energy of the carbonyl stretch is a Stark shift from local electrostatic effects. The simple nature of the carbonyl stretch is an advantage here (in contrast to more complex vibrations, such as ring-breathing modes in nucleic acids) because interpretation of the Stark shift becomes more straightforward. The carbonyl frequency will shift in a manner proportional to the projection of the local electric field of the system onto the carbonyl bond axis (or slightly off the carbonyl bond axis in the case of amide groups) [75], i.e. ωi = ω0 + cµi · E

(2.24)

Determination of this proportionality constant is one of the key goals of spectroscopic “maps” [76], which take an input such as a MD simulation and endeavour to generate a spectroscopic Hamiltonian from that simulation, such as in §5. The formation of a hydrogen bond to a carbonyl creates a large local electric field (≈ 10-20 MV/cm), which in turn results in a significant frequency shift (≈ 16 cm−1 ) to lower frequencies. In many biomolecules, the presence of absence of a hydrogen bond is the primary determinant of the site energy. Several early models of amide I spectroscopy captured many details of the system considering only whether the carbonyl was involved in a hydrogen bond when determining site energies [77]. Several systems in this thesis (§5 and §7) involve proximity of the carbonyl groups to charged residues or ions, so the Stark shift must be explicitly considered. The Stark tuning rate for carbonyls varies between 0.51.5 cm−1 /(MV/cm) [53]. For amide I vibrations, a spectroscopic map has been developed which can quantitatively predict site energies from MD simulations [76].

2.5.2

Couplings

Models for how carbonyl stretches couple to each other are less well-developed than siteenergy maps. The simplest model considers the carbonyl stretches interacting simply as

66


two dipoles. This gives rise to coupling interactions of the form

(µi · rij )(µj · rij ) µi · µj −3 Jij = A 5 |Rij | |Rij |3

(2.25)

where µi is the transition dipole of site i, rij is the unit vector connecting sites i and j and Rij is the distance between sites i and j. Values for the coefficient A of 400-600 cm−1 have been used to obtain satisfactory agreement between simulated and experimental spectra [55]. This form of coupling is typically referred to as point-dipole coupling (PDC), transitiondipole coupling (TDC) or the floating-oscillator model. One known deficiency of this model is that it does not capture “mechanical” effects. Carbonyls connected by only a small linker (such as neighbouring peptide units) will interact through the linker, not just through dipole interactions. Density-Functional Theory (DFT) calculations of small peptides in different conformations have been used to determine values for the coupling of nearest-neighbor groups as a function of the Ramachandran angles φ and ψ [78].

2.5.3

Constructing a Hamiltonian

Armed with expressions for the site energies and coupling between carbonyl groups, it is possible to construct a spectroscopic Hamiltonian, and use this to determine the eigenstates of the system. The one-quantum Hamiltonian is straightforward to construct

Hij = ωi δij + (1 − δij )Jij .

(2.26)

Combined with site-basis dipole vectors, this Hamiltonian is sufficient to interpret (or calculate) FTIR absorption spectra. To understand 2D IR spectra, however, consideration of the two-quantum Hamiltonian is required [64]. Combination bands (such as those which give rise to crosspeaks) cannot exist in a one-quantum description. Construction of this two-quantum Hamiltonian requires some additional assumptions. The first of these is that there is no coupling between the one- and two-quantum manifolds - the Hamiltonian is

2.6 Conclusions

67

block diagonal. For the two-quantum states, the site energies are taken to follow the form ωmn = ωm +ωn −δmn ∆, where ∆ is the anharmonicity (found to be consistently 16 cm−1 for amide I vibrations). Coupling in the two-quantum manifold is taken to be of the form √ Jmn,mq = 1 + δmn Jnq . Thus for a model two-oscillator system, the Hamiltonian becomes   0      ωa Jab        Jab ωb  . H0 =  √     2J 2ω − ∆ 0 a ab     √  0 2ωb − ∆ 2Jab      √ √ 2Jab 2Jab ωa + ωb

(2.27)

Diagonalizing this Hamiltonian yields the eigenstates used in the response functions. If these parameters are obtained from a simulation, a Hamiltonian may be constructed at each time step, yielding a time-dependent Hamiltonian. A time-dependent Hamiltonian is necessary to capture effects such as CLS decay and energy transfer [65].

2.6

Conclusions

The response of a system to incident infrared electric fields has the potential to describe not only the structure of the system, but also the dynamics. In order to interpret the linear and nonlinear IR spectra, however, an understanding of the factors which give rise to the eigenstates of the system is necessary. For the amide I vibration in proteins, these eigenstates can be calculated by considering the system as a set of coupled oscillators, each with a site energy which reports on the local environment of that oscillator.

68


2.A

Lineshape Functions

The general expressions for the lineshape functions (without vibrational relaxation), F abcd may be expressed in terms of the frequency-frequency correlation functions of the system. − ln F1abcd =hbb (τ3 ) + hcc (τ2 ) + hdd (τ1 ) + hbc (τ3 + τ2 ) − hbc (τ3 ) − hbc (τ2 ) +hcd (τ2 + τ1 ) − hcd (τ2 ) − hcd (τ1 ) + hbd (τ3 + τ2 + τ1 ) − hbd (τ3 + τ2 )

(2.28a)

−hbd (τ2 + τ1 ) + hbd (τ2 ) (2.28b) − ln F2abcd =hcc (τ3 ) + hbb (−τ2 ) + hdd (τ1 + τ2 + τ3 ) + hbc (−τ2 − τ3 ) − hbc (−τ2 ) +hbc (−τ3 ) + hcd (τ1 + τ2 ) − hcd (τ3 + τ2 + τ1 ) − hcd (−τ3 ) + hbd (τ1 )

(2.28c)

−hbd (−τ3 − τ2 ) − hbd (τ2 + τ1 ) + hbd (−τ3 ) (2.28d) − ln F3abcd =hbb (−τ3 ) + hcc (τ3 + τ2 ) + hdd (τ1 ) + hcd (τ3 + τ2 + τ1 ) − hcd (τ3 + τ2 ) −hcd (τ1 ) − hbc (τ3 + τ2 ) − hbc (−τ3 ) + hbc (τ2 ) − hbd (τ3 + τ2 + τ3 )

(2.28e)

+hbd (τ1 + τ2 ) + hbd (τ3 + τ2 ) − hbd (τ2 ) (2.28f) − ln F4abcd =hcc (τ3 ) + hdd (τ2 + τ1 ) + hbb (−τ3 − τ2 ) + hbc (−τ2 ) −hbc (τ3 ) − hbc (−τ2 − τ3 ) + hcd (τ1 + τ2 + τ3 ) − hcd (τ1 + τ2 ) − hcd (τ3 )

(2.28g)

+hbd (τ1 ) − hbd (τ3 + τ2 + τ1 ) − hbd (−τ2 ) + hbd (τ3 ) where Z hpq (t) =

t

Z

0

s2

ds1 Cpq (s2 − s1 )

ds2

(2.29)

0

and the frequency-frequency correlation function Cpq is defined as

Cpq (s2 − s1 ) = hωp (s2 )ωq (s1 )i.

(2.30)

2.1 Lineshape Functions

69

If the system of interest is only weakly anharmonic, then the harmonic scaling rules may be used to reduce the number of expressions:

Ca,2a = 2Ca,a C2a,2a = 4Ca,a Ca,ab = Ca,a + Ca,b Cab,ab = Ca,a + Cb,b + 2Ca,b

For a six-level system, these relations may be used to reduce the number of lineshape terms: Rephasing F30,a,0,a = F40,a,0,a = F20,a,2a,a F30,b,0,b = F40,b,0,b = F20,b,2b,b

GSB/ESA diagonal, ω1 = ωa , ω3 = ωa GSB/ESA diagonal, ω1 = ωb , ω3 = ωb

F30,b,0,a = F20,a,ab,a

GSB/ESA crosspeak τ2 population, ω1 = ωa , ω3 = ωb

F30,a,0,b = F20,b,ab,b

GSB/ESA crosspeak τ2 population, ω1 = ωb , ω3 = ωa

F40,b,0,a = F20,a,ab,b

GSB/ESA crosspeak τ2 coherence

F40,a,0,b = F20,b,ab,a

GSB/ESA crosspeak τ2 coherence

NonRephasing F20,a,0,a = F10,a,0,a = F40,a,2a,a F20,b,0,b = F10,b,0,b = F40,b,2b,b

GSB/ESA diagonal, ω1 = ωa , ω3 = ωa GSB/ESA diagonal, ω1 = ωb , ω3 = ωb

F20,a,0,b = F40,b,ab,a

GSB/ESA diagonal τ2 coherence, ω1 = ωa , ω3 = ωa

F20,b,0,a = F40,a,ab,b

GSB/ESA diagonal τ2 coherence, ω1 = ωb , ω3 = ωb

F10,b,0,a = F40,a,ab,a

GSB/ESA crosspeak, ω1 = ωa , ω3 = ωb

F10,a,0,b = F40,b,ab,b

GSB/ESA crosspeak, ω1 = ωb , ω3 = ωa

70


Chapter 3

2D IR Spectrometer 3.1

Introduction

A 2D IR spectrometer is comprised of many different elements, which can be assigned to one of three functions - pulse generation, manipulation or detection. Mid-IR pulses must be created, then spatially and temporally overlapped in such a way as to generate signal, which must then be detected. This is summarized in Figure 3.1, which shows schematically the layout of the 2D IR spectrometer used to collect the data presented later. Each of these stages is discussed in detail later, but in brief, a Regenerative Amplifier produces 90 fs pulses of light centered at a wavelength of 800 nm, which are in turn used to pump an Optical Parametric Amplifier (OPA). The outputs of this OPA are used to generate mid-IR light at the difference frequency of the two OPA outputs in a Difference Frequency Generation (DFG) setup. The mid-IR pulse is then sent through a 5-beam interferometer which generates three pulse replicas in a boxcar geometry, which can be focused down to a single spot using an off-axis parabolic mirror. The time-delays between these pulses are controlled with retroreflective optics and motorized stages. The remaining two beams from the five beam interferometer are used for heterodyne detection and alignment purposes respectively. After the beams are focused in the sample, a signal field is emitted in a unique direction. 71

72

Chapter 3. 2D IR Spectrometer

Figure 3.1: Schematic layout of the transient 2D IR spectrometer. The output of a regenerative amplifier is converted to mid-IR frequencies by a two-stage downconversion process involving Optical Parametric Amplification (OPA) and Difference-Frequency Generation (DFG). This signal field is overlapped with a weak reference beam (known as the Local Oscillator, LO) for heterodyne detection. The heterodyned signal is then dispersed by a grating on to a Mercury-Cadmium-Telluride (MCT) array detector. The spectrometer can be extended to study transient triggered phenomena by coupling an additional laser into the sample. The most common implementation of this is using an intense nanosecond 2 µm wavelength pulse to rapidly increase the temperature of the sample.

3.2 3.2.1

Mid-IR Generation Regenerative Amplifier

A Regenerative Amplifier is a system which amplifies a weak input pulse by passing it through an excited medium many times, controlling the number of passes with Pockels cells and polarizers. The seed for the amplifier is produced by an oscillator (Vitesse, Coherent Inc.) which outputs a pulse train of ωs > ωi ) we must begin with an expression for the field which is a superposition of these fields:

E(z, t) = Ep ei(ωp t+kp z) + Es ei(ωs t+ks z) + Ei ei(ωi t+ki z) + c.c.

(3.7)

Note that this expression for the field has not introduced any artificiallity into our solution; if no signal field is present, the amplitude of this component will be zero. The term E 2 contains 36 separate terms. We will show later that with careful experimental design, we can minimize the number of these terms which are observed in an experiment. For now, we will consider the second-order polaization to contain only the terms:

P (2) = ε0 def f Ep Ei∗ ei(ωs t−(kp −ki )z) + Ep Es∗ ei(ωi t−(kp −ks )z) + Es Ei ei(ωp t−(ks +ki )z) + c.c (3.8)

where we have used the relation ωp = ωs + ωi to identify the terms which oscillate at the frequencies of our three mixed fields. Inserting Equation 3.8 into Equation 3.6 and collecting terms yields a set of three coupled equations:

def f ωs ∂Es = −i Ep Ei∗ e−∆kz ∂z ns c def f ωi ∂Ei = −i Ep Es∗ e−∆kz ∂z ni c def f ωp ∂Ep = −i Es Ei e∆kz ∂z np c

(3.9) (3.10) (3.11)

where ∆k = kp −ks −ki . These equations show how the different frequencies are nonlinearly coupled; if two fields are injected into the nonlinear medium, the third will be generated. We can solve these equations under the (usually reasonable) assumption that no signif-

3.2 Mid-IR Generation

77

icant depletion of the pump field occurs (Equation 3.11 is equal to zero). The solutions to these equations are:

Is (z) = Is (0) cosh gz Ii (z) =

k1 Is (0) sinh gz k2

(3.12) (3.13)

where g=

p Γ2 − (∆k/2)2

Γ2 =

d2ef f ωs ωi c2 ns ni

|Ep |2

(3.14) (3.15)

We have no idler field and some arbitrary level of signal field as our initial conditions. The difference between OPG, OPA and DFG lies in the magnitude of the initial signal field present. OPG has no signal input initially, only the pump field. Though Equations 3.12 and 3.13 would suggest this should result in no signal and idler generations, effects beyond the scope of our classical description (such as vacuum fluctuations) give rise to some small initial effective Is (0). OPG is generally only practical for short pulses which are able to produce very high instantaneous electric fields, and finds more application in quantum optics (where low-photon number, entangled sources are of great interest) than in spectroscopy. OPA has a weak signal field at input, which acts as a seed to convert energy from the pump field into the signal and idler fields. The practical aspects of the generation of this seed field are discussed in Section 3.2.3. DFG is the case where the pump and signal fields have comparable intensity initially. This has the effect of slightly amplifying the signal field and greatly increasing the idler field magnitude. Often DFG schemes (such as the one presented in Section 3.2.4) use the signal and idler outputs from one OPA stage as the pump and signal inputs for the DFG stage.

78

Chapter 3. 2D IR Spectrometer One of the key terms in the above expressions is ∆k - this is known as the phase-

mismatch. To illustrate the importance of this term for experimental design, consider the reverse process of OPG, Sum-Frequency-Generation. Under the assumption of undepleted signal and idler fields, Equation 3.11 can be integrated between z = 0 and z = L to yield:

Ip (L) =

d2ef f ωp2 n2p c2

Is Ii L2 sinc

∆kL . 2

(3.16)

Clearly the optimal case for SFG (and other nonlinear processes) is ∆k = 0. This condition is challenging to meet in isotropic media. In a birefringent crystal, the multiple axes with different refractive indices can be used to tune the degree of phase-mismatch by varying the angle of the crystal with respect to the direction of propagation if the the three fields do not have the same polarization. This phase-matching requirement justifies our selection of only a few of the 36 terms in the second-order polarization - by careful choice of crystal and polarization conditions, many of these terms may be neglected. Phase-matching is also exploited in 2D IR experiments where non-collinear beams incident on the sample result in signals corresponding to different processes being generated in different directions. At this point, we risk overlooking another key experimental design principle because of our assumption earlier of monochromatic plane waves. Our experiments use pulsed sources; the phase-matching constraint is typically only strictly satisfied for one frequency; frequency-components around this central frequency will experience a degree of phasemismatch, which limits the bandwidth of the pulse which can be generated by OPA. Equation 3.16 shows that this effect will become more pronounced with longer crystals.

3.2.3

Optical Parameteric Amplifier

The experimental implementation of the OPA process discussed above involves three sequential nonlinear stages. First, a weak continuum pulse is generated to use as a signal seed. This is achieved by focusing a small amount of the 800 nm wavelength output from the Regenerative Amplifier into a Sapphire plate. The field undergoes self-phase modulation


79

Figure 3.2: Schematic layout for the TOPAS OPA. Focusing optics have been omitted for clarity. Linewidth represents the relative intensity of each beam. Note that before whitelight generation, a waveplate is present. This results in the output signal polarization being orthogonal to the pump polarization.β-BBO - β-Barium Borate crystal, BS - Beamsplitter, DCM - Dichroic mirror, WP - Waveplate. (a χ(3) process which utilizes the intensity-dependence of the refractive index) to generate a very small amount of light with wavelengths between 1300 nm and 1500 nm. Sapphire is chosen for this purpose for its mix of large nonlinear susceptibility and high-damage threshold; other materials (such as CaF2 ) have significantly higher susceptibilities, but are prone to damage. This weak continuum is overlapped, temporally and spatially, with 800 nm light in a β-Barium Borate (BBO) crystal to generate ≈3-5 µJ of signal. β-BBO is one of the most commonly used birefringent crystals used in nonlinear optics, particularly in the visible and near-IR frequency range. This is because of its high def f for a number of processes and high damage threshold. The final stage is the power-amplification stage, where the µJ signal seed is again overlapped with 800 nm light in a β-BBO crystal to produce > 200 µJ of signal and idler. The center frequencies of the output fields may be tuned by varying the angle of the β-BBO

80


crystals. The amplification stage is one area where experimental designs differ; two competing OPA designs use either the same β-BBO crystal for the pre-amplification and amplification stages in a double-pass geometry, or use two different β-BBO crystals for these stages. The double-pass geometry requires fewer optics and crystals; however, alignment in this geometry is significantly more challenging. Using multiple crystals requires more careful design, but offers more degrees of freedom for alignment. The OPA here utilizes TypeII phase-matching, where the signal and idler have orthogonal polarizations (in contrast to Type-I, where signal and idler have the same polarization, orthogonal to the pump polarization). A schematic of the layout of the commercial OPA used in these experiments (TOPAS, Light Conversion Inc.) is shown in Figure 3.2

3.2.4

Difference Frequency Generation

The signal and idler from the OPA are lowered in height with a periscope, which has the effect of rotating the polarization such that the signal is p polarized and the idler is s polarized. The signal and idler are aligned into a modified Michaelson interferometer where a dichroic mirror is used to separate the signal and idler into different arms of the interferometer. The interferometer is misaligned such that the reflected signal and idler are still collinear, but exit the interferometer displaced from the input beam. This allows the time delay between the signal and idler fields to be controlled; the time delay between these fields is key to the stability of the output of the DFG process [80], and different group velocities can result in a time delay between the signal and idler fields. The layout of the DFG setup is shown in Figure 3.3 The re-timed collinear signal and idler fields are combined in a AgGaS2 crystal. Using this crystal, 5-15 µJ (depending on wavelength) of mid-IR light can be generated. One aspect of DFG which is rarely explicitly discussed is that the process generates light at the input idler wavelength (here, ≈ 2 µm) and mid-IR light. Though from an energy conservation perspective, this may seem trivial, it is nonetheless important to recognize the


81

Figure 3.3: Schematic layout for the DFG apparatus and visible tracer overlap. AGS AgGaS2 crystal, DCM - Dichroic mirror, Ge - AR-coated Germanium plate. presence of additional idler light because of the experimental challenges of separating the near-IR and mid-IR frequencies. The leftover near-IR light is filtered from the mid-IR by an anti-reflection coated Germanium plate. Ge transmits a significant amount of leftover idler light; however, if the Ge plate is not aligned normal to the incident beam, the idler and mid-IR will be spatially displaced. Ge has positive Group-Velocity Dispersion (GVD) in the mid-IR, in contrast to the other optical materials used in the experiment. Thus, Ge is able to pre-compensate for some material GVD to maintain a short pulse. The Ge plate also serves as a convenient way to overlap a visible CW laser (HeNe, 633 nm wavelength) with our mid-IR light. This visible tracer is essential for effective alignment of the interferometer. The HeNe tracer is reflected from the Ge plate such that it propagates colinearly with the mid-IR pulses. The interferometer is first aligned using the visible tracer, then the mid-IR beam is walked to the same position as the visible tracer by optimizing overlap in the near-field (≈ 30 cm from the source) and far-field (> 9.5 m from the source). Before entering the interferometer, the beam passes through a reflective telescope to expand and collimate the beam. Reflective optics are used to minimize the effect of dispersion on the pulses. The first optic in the telescope is a plano-convex mirror, followed

82


by a plano-concave mirror. This ensures the beam never comes to a focus; from practical experience, using a focusing geometry results in significant pulse distortion, presumably from self-phase modulation (and other nonlinear effects) in air, involving the mid-IR light, and any leftover near-IR light. The entire DFG and HeNe overlap optics are enclosed in a box held under slight positive pressure of dry air. All optics after the DFG stage are also kept in a dry air environment. This enclosure serves two purposes. First, the absorptions from water vapor can have deleterious effects on the mid-IR pulse. The reduction in intensity from absorption is actually one of the least disruptive aspects of water vapor. The refractive index around these absorptions changes significantly with wavelength, giving rise to significant temporal dispersion in the pulse. Finally, small variations in the density of water vapor across the cross-section of the beam can degrade the spatial coherence of the pulse. This will result in much looser focusing at the sample, and a more complicated temporal profile of the signal. The second role of the enclosures is to isolate the optical path from air currents in the room. This has a significant effect on stability. The entire optical path, from regenerative amplifier to signal detection, is enclosed. Care must be taken with the boxes purged with dry air to ensure that the purging itself does not introduce turbulent air currents into the beam path. For this reason, the airflow is kept at a very slight positive pressure, and regulators and diffusers are attached to the purge lines to limit the flow rate and minimize turbulent air currents.

3.3

Boxcar Interferometer

At the entrance to the interferometer, a periscope is used to rotate the p polarized midIR light to s polarization. The beam is then split into five different beams, as shown in Figure 3.4. ZnSe beamsplitters are used to ensure that ≈99% of the intensity is found in the three beams used to generate the nonlinear signal, labelled α, β and γ. (These beams are sometimes referred to as k1 , k2 and k3 , respectively. This notation is avoided here to

3.3 Boxcar Interferometer

83

prevent confusion with the wavevectors of the time-ordered pulses used in the description of phase-matching in 2D IR). The remainder of the intensity is split between a beam to be used as the LO and a beam which traces the signal beam path (only used for alignment purposes). ZnSe is used for the beamsplitter material because of its low dispersion for light 6 µm wavelength. Despite this, ZnSe compensation plates are added to each arm of the interferometer to ensure each beam passes through the same amount of material. Polarization is controlled independently for α, β and γ and the LO by a waveplatepolarizer (Molectron, > 100 : 1 extinction ratio) combination in each arm of the interferometer. Computer-controlled motorized rotation mounts (PRM1Z8, Thorlabs) are used in the α and β beams, which enables switching between different polarization conditions to be done automatically. The intensity of LO is controlled by rotating the waveplate while keeping the polarizer fixed. Only the projection of the rotated polarization vector is transmitted. Note, however, that waveplates are generally designed for a specific wavelength; unless there is a fortuitous variation in the refractive indices of the fast and slow axes with wavelength, a half-waveplate will only produce linearly polarized light at one wavelength other elements of the spectrum will be elliptically polarized. Attempting to attenuate the LO intensity too much with this waveplate-polarizer combination can lead to significant spectral distortions in the LO. The time delays between beams are controlled by the use of motorized stages (ANT-90, Aerotech), on which corner cube retroreflector mirrors are mounted. Beams α, β and LO are controlled by motorized stages, while γ is used as a reference point, so it is mounted on a manual translation stage instead. The time delay between pulses is determined by scattering pairs of beams through a 50 µm diameter pinhole. The maximum of the envelope of the interferogram shows the point where the delay between the two pulses is zero. The corner cube retroreflectors also serve as a convenient way to vertically displace beams α and γ from β and the tracer. These beams are focused into the sample by an off-axis parabolic mirror. The approximate size of the focused beam can be obtained by determining the intensity throughput through several different sized pinholes (here, 50 µm,

84


Figure 3.4: Schematic layout for the boxcar interferometer and sample area, including the optional coupling of an external nanosecond laser for relaxation experiments. BS - ZnSe beamsplitter, CP - ZnSe compensation plate, Pol. - Polarizer, WP - Half-waveplate 75 µm and 100 µm diameter pinholes) and using the relation for the transmission T of a Gaussian beam with waist w through a circular aperture of radius r: T = 1−exp (−r2 /2w2 ). The diameter of the beam (twice the waist) at the focus is ≈75 µm. Beams α, β, γ and tracer are incident on the parabolic mirror in the boxcar configuration. This geometry has the four beams propagate parallel to each other as if each beam defined the corner of a “box”. When focused to a common point by the parabolic mirror, each of these beams has a unique wavevector; all third order signals originating from some combination of α, β and γ will be generated in a unique, background-free direction. The time ordering of α, β and γ determines which signal will be observed in the same direction as the tracer beam (i.e which signal will be aligned into the detector). The Rephasing spectrum (see later) is detected when the pulses are ordered α − β − γ, while the NonRephasing spectrum is observed for the ordering β − α − γ. In this way, both components of the 2D IR correlation surface may be detected simply by switching the ordering of α and β. The beams focus to a waist of ≈40µm, determined by the throughput of the incident

3.3 Boxcar Interferometer

85

beams through a series of pinholes of various sizes. The LO is also incident on the parabolic mirror, and is focused into the sample region in the same manner as the other beams. The LO is incident at the sample ≈ 100 ps before the other beams. This ensures the LO pulse has the same spectral characteristics as the signal, while also ensuring the LO plays no role in signal generation. The detection scheme outlined in 3.4 requires that the LO and signal field be timed up at the detector. This requirement that the LO be temporally-separated from the other beams at the sample, but arrive at the same time as the signal at the detector precludes the LO from taking a path parallel to α, β and γ. To achieve this non-parallel path, the LO beam enters the parabolic mirror in a non-parallel fashion, exaggerated in Figure 3.4. This results in a mild astigmatism of the LO beam. The primary disadvantage of this is in the collection of pump-probe spectra (used for phasing the data, see §4), where β and LO are used as the pump and probe, respectively. A chopper is placed in the β beam, which chops the beam at 500 Hz. On-thefly subtraction of alternating shots yields only signals that include an interaction with β, and are detected in the correct phasematched direction.

One particular con-

cern for studying membranes and membrane protein samples is light scattered from the intense beams α, β and γ reaching the detector. Scattered light from any of these beams individually is not a concern Figure 3.5: Balanced detection scheme. The (unless it is so large as to saturate the designal and LO are combined on a beamsplitter AR coated on one side. tector) since we are primarily interested in signals which oscillate as a function of the delay time between α and β. Scattered light from α, β and γ which interferes with other sources of scatter is a greater concern. By chopping,

86


we eliminate our sensitivity to scatter between α and γ. Scatter between β and γ is also not a concern; the delay between these pulses is fixed during a measurement, so this contribution is static. Scatter between α and β, however, is of greater concern. This scattered light appears as an intense series of peaks along the diagonal of the 2D IR spectrum. There are several strategies for minimizing scattered light between α and β. The most commonly employed of these is to use a ZZYY polarization scheme - α and β are polarized perpendicularly to γ and LO. By placing a polarizer before the monochromator, α and β are attenuated by a factor of > 100. Lipid membrane samples, however, are able to rotate the polarization of incident fields; to minimize scatter in lipid samples, careful sample preparation is a more effective route. This is discussed more in later chapters for specific types of sample.

3.4

Signal Detection

The third-order signal generated from the sample is collimated with an off-axis parabolic mirror, overlapped with the LO beam and dispersed by a grating monochromator onto a HgCdTe (MCT) array. The effect of the grating is to perform a cosine transform of the signal field by encoding the frequency information of the pulse onto a spatial dimension. In contrast to other 2D IR experimental geometries, the boxcar geometry requires that the signal field be overlapped with an additional field for heterodyne detection [81]. Though this adds an additional stage of alignment, it also allows a balanced detection scheme to be utilized. Here, the signal and LO fields are overlapped on a specially coated beamsplitter so that the signal field is either reflected from the front face of the beam splitter, or transmitted, and the LO is either reflected from the back face or transmitted. The phase shift of a beam on reflection depends on whether the beam is reflecting from a medium of higher or lower index. When light propagating through medium 1 (with refractive index, n1 ) reflects off an interface with medium 2 (refractive index n2 ) the beam will experience a phase shift of π if n2 > n1 , or 0 if n2 < n1 .

3.5 Temperature-jump Experiment

87

The reflection of the signal field off the front face of the ZnSe beamsplitter (medium 1, air, n1 ≈ 1, medium 2, ZnSe, n2 = 2.4) will experience a π phase shift. The reflection of the LO off the back face of the ZnSe beamsplitter (medium 1, ZnSe, n1 = 2.4, medium 2, air, n2 ≈ 1) will experience no phase shift. This is shown schematically in Figure 3.5. Thus, each arm of the balanced detection will contain a different combination of signal and LO fields. These two heterodyned fields can be displaced vertically and detected on the dual-stripe MCT array. The detected signal will be:

2 I1 = |ELO + Es |2 = Es2 + ELO + Es ELO

(3.17)

2 I2 = |ELO − Es |2 = Es2 + ELO − Es ELO

(3.18)

Subtracting these signals on a shot-to-shot basis yields 2Es ELO , and has the advantage of significantly suppressing correlated noise. This detection scheme has the additional benefit of avoiding the (sometimes questionable) assumption that the homodyne signal Es2 is small enough to be neglected; strong dipoles such as metal carbonyls, or even carboxylic acid stretches can give rise to large homodyne signals.

3.5

Temperature-jump Experiment

One extension of the 2D IR experiment described previously is to perturb the system somehow, and then use the 2D IR pulse sequence to probe the response of the system to this perturbation. One of the key challenges in these experiments is the time taken to perturb the system; the perturbation should be rapid relative to the timescales of interest in the system. The method of introducing this perturbation should also be specific to the perturbation type of interest. One of the earliest methods for introducing a rapid temperature-jump (T-jump) in to the system was to use the discharge of a capacitor in solution [82]; this produced a rapid increase in temperature of the system, but also introduced significant transient

88


electric fields. This made the interpretation of results where the system may respond to either perturbation (such as lipid vesicles) challenging. A sub-10ns T-jump can be achieved by exciting the overtone of the OD stretch of the D2 O solvent. Vibrational relaxation in the isotopomers of water is rapid, so this vibrational excitation relaxes to low frequency modes within the pulse duration. The stretch overtone is a forbidden transition in a harmonic system, so the absorption cross-section for this transition is very small. This ensures there is no significant depletion of the T-jump pulse across the sample, which creates a uniform temperature increase in the sample. This small cross-section also necessitates high-energy pulses; typically, less than 10% of the pulse energy is absorbed by the sample. To produce high-energy nanosecond pulses with a wavelength of 1.9µm the second harmonic of a Nd:YAG laser is used to pump an Optical Parametric Oscillator (OPO), which produces idler pulses at the required wavelength. Pulse energies of approximately 10 mJ are required at the sample to produce a T-jump of 10-20 K. To achieve the required pulse energies, a flashlamp-pumped Nd:YAG with oscillator and amplifier stages (Quantel) is used. This provides pulses of 500 mJ at a center wavelength of 1064 nm. Flashlamp pumped Nd:YAG lasers typically have poorer transverse beam modes and are less stable than diodepumped solid state (DPSS) alternatives, but are able to produce much higher output energies. The output from the laser is doubled by a temperature-controlled KDP crystal. Residual 1064 nm light is filtered out by use of dichroic mirrors. The 532 nm light is used to pump an OPO. An OPO operates in a similar manner to an OPA; both utilize a seeded OPG process to produce signal and idler fields. The main difference with an OPO is that there is no external seed; the nonlinear crystal is placed in a cavity with an output coupler which reflects a small amount of signal field back to the crystal. In this way, the initial signal and idler are generated by OPG, but this signal then acts as a seed for subsequent processes. A β-BBO-based OPO (Opotek) is used to generate an idler field at 1.9 µm. This near-IR light is focused into the sample to a beam-waist of ≈300 µm. This is significantly larger than the focus of the mid-IR pulses to ensure these pulses probe a

3.6 pH-jump and Other Perturbations

89

uniform region of the sample. The T-jump laser operates with a repetition rate of 20 Hz which is usually sufficient for the sample to have returned to equilibrium between subsequent shots. In the case where millisecond (and longer) time constants are present in the system, repetition rate of the T-jump laser can be lowered, at the expense of longer acquisition times. The delay between the T-jump laser and the mid-IR pulses is controlled electronically using the oscillator in the regenerative amplifier as a master clock for both systems.

3.6

pH-jump and Other Perturbations

Temperature is a widely used perturbation in so-called “relaxation” experiments. Two of the key advantages temperature has is the variety of methods to introduce a temperature change in the system (optical excitation, capacitive discharge) and the broad applicability of the perturbation. Consider the Van ’t Hoff equation: d ln Keq /dT = ∆H/RT 2 . Few equilibria have ∆H = 0. Changing the temperature is therefore a convenient way to alter the equilibrium of a system. For biophysical studies, however, rapid temperature changes represent a very unnatural perturbation. Useful information can be extracted from these experiments, but relating the observed dynamics to biological function becomes more challenging. To study the functional dynamics of a biomolecule, it would be desirable to perturb the element of the system (pH, ion concentration, membrane potential) that the biomolecule is directly sensitive to. The capacity to transiently increase the concentration of protons in the system would enable many interesting pH-gated processes to be studied. As one example, the conductivity of the bacterial potassium channel KcsA is regulated by the local pH [83]. Crystal structures are known for both the open and closed states of the protein [30], but the timescales, and the mechanism, of the transition between these states remains mysterious. Creating a transient population of protons in the system is fundamentally different to a T-jump experiment; the latter involves coupling energy from a laser into the system, while the former requires that

90


Figure 3.6: Structures and lifetimes of the three photoacids mentioned in the text - HPTS (left), NO2nH (middle) and o-NBA (right). the composition of the system is somehow rapidly changed. This rapid composition change can be achieved by the use of a photoacid - a molecule which, upon optical excitation, will release a proton. There are two classes of photoacids, each with unique advantages and disadvantages. The photoacids referenced in this section are summarized in Figure 3.6.

Reversible Photoacids Reversible photoacids operate on a common principle - they are molecules which contain an ionizable group (such as a hydroxyl group) which has a high pKa in the electronic ground state, but which have an electronic transition which results in a significantly lower excitedstate pK∗a . One of the most commonly studied reversible photoacids is 8-hydroxypyrene 1,3,6-trisulfonate (HPTS) which has pKa = 8.1, pK∗a = 1.6. The excited state lifetime


91

of HPTS is only a few nanoseconds [84]. This short lifetime makes HPTS unsuitable as a photoacid for studying biophysical systems. A longer-lived reversible photoacid is 1-(2nitroethyl)-2-napthol (NO2nH) which relaxes to a metastable intermediate, and can give transient pH changes which persist for many milliseconds [85]. However, this photoacid is only sparingly soluble in aqueous solution.

Irreversible Photoacids Irreversible photoacids undergo some sort of irreversible chemical change on optical excitation. The most commonly used irreversible photoacid is ortho-nitrobenzaldehyde (o-NBA). On UV excitation, this undergoes an isomerization to ortho-nitrosobenzoic acid, which can ionize to release a proton [86]. Irreversible photoacid are not subject to the same lifetimelimited pH-jump as reversible photoacids, but require either single-shot sensitivity or large volumes of sample for averaging.

3.6.1

Experimental Implementation

Both NO2nH and o-NBA have electronic transitions which lie in the UV wavelength range. The Nd:YAG laser used for the T-jump experiments can be modified to generate this excitation pulse. The doubled output of the T-jump laser (λ = 532 nm) can be doubled again in a β−BBO crystal (θ = 47.7◦ , 6 mm thickness, CASTECH Inc.) to produce light with a wavelength of 266 nm. This doubling process has an efficiency of ≈ 25% - for 50 mJ input energy, 12.5 mJ UV light is generated. Coarse optimization of the crystal angle with respect to the input beam can be achieved by eye. If a card covered with yellow highlighter ink is viewed through laser goggles with high OD at λ = 532 nm, UV-induced fluorescence on the card provides a convienient method for visualizing the generated UV. Overlap between the UV and IR pulses (both temporal and spatial) is achieved with a thin Si wafer. Si is transparent to mid-IR light; however, on excitation by the UV pulse, excited state intraband transitions reduce the transmission of the mid-IR light. Careful attention must be paid to the fluence of the UV beam; burned spots appeared rapidly on

92


the silicon wafer at fluences > 10 mJ/cm2 .

Limitations Both the reversible photoacid NO2nH and the irreversible photoacid o-NBA were tested as candidates for a pH-jump experiment. Each photoacid has significant obstacles which still remain to be overcome before they can be used with IR spectrscopy as a tool for studying biophysical systems. NO2nH is limited primarily by the solubility of the photoacid. Typical concentrations for transient 2D IR experiments lie in the 10 mM range (where concentration refers to concentration of oscillators, not molecules). Even after including a small amount of methanol in the sample as a solubilizing agent, the solubility of NO2nH is only 2 mM [85]. The quantum yield of the proton-releasing channel of NO2nH is < 0.5 - to see any significant changes in the 2D IR spectrum, UV fluences of ≈ 0.3 J/cm2 are required. (This estimation assumes a molar absorption coefficient ε266 = 3.5 × 103 M−1 cm−1 and that the system is initially near the pKa of the target, so all released protons protonate these groups.) At these fluences, damage to CaF2 (via UV absorption by impurities) and two-photon ionization of water to produce solvated electrons are significant problems. Solvated electrons give rise to a large response in the mid-IR. A related problem is the dissipation of the UV excitation - the fluences required for NO2nH would result in an increase in sample temperature of > 10 K. Disentangling the effect of this temperature increase from the effect of the released protons would present a significant barrier to interpretation of any observed signal. An additional problem is the UV absorption of the biomolecules themselves - both proteins and DNA have significant UV absorptions. Excitation of these UV transitions is often deleterious to proper function of the biomolecule. o-NBA is significantly more soluble than NO2nH, saturating at 8 mM with no solubilizing agent. Indeed, several studies exist probing the effect of rapid proton release from o-NBA with mid-IR pulses, even progressing as far as elucidating the mechanism of pHinduced folding in the leucine zipper [87]. Two significant problems remain to be overcome.


93

First, the irreversible nature of the photoacid requires a significant volume of sample. Even an impossibly-optimistic calculation where only a 3 nL volume (the irradiated volume) is refreshed at a 20 Hz rate (the laser repitition rate) would still require a 2 mL volume of sample for a single experiment (≈ 8 hour run time). The practicalities of building a flow cell, and ensuring fresh sample for each shot would increase this sample volume by at least an order of magnitude. For reference, the entire Lysozyme content of an egg would yield only ≈ 10 mL of sample under the concentrations necessary for IR spectroscopy. The final problem, common to all photoacids, is that even both the number of protons released and the thermal heating of the sample will be proportional to the number of protons released, even for irreversible photoacids. For more complex biophysical systems, the protein itself will act as a buffer (sidechains not involved in the functional change may still be protonated or deprotonated). Effecting a detectable change in the protein conformation while avoiding thermal artifacts may not be possible with current photoacids for very many proteins.

Application to Difference Spectroscopy Though the use of photoacids for time-resolved studies of biophysical systems still has several experimental issues to overcome, irreversible photoacids are useful tools for acquiring highquality difference spectra. Preparing a sample containing o-NBA and collecting a spectrum (either FTIR or 2D IR) before subjecting the sample to mild irradiation with UV light (thus avoiding the thermal artifacts and two-photon ionization of water) and acquiring a new spectrum enables high-quality pH-difference spectra to be acquired. The use of photoacids avoids issues with preparing separate samples, such as changes in path length or sample concentration. Figure 3.7 shows that illumination of a photoacid can trigger the color change in a pH-indicator. Figure 3.7 also shows the pH-difference 2D IR spectrum of diglycine. An extension of this idea would be to continuously collect spectra of the sample under low, uniform illumination. As more and more of the photoacid dissociated, the pH of the

94


Figure 3.7: Left, the pH changing capability of o-NBA is demonstrated using the pH indicator bromophenol blue. On exposure to a low intensity UV source (λ = 366 nm), the indicator reports a significant pH change. Right, the pH-difference 2D IR spectrum of diglycine using a photoacid. Protonation of N-terminus results in a shift in the carbonyl peak frequency. The contours of 2D IR spectrum of the protonated N-terminus are plotted on top of the difference spectrum. system will lower; this may provide a convenient, high-quality alternative to the laborious task of collecting pH-series data manually.

3.7

Pulse Characterization and Nonresonant Effects

It is clear from the preceding sections that the characteristics of the input fields play a significant role in nonlinear optics. Determining the spectrum of the pulse is relatively straightforward; either a grating and a camera, or an interferometer and a single-channel detector can yield the pulse spectrum. Determining the temporal profile of the field, however, is significantly more challenging. Fundamentally, this difficulty arises from the timescales involved. In order to measure femtosecond pulses, any detector would have to have a response time also on the femtosecond (or shorter) timescale. Furthermore, direct measurement of the detector response would require electronics with many tens of THz bandwidth. Both of these requirements are

3.7 Pulse Characterization and Nonresonant Effects

95

Figure 3.8: Transient-Grating FROG measurement with the PG-FROG time-ordering. The pulse shows minimal linear chirp, and yields a value of 97fs for the FWHM of the electric field. Over the region of interest (≈1500-1750 cm−1 ), the maximal difference in arrival time of different frequencies is 4 fs. beyond current technology. An alternative to direct electronic measurement of the pulse is to use an optical gate to measure the pulse length. The idea is to use some nonlinear optical process to take “snapshots” of the electric field of interest. If the resulting pulse is frequency resolved, the measurement is known as Frequency-Resolved Optical Gating (FROG) [88]. To consider the type of information this measurement can provide, consider the time-domain expression for the electric field:

E(t, z) = Eei(ω0 t+φ(t)−kz) .

(3.19)

From this, we can define the instantaneous center frequency of the pulse ω(t) = ω0 +

dφ dt .

A quadratic temporal phase gives rise to a time dependent instantaneous frequency; the leading edge of the pulse has a different frequency to the trailing edge. This is what is usually referred to as “Chirp”. Higher order terms of the temporal phase can give rise to more complicated effects. The effects of chirp in nonlinear spectroscopy can be significant. One of the most concerning, and often ignored, effects is a complicated frequency-dependent phase-shift

96


in the 2D IR spectrum. For a particular mechanical delay between pulses, the time delay between the same frequency in each pulse (e.g. red-red time delay) will be different from the time delay between different frequencies in each pulse (e.g. red-blue time delay). Thus, for a given mechanical delay, different pathways will contribute to the spectrum will have different time delays, which will introduce significant lineshape distortions into the spectrum. Clearly, before any detailed interpretation of a nonlinear spectrum can take place, the temporal characteristics of the pulses must be determined. FROG measurements offer a robust (and often convenient) way to characterize the spectral and temporal characteristics of a pulse. The detected signal of a FROG measurement is:

Z

∞

IS (ω, τ ) = |

g(t)E(t − τ )eiωt dt|2

(3.20)

−∞

where g(t) is the gating function, which depends on the type of the choice of nonlinear interaction. There are many variations of FROG which vary in the nature of the gate - these range from Second-Harmonic Generation (SHG-FROG), which uses SHG in a birefringent crystal as the gate, to Polarization Gating (PG-FROG), which uses the Kerr rotation of a probe pulse as a gate. Here, Transient-Grating FROG is used to determine the pulse characteristics. This type of measurement uses the non-resonant response of a polarizable material (here, ZnSe) to act as the gate; the signal is only generated when the pulses are overlapped. The main disadvantage of this measurement type is the requirement that the signal be emitted in a unique direction, necessitating a more complicated experimental geometry (such as a boxcar interferometer). However, since this is already implemented in the interferometer geometry for the 2D IR experiment, the TG-FROG measurement is straightforward. Some advantages of TG-FROG measurements are the unambiguous determination of chirp direction (because TG-FROG utilizes a χ(3) process rather than χ(2) for SHG-FROG, etc.) and the background-free signal detection. Experimental TG-FROG traces are shown in Figure 3.8 One subtlety in a TG-FROG experiment is which time delay to scan; this choice will

3.7 Pulse Characterization and Nonresonant Effects

97

alter the signal observed. If α is scanned, the observed signal becomes:

Z IS (ω, τ ) ∝

∞

∗

2

E (t − τ )E (t)e

2 dt

(3.21)

2 dt .

(3.22)

iωt

−∞

while if β is scanned, the signal becomes

Z IS (ω, τ ) ∝

∞ 2 iωt

E(t − τ )|E(t)| e

−∞

These signals are equivalent to either self-diffraction (SD-FROG) or PG-FROG measurements respectively. An experimental TG-FROG trace is shown in Figure 3.8, taken in the PG-FROG time-ordering. This trace shows minimal linear chirp, and only a very small amount of higher order chirp (evident in the curvature of the arrival time of each frequency). Analysis of this trace yields a pulse envelope of FWHM 97 fs, which compares favorably with the 90 fs pulses output from the regenerative amplifier. The reality of TG-FROG measurements can be more complicated than the discussion presented here. In particular, it is important to be confident that no higher order processes are occurring; these will complicate the interpretation of the FROG trace. ZnSe, in particular, will display unusual behaviors at high incident powers, such as an intense signal which grows in on a ps timescale and persists for ns, and visible (blue) emission from the ZnSe plate. Though non-resonant response is useful for characterizing the pulse, it can also be a hindrance in experiment. At short delay times between pulses, it can sometimes overwhelm small signals. For this reason, where other properties are comparable (e.g. dispersion), it is often desirable to use materials with low polarizability as sample windows (e.g. CaF2 ). Even using extremely short pulses does not fully alleviate the issue of non-resonant response; shorter pulses give rise to a more intense non-resonant response. Since the metric of interest is not the FWHM of the non-resonant response, but rather the time at which the nonresonant response is negligible compared to the signal, shorter pulses can actually result in more signal contamination from non-resonant signals. Non-resonant response is

98


not the only factor to consider when setting the waiting time in a 2D IR experiment, however. One of the primary assumptions in §2 was strict time-ordering of the interactions with the pulses. If the waiting time delay is significantly shorter than the convolution time of the pulses (≈ 1.44×FWHM) then significant signal contamination may occur - for example, signals arising from double-quantum coherence pathways may now be emitted in the same direction as those arising from Rephasing or Non-Rephasing pathways. This signal contamination makes lineshape interpretation challenging.

Chapter 4

Processing and Analysis in Transient Spectroscopy 4.1

Introduction

Raw data collected in a 2D IR experiment does not yield itself readily to interpretation. Many stages of processing are required to transform the raw data into a form which may be analyzed. These processing challenges are distinct from the challenge of interpreting the data. In temperature jump (T-jump) experiments, the challenges associated with 2D IR spectroscopy are compounded by the non-equilibrium nature of the system being probed. Tools for processing and analyzing the data fall into three broad categories. The first of these categories covers processing steps which correct for experimental deficiencies. One example of this is correcting for position error from the translation stages in the system; in principle, a stage with sufficient accuracy and precision would negate the need for this step. However, such arbitrarily accurate and precise stages do not exist, so additional processing steps are needed to compensate for this. The second category consists of operations which do not change the information content of the data, but do alter its presentation. The key example of this in these experiments is the transformation of the signal from the time-domain (experimentally obtained) to the 99

100

Chapter 4. Processing and Analysis in Transient Spectroscopy

frequency-domain. The frequency-domain is a more intuitive way to analyze the signals, specifically how different vibrational frequencies are correlated. Associated processing (such as windowing of the signals) does not strictly conserve the information content of the original raw data, but is included in this category anyway. The final (and broadest) category covers analysis tools. These are manipulations of the processed data designed to bridge the often-formidable gap between data and interpretation. Such is the importance of these tools for 2D IR that they have received extensive review in the literature (such as robust extraction of CLS [70] etc.), and are not the focus of this chapter. Rather, the less-explored area of analysis tools for mid-IR relaxation experiments is discussed. The remainder of this chapter is structured as follows: relevant processing and analysis techniques for 2D IR experiments are discussed, then the underlying principles of relaxation experiments are outlined, before finally, analysis tools suited to relaxation experiments are presented.

4.2 4.2.1

Processing Tools for 2D IR Spectra Phasing

Uncertainty in the positions of the translation stages used to control the delay between the first two pulses in the experiment is one of the most common sources of error in multidimensional optical experiments. Current state-of-the-art motion controllers are precise to sub-micron length scales, but not accurate! Furthermore, an optical cycle at the frequencies studied in this thesis is on the order of 20 fs. A stage positioning error of just 750 nm would give rise to a timing error equivalent to a quarter of a cycle (5 fs). Clearly, there is such little tolerance with respect to position error that some sort of post-processing is desirable. The first issue to address is what deleterious effects on the spectrum this timing error may have. All the response functions outlined in §2 have terms which contain products with exp(iωa0 t1 ). The absorptive (or correlation) spectrum is the real part of the Fourier

4.2 Processing Tools for 2D IR Spectra

101

transform of these signals. However, if there is some timing error in t1 , such that t1 → t1 + δt, then this will appear as a frequency dependent phase shift in our response function exp(iωa0 t1 ) → exp(iωa0 t1 ) exp(iφ(ω)). Thus, the real and imaginary parts of the “true” signal will become mixed. Uncertainty in δt renders lineshape analysis of absorptive spectra all-but-impossible. Figure 4.1 shows how even a small error in stage position can introduce significant spectral distortions. Determination of δt may seem at first impossible in post-processing - the “true” spectrum and δt are unknown. However, by considering the relation between a pump-probe experiment and 2D IR experiments, some progress may be made. Both are third-order experiments, the principal difference between the two being frequency resolution in the pump axis. Thus, because the first two interactions with the system are from the same pump pulse, t1 (and also δt) are necessarily equal to zero. The pump probe spectrum can then act as a reference point; when δt has been corrected for then

S 2D (t1 = 0, ω3 ) = S P P (ω3 ).

(4.1)

This is often expressed in the frequency domain as Z

∞

S 2D (ω1 , ω3 )dω1 = S P P (ω3 )

(4.2)

−∞

and referred to as the projection-slice theorem [89]. It is easy to demonstrate these two expressions are equivalent; the frequency resolution of the signal is related to the maximum delay time sampled, so when t1 = 0, there is no frequency resolution at all. This is equivalent to convolving the frequency-domain response with an infinitely-broad Gaussian function, which yields the latter expression. In the boxcar geometry, acquisition of the pump-probe spectrum is often the limiting factor with regards to sample concentration. To correct for the timing error in the boxcar geometry, the signal in the range ±20 fs from the translation stage t1 = 0 is compared to the pump-probe spectrum to determine where the “true” t1 = 0 is. The data are collected

102


Figure 4.1: Timing error corrected (bottom left) and uncorrected (bottom right) 2D IR spectrum of diglycine. Pump-probe signal and Rephasing signal at t1 =0 (upper left). The uncorrected spectrum has a timing error of 6 fs, which corresponds to a position error of < 1 µm.


103

in 4 fs increments; to enable finer determination of δt, the data is interpolated to arbitrary resolution using a cubic spline interpolation. For this interpolation to be effective, the oscillations must be well sampled - undersampling of the free-induction decay introduces lineshape distortions into the data. In the pump-probe geometry [81], the spectrum is often instead integrated over ω1 and a spectral phase shift applied to maximize agreement between the pump-probe spectrum and the projected 2D IR spectrum. This route is not applicable in the boxcar geometry each of the surfaces (Rephasing and Non-rephasing) contains only some of the pathways present in the pump-probe signal.

4.2.2

Fourier Transforms

With the data corrected for position error from the translation stage, the next processing step is to transform the 2D IR spectrum from the mixed time-frequency domain to the frequency-frequency domain. Though the Fourier transform is conceptually straightforward in this case (transforming an oscillating signal in time to a peak in the frequency domain), the application of numerical methods (specifically the Fast-Fourier Transform, FFT, algorithm) must be done with care to avoid artifacts. The first point to note is that the FFT algorithm requires regularly spaced points. The spacing of these points (∆t) defines the upper bound of the transformed frequency axis, according to the Nyquist sampling rate fmax = 0.5/∆t. The raw experimental data is not uniformly sampled - this is another consequence of the precise, but not accurate, translation stages used. The data should be interpolated to uniform spacing; this can be combined with the interpolation in the phasing step for computational efficiency. Another requirement for correct implementation of the FFT algorithm is that the data go smoothly to zero. Otherwise, this appears to the algorithm as a discontinuity, which gives rise to “ripples” in the frequency domain. Even the tedious exercise of scanning the delay between pulses to times much longer than the dephasing time of the signal does not remedy this problem - there is always some amount of noise from the detector. The solution lies in

104


Window Name

Functional Form

Notes

Gaussian

exp(−t2 /τ 2 )

Long tail, but least spectral distortion

Super-Gaussian

exp(−tn /τ n ) πt 1 + cos τmax

Sharper transition to zero than Gaussian

Hann

Hamming

α − β sin

Cosine

cos

πt τmax

πt 2τmax

Ripples are low amplitude

α and β can be adjusted to minimize ripples

Derivative of window is non-zero at edge

Table 4.1: Some common window functions used in time-domain spectroscopy and signal processing. post-processing: multiplying the raw signal by some function which (in principle) preserves the information in the signal, while enforcing a smooth decay to zero. These are referred to as window functions. There are very many window functions used in time-domain spectroscopies (and more generally, in signal processing), a few of which are summarized in Table 4.1. Each window will have a unique effect on the spectrum; the resulting spectrum will be a convolution of the Fourier transform of the raw spectrum and of the Fourier transform of the window. The most appropriate window function will vary on a case-to-case basis. The Gaussian window gives the smoothest spectrum, though the long tail of the Gaussian results in the most frequency-resolution loss. The Super-Gaussian function alleviates this long tail problem, though at the expense of some rippling. The Hann window is the most versatile; frequency resolution degradation is less than with a Gaussian window, and ripples are < 2% of the intensity of the primary peak. Note that one advantage of the time-domain implementation of 2D IR is the Fourierfilter incorporated into the processing. Otherwise deleterious signals (such as background scatter) which do not oscillate do not appear in the final processed spectrum; this is not the case for the narrowband-pump (holeburning) implementation of 2D IR [68].


105

The converse case may apply when collecting pump-probe measurements. These signals are expected to decay smoothly with delay time (or in some cases, show beating at the difference frequency between modes in the system). However, scattered light may interfere with this measurement; this scattered light will often result from interference between the pump and probe beams, and thus oscillate with the optical period of the light used. In this case, it may be advantageous to Fourier transform the time-domain data, apply a window in the frequency domain to select only zero- and low-frequency components, and then transform back to the time domain.

The question of resolution arises when considering the effect of the windowing functions. All windows necessarily broaden the spectrum; the only window function which would not result in loss of frequency resolution would be a delta function in the time domain. This would correspond to the trivial window function of multiplication by a constant. The frequency resolution of the raw data depends on the maximum time delay between excitation pulses measured, ∆ω1 = (3 × 10−5 τmax )−1 where ∆ω1 is in cm−1 and τmax is in fs. In the experiments presented later, a maximum delay of 3000 fs, corresponding to ∆ω1 = 11 cm−1 was used.

A related concept worth noting is “zero-padding”. In this procedure, elements of value zero are appended to the time-domain signal, effectively increasing the scan time. The output of the FFT algorithm for this zero-padded data appears to be at much higher frequency resolution than the raw data. However, care should be taken with such an approach, as the resolution of the data is not actually higher than the data without zeropadding. The effect of zero-padding in the time domain is equivalent to an interpolation in the frequency domain. These may be useful tools for presenting data, or comparing different datasets, but should not be confused with a true increase in frequency resolution.

106


4.3

Relaxation Experiments

Relaxation experiments, pioneered by Manfred Eigen in the 1950s, use a rapid perturbation to shift the equilibrium of the system under investigation, and then monitor the relaxation of the system to this new equilibrium state [90]. In order to interpret the relaxation rates measured, it is first necessary to develop a theoretical framework of how the system responds to a perturbation [82]. ka

− * To begin, consider the simple chemical equilibrium A − ) − − B. We can write a differential kb

equation which describes the time evolution of the concentration of these species, cα ,

c˙a = kb cb − ka ca

(4.3)

c˙b = ka ca − kb cb

(4.4)

The definition of equilibrium requires that c˙α = 0. Equation 4.3 can then be rearranged to relate the equilibrium constant of the system to the rates present in that system. However, relaxation experiments are non-equilibrium experiments, so c˙α 6= 0. The time-dependence of this concentration is what relaxation experiments seek to determine.

One of the main challenges in the discussions of relaxation experiments is notation. Here, c¯βα will be used to denote the equilibrium concentration (the bar signifies equilibrium) of species α at either the final (after perturbation, β = f ) or initial (before perturbation, β = i) conditions. Conservation of mass requires that the following hold true: ca (t) = c¯fa + x(t) cb (t) =

c¯fb

(4.5)

− x(t)

where x(t) is some time-dependent change in concentration. Combining Eqs. 4.3 and 4.5 yields

4.3 Relaxation Experiments

x˙ = kb (¯ cfb + x) − ka (¯ cfa − x)

107

(4.6)

= −(ka + kb )x

Where the equilibrium condition ka c¯a −kb c¯b = 0 has been used to simplify the expression. This equation then yields exponential solutions with a characteristic time constant τ −1 = ka + kb . Similar expressions can be obtained for more complex reaction schemes, such as kab

−− * A+B) − − C. kc

However, an implicit assumption has been made in this treatment that the system is instantaneously perturbed, i.e. c¯˙fα = 0. This may be a good approximation in some cases (such as studying the µs response of a system to a < 10ns perturbation), but for observed timescales of the same magnitude as the perturbation, this is not true. Additionally, it enforces a step-function form on the perturbation; other functional forms of the perturbation may be advantageous, as discussed later. To explicitly consider cases where the equilibrium concentrations are time-dependent, a time-independent reference state, c∗α is defined. Typically this is taken as the equilibrium before the perturbation is applied, ciα . With this, it is then possible to define the following: ∆cα (t) = cα (t) − c∗α

(4.7)

∆¯ cα (t) = c¯α (t) − c∗α ,

(4.8)

cα = c¯α (t) + ∆cα (t) − ∆¯ cα (t)

(4.9)

which can be rearranged to yield

Defining x(t) = ∆ca (t) = −∆cb (t) as before, introducing x ¯(t) = ∆¯ ca (t) = −∆¯ cb (t), and noting that the definition of equilibrium still holds (ka c¯a (t)−kb c¯b (t) = 0) yields the following

108


differential equation: 1 ¯) x˙ = − (x − x τ

(4.10)

This first-order differential equation has the solution −I

x(t) = x0 e

where I =

R

+e

−I

Z

1 x ¯(t)eI dt τ

(4.11)

τ −1 dt. If τ˙ 6= 0, interpreting the experimental data becomes a herculean task.

First, a functional form for the rates must be assumed. Even if a simple Arrhenius model is assumed, detailed knowledge of the potential energy surface of the system is required. Assuming the rates in the system are time-independent then Equation 4.11 becomes

− τt

x(t) = x0 e

1 t + e− τ τ

Z

t

x ¯(t)e τ dt.

(4.12)

The functional form of x ¯(t) will vary between experiments. A few example cases with tractable solutions are discussed below.

4.3.1

Instantaneous Perturbation

This is the “ideal” case for many experiments, where the perturbation is so rapid relative to the dynamics of interest that it may be approximated as a step function, x ¯(t < 0) = 0 (4.13) x ¯(t ≥ 0) = x∞ . If the initial state is chosen as the reference state (c∗α = ciα ), then x0 = 0 and t x(t) = x∞ 1 − exp (− ) . τ

(4.14)

This is the simplest transient experiment to interpret; the measured time constant reflects exactly that of the system.


4.3.2

109

Exponential Perturbation

Prior to the development of optically-induced T-jump methods, Joule heating (capacitive discharge) was a commonly-used method of rapidly increasing the temperature of the sample. This causes a T-jump of the form: x ¯(t < 0) = 0 x ¯(t ≥ 0) = x∞

(4.15)

t 1 − exp (− ) . τd

The solution to the full relaxation equation (using the initial state of the system as the reference state) is then x(t) = x∞ 1 +

t τ t τd exp (− ) − exp (− ) τd + τ τd τd + τ τ

(4.16)

It is straightforward to see that as τd → 0, the expression for an instantaneous perturbation is recovered. Equation 4.16 highlights an important detail ignored when considering only an instantaneous perturbation - the perturbation itself can influence the measured timescale. For values of τ on the order of τd , it becomes extremely difficult to accurately recover the two timescales, or indeed identify that multiple timescales are present in the system. Even in model data with only 1% additive white noise (which would be extremely high quality experimental data), it is not possible to distinguish the presence of timescales when τ ≤ 3τd . This is one of the most significant drawbacks of the Joule heating method for T-jump; the discharge typically occurs with a µs timescale, which limits this method to studying processes with timescales > 10µs. Other drawbacks include the necessity of high ionic strength conditions (required for the discharge to produce a significant increase in temperature) which may alter the equilibrium of the system (such as in the case of DNA dimerization), and the high electric fields produced during the discharge, which may act as a secondary perturbation (as with poration in lipid membranes). However, Joule heating is

110


significantly simpler (and less expensive) to implement than optical methods; Joule heating spectrometers are commercially available.

4.3.3

Optical Excitation with a Gaussian Pulse

The most relevant type of perturbation to the experiments presented in this thesis is optical excitation of the sample. Assuming a Gaussian form for the excitation pulse, the system perturbation is given by

Z x ¯(t) =

−

exp

x∞ t2 dt = 1 + erf(t/τ ) . τp2 2

(4.17)

This yields (after tedious integration) − τt

x(t) = x0 e

2 τp τ τp x∞ t t − τt 4τ 2 − τt 2τp + 1 + e e erf − , + erf −e e 2 2τ τp τp

(4.18)

where x∞ x0 = 2

exp

τp 2τ

− exp

τp2 4τ 2

.

(4.19)

As before, as τp → 0, the expression for an instantaneous perturbation is recovered. Using this expression it is possible to estimate the error in the recovered timescale as a function of the separation of τ and τp . The error in the recovered timescale is < 1% when τ ≥ 6τp . For the 10 ns FWHM pulses used in the experiments presented here, this puts a lower limit of 35 ns on timescales which can be reliably measured without explicitly considering the form of the perturbation. One of the principal advantages of optical T-jump methods is the rapid perturbation; in some cases perturbations on a sub-ns timescale are possible. The principle disadvantage of this technique are the experimental challenge of coupling the optical excitation into the sample; a strongly absorbing sample will result in a temperature gradient in the direction of propagation of the beam, while a weakly absorbing sample will require a very intense beam.


4.3.4

111

Oscillatory Perturbation

The full relaxation equation (Equation 4.11) suggests that instantaneous perturbations are not the only approach to relaxation experiments: if x ¯(t) can be controlled, this may provide route to accessing the dynamics. Consider an oscillatory perturbation, x ¯(t) = A cos(ω0 t). For times t >> τ , then the response of the system to this perturbation is

x(t) =

Aω0 τ A cos(ω0 t) + sin(ω0 t). 2 2 1 + ω0 τ 1 + ω02 τ 2

(4.20)

This form of perturbation is used extensively to measure the response of systems to pressure and electric field perturbations (in the form of ultrasound and dielectric relaxation experiments). Experiments of this type are limited by the highest frequency perturbation they can measure. However, there are some significant advantages to this approach to determining τ . Notably, it is straightforward to incorporate a Fourier-filter into the data analysis. Fourier transforming x(t) yields =(˜ x(ω)) = ω0 τ. 500 µs becomes challenging with this perturbation. However, mod-

4.4 Extracting Rates from Transient Spectra

113

eling the response of the system using the experimentally-determined temperature profile offers a route to estimating these timescales. Figure 4.3 shows the calculated transient response for timescales τ , where the pertur c bation is x ¯ = exp − (t/τr ) , τr = 2 ms c = 0.65. These simulations reveal that three parameters of the observed response are sensitive to the input timescale: the rise time of the response, the relaxation of the transient signal, and the functional form of the relaxation. For slower rates, there is a significant lag time between the relaxation of the temperature perturbation and the relaxation of the transient signal, and the form of the relaxation varies smoothly from the stretched exponential form of the perturbation to the pure exponential form of the system kinetics. The temporal temperature profile of the system may be obtained by monitoring the change in transmission of the solvent after the T-jump pulse. Using this measurement and the three constraints demonstrated by the simulation of the transient response, it may be possible to accurately estimate millisecond timescales in the system.

4.4

Extracting Rates from Transient Spectra

The simple reaction scheme used in Section 4.3 circumvents one of the major challenges of interpreting experimental data: how many timescales are actually present in the system? Fitting the data to more timescales might improve the quality of the fit, but do these timescales have any physical meaning? A natural place to begin when tackling this issue is a microscopic picture of the system. There are three primary sources of multiple timescales in the system: sample inhomogeneity, the presence of intermediates in a reaction scheme, and multiple reaction pathways - these are summarized in 4.4. Sample inhomogeneity is when a distribution of states is initially present, each of which will proceed along a common reaction coordinate with different rates (e.g. if a distribution of protein conformational states is initially present, the unfolding might not be described

114


Figure 4.3: Calculated transient response for a perturbation x ¯ = exp −(t/τr

)c

, τr = 2 ms,

c = 0.65, for a variety of input timescale τ . Upper left panel shows the extracted rise time of the response as a function of input timescale. Upper right shows the decay time as a function of input timescale. Bottom right panel shows the exponent of the stretched exponential as a function of input timescale. Bottom left panel shows calculated responses for a variety of input timescales. Dashed horizontal lines represent the timescale or stretch coefficient of the perturbation.


115

Figure 4.4: Three scenarios which can give to multiple observed timescales in relaxation experiment. by a single rate). −− * −− * More complex reaction schemes, such as A ) − − B ) − − C may give rise to multiple observed timescales. However, this will be highly dependent on whether the states are distinct with respect to the method of probing. For example, if states B and C have the same IR absorption signature, then an IR probe would only see one timescale for the system, −− corresponding to the A ) −* − B0 equilibrium, where B0 is some combination of B and C. Multiple reaction pathways (such as multiple dissociation mechanisms of a dimer) will also give rise to multiple timescales. All of these effects may be considered together by defining a probability distribution for rates in our system, P (λ). In keeping with Section 4.3, it is assumed that the rates are time-independent (and thus P (λ) is time-dependent). The observed relaxation signal can then be written as:

Z I(t) =

∞

P (λ)e−λt dλ.

(4.23)

0

This is simply the Laplace transform of P (λ); determining the distribution of rates in our system must be as straightforward as performing an inverse Laplace transform on the experimental data. Unfortunately the inverse Laplace transform (iLt) is an example of an ill-posed problem. In this case, the problem is ill-posed because small changes in the intial conditions (here, noise in experimental data) result in very large changes in the solution. However, there are strategies to solve ill-posed problems which may be of use here. These are known as regularization methods; by providing additional information on what the

116


solution should look like, it is possible to find robust solutions to iLt. Different regularization methods specify different constraints on the solution. Here, two different constraints are discussed - Tikhonov regularization (which biases the solution to have small norms) [92] and Maximum Entropy regularization (which biases the solution to minimize an informational entropy term) [93].

4.4.1

Tikhonov Regularization

Tikhonov regularization is one of the most common regularization methods. In spectroscopy, Tikhonov regularization has found widespread use in the reconstruction of distance distributions in Double-Electron-Electron Resonance (DEER) experiments [92]. The primary advantages of Tikhonov regularization are the computational efficiency of the procedure and the well-characterized behavior (discussed later) of this regularization method. Regularization methods, in general, aim to find a stable solution x to the problem Ax = b. In this case, x is a vector which describes the rate distribution, b is the experimental data and A is the Laplace transform. Tikhonov regularization returns the solution which minimizes the function ||Ax − b||2 + µ||x||2 . The solution which minimizes this function is

xT ik = (A> A + µ2 I)−1 A> b.

(4.24)

However, a more useful form of the solution may be calculated using Singular-Value Decomposition (SVD). If A = U ΣV > , then the solution may instead be expressed as

xT ik = V DU > b

(4.25)

where Dii = σi /(σi2 + µ2 ) and σi are the diagonal elements of the Σ matrix. This form is computationally efficient to calculate, and straightforward to implement. Additionally, comparison with the Pseudo-Inverse matrix provides some insight into what the regularization procedure actually does. The Pseudo-Inverse of A is defined as A+ = V Σ+ U > , where Σ+ is calculated by inverting all of the non-zero elements of Σ and transposing the matrix.


117

Comparing Equation 4.25 to the Pseudo-Inverse reveals that the role of the regularization is to re-weight components of the SVD before calculating a “regularized” Pseudo-Inverse. The only remaining element to determine prior to reconstruction is the value of the regularization parameter, µ. Here, the properties of Tikhonov regularization are advantageous; by systematically varying µ and calculating the residual, ||Ax − b||, and the norm, ||x||, it is possible to determine an optimum value for µ. Plotting ||xµ || against ||Axµ − b|| yields a curve with a distinct “L” shape [94]. Examples of this are shown in Figure 4.5. The steep part of this curve is dominated by noise contamination - the parameter µ is not sufficiently large to recover the regularized distribution. The flat part of the curve arises from regularization error - the parameter µ is so large that the recovered distribution no longer bears any resemblance to the experimental data. The corner of this curve provides a parameter µ which satisfies the Tikhonov condition (minimize ||x||) while reconstructing the experimental data. The shape of this curve also provides a robust criterion to determine µopt - the value of µ where the curvature of the L-curve is maximized. Figure 4.5 shows the Tikhonov regularization reconstruction for a noisy biexponential time trace. In the noisiest case, it is not clear by eye how many timescales are involved. The reconstructed rate distributions clearly highlight the biexponential nature of the data, and are able to recover the timescales. The width of a peak in the rate distribution is related to the noise in the data. This is an important point; regularization methods cannot arbitrarily compensate for noisy data - for very noisy data, the peaks in the rate distribution may be so broad as to provide no useful information. The rate distributions in Figure 4.5 highlight one of the disadvantages of using the Tikhonov regularization condition - low amplitude artifact peaks are present. In many cases, the “true” reconstructed peaks are large enough to be readily distinguished from these artifactual peaks. However, if a low amplitude component of the rate distribution is of interest, Tikhonov regularization may not be an appropriate method. An example script for implementing Tikhonov regularization is outlined in Appendix 4.A.1. One important point to note is that the regularization condition in this case (min{||x||})

118


Figure 4.5: Model biexponential decay, I = exp(−t/τ1 ) + exp(−t/τ2 ) where τ1 = 100 µs and τ2 = 1 µs, with 5% noise level (upper left, dark blue) and 15% noise level (upper right, light blue). Bottom left panel shows the L-curve used to determine the optimum value of the regularization parameter, µ. Bottom right panel shows the optimized rate distribution reconstructions, with the ideal rates denoted by black lines.


119

permits negative values of x. The sign of the rate distribution amplitude reflects the sign of the equivalent coefficient in a sum of exponentials e.g if the data followed the form I(t) = a(exp(−λ1 t) − exp(−λ2 t)), the rate distribution would show two peaks of equal magnitude (at λ1 and λ2 ) with opposite signs.

4.4.2

Maximum Entropy Regularization

An alternative regularization method is to use a Maximum Entropy constraint. Maximum Entropy (MaxEnt) regularization seeks to maximize the entropy of the distribution, defined P as S = − i xi ln(xi /Xi ) − 1 , where X is some model reference distribution (often referred to as a “Prior”) which represents the information already known about the system [95]. This is typically taken to be a flat, low amplitude distribution. MaxEnt regularization is less prone to the low amplitude artifacts found in Tikhonov rate distributions, however it is significantly more computationally taxing. To implement the MaxEnt reconstruction, the function F is first defined as

F = ηχ2 − S

(4.26)

where χ2 =

N 1 X (Iexpt (tk ) − Icalc (tk ))2 . N σk2

(4.27)

k

The term which characterizes the noise, σk , is typically assumed to be independent of time, and can thus be taken outside the sum. The “true” solution should have a residual which is only the noise in the data, thus χ2 = 1 for the true solution. It may appear that MaxEnt has two free parameters, η and σ, but this is not infact the case. Rather, η is a Lagrangian multiplier which ensures that χ2 = 1. The only unknown term to determine in the MaxEnt expression is the experimental noise σ. This may be estimated from transient data by probing time points before the perturbation; there should be no transient response at this point. The response observed here gives a reasonable estimate of the noise level, though it often underestimates the value of σ.

120

Chapter 4. Processing and Analysis in Transient Spectroscopy Finding the distribution x which minimizes F is not straightforward; empirically, for

T-jump data, the most robust method is to use a Newton-Raphson minimization scheme. One procedure for finding xM E is: 1. Set η to an initial value of ≈ σ, and assume a uniform, small magnitude Prior, X. Typical values for X are 1 × 10−5 to 1 × 10−4 . 2. Minimize F with the Newton-Raphson algorithm, using X as the initial guess. 3. Evaluate χ2 . If χ2 ≤ 1, terminate algorithm. If χ2 > 1, continue. 4. Increase η by a predefined increment. 5. Minimize F with the Newton-Raphson algorithm, using the previously optimized distribution as the initial guess. 6. Evaluate χ2 . If χ2 ≤ 1, terminate algorithm. If χ2 > 1, return to Step 4 and repeat. An example script to implement this algorithm is given in Appendix 4.A.2. The form of the entropy term precludes negative coefficients in the rate distribution. To address this, two entropy terms (each corresponding to to the positive and negative distributions) are calculated.

4.4.3

Application of Regularization Methods

To demonstrate the strengths and weaknesses of the regularization methods discussed above, they are applied to real experimental data. T-jump experiments can be used to induce the dissociation of DNA duplexes; the dissociation of DNA duplexes may proceed with multiple timescales. The self-complimentary DNA sequence 5’-ATATGCATAT-3’ shows a frequencydependence to the transient response [96]. Figure 4.6 compares the rate distributions obtained from Tikhonov and MaxEnt regularization with fits to sums of exponentials. Both the MaxEnt and Tikhonov regularization methods are able to recover two time constants for the temporal response of the A


121

Figure 4.6: Upper left, transient-HDVE data for the 5’-ATATGCATAT-3’ DNA oligomer. Upper right, temporal profiles of vibrations corresponding to the G and A ring modes. Bottom row shows the reconstructed rate distributions using either the MaxEnt (left) or Tikhonov (right) regularization scheme. Solid black line indicates rates obtained by fitting the data directly. Dashed black line separates timescales associated with thermal relaxation of the sample.

122


nucleobase, and only one timescale for the G nucleobase response, consistent with the results obtained from fitting the data directly. The rate distributions highlight some of the advantages and disadvantages of each regularization method. Tikhonov regularization is computationally efficient, allowing many points in the rate distribution can be determined. However, it is subject to “ripple” distortions in the distribution, which may be difficult in some cases to distinguish from true peaks. The MaxEnt method is less subject to these artifacts, and thus provides a more reliable distribution. However, the MaxEnt method is much more computationally taxing, and so is restricted to much coarser sampling of the rate space. Of course, these methods are not mutually exclusive. A robust analysis of the data should utilize both regularization methods; if the MaxEnt and Tikhonov regularization methods reach a consensus on the number of timescales, and the shape of the distributions, present in a system, the data may be then be fit to this functional form to determine the timescales to a precision beyond either regularization method.

4.5

Summary and Outlook

Transient experiments (either on a ps timescale, like 2D IR, or ns-ms T-jump measurements) require significant processing and analysis to extract useful information. Sub-ps experiments must be corrected for the limitations of the equipment - errors in stage motion lead to lineshape distortions, and noise on the detector requires processing to be compatible with numerical FFT methods. Transient experiments on a longer timescale pose different challenges, such as how to map measured timescales to the true timescales in the system. The proper method of analysis of these longer-timescale transient experiments depends on the form of the perturbation. This opens new possibilities for studying the dynamics of biomolecular systems. For example, fast timescales may be best captured by a sudden perturbation to the system, while slower dynamics may be better observed with an oscillating perturbation. Using multiple

4.5 Summary and Outlook

123

implementations of transient experiments has the potential to greatly extend the range of timescales accessible.

124

4.A 4.A.1 1 2 3 4 5 6 7

% % % % % % %


Example Matlab Scripts Tikhonov rates.m

Tikhonov_rates.m - Paul Stevenson 2017 An example script for the implementation of Tikhonov regularization of the inverse Laplace transform of input data This script demonstrates how SVD can be used to efficiently calculate the regularized rate distribution

8 9

%%

10 11 12 13

close all; clear all; clc;

14 15

%% Generate Test Data

16 17 18 19

lowlt = 1e-8; % earliest time point hilt = 1e-1; % latest time point nptst = 201; % number of time points

20 21 22

t = linspace(log(lowlt),log(hilt),nptst); timeax = exp(t); % generate model time axis

23 24 25 26 27

decay_trace = exp(-(timeax./1e-6))’ + exp(-(timeax./1e-4))’; % generate % model T-jump data decay_trace = decay_trace + 0.15.*randn(size(decay_trace)); % add noise to % the model data

28 29

%% Build Rate Axis

30 31 32 33

lowl = 1e0; hil = 5e9; npts = 201;

34 35 36

x = linspace(log(lowl),log(hil),npts); rateax = exp(x);

37 38

Amat = exp(-timeax’*rateax); % this is the lapace transform matrix

39 40 41

%%

4.A Example Matlab Scripts 42

125

mu_val = logspace(-5,2,201); % values of the Tikhonov weighting factor

43 44 45 46

resid_val = zeros(numel(mu_val),1); % initialize matrices before loop norm_val = zeros(numel(mu_val),1); Tik_all = zeros(numel(rateax),numel(mu_val));

47 48

for k = 1:numel(mu_val); % loop through many test values of mu

49 50

%%

51 52

[U,S,V] = svd(Amat);

53 54 55

Dmat_temp = diag(diag(S)./(diag(S).ˆ2 + mu_val(k).ˆ2));

56 57 58

Dmat = zeros(npts,numel(decay_trace)); Dmat(1:npts,1:npts) = Dmat_temp;

59 60 61

Tik_Opt = V*Dmat*(U’*decay_trace); % this is the solution of the % Tikhonov regularization using the SVD implementation

62 63 64 65 66

resid_val(k) = norm(Amat*Tik_Opt - decay_trace,2); % store value of the % residual norm_val(k) = norm(Tik_Opt,2); % store value of the norm of the rate % distribution vector

67 68

Tik_all(k,:) = Tik_Opt;

69 70

end

71 72

%% Find optimum value of mu

73 74 75

rho = log(resid_val); eta = log(norm_val);

76 77 78 79 80

rho_p = interp1(mu_val(2:end)-0.5.*mean(diff(mu_val)),diff(rho),mu_val,... ’spline’,’extrap’); eta_p = interp1(mu_val(2:end)-0.5.*mean(diff(mu_val)),diff(eta),mu_val,... ’spline’,’extrap’);

81 82 83 84 85 86

rho_pp = interp1(mu_val(2:end)-0.5.*mean(diff(mu_val)),diff(rho_p),... mu_val,’spline’,’extrap’); eta_pp = interp1(mu_val(2:end)-0.5.*mean(diff(mu_val)),diff(eta_p),... mu_val,’spline’,’extrap’);

126 87 88


kappa_val = 2.*(rho_p.*eta_pp - rho_pp.*eta_p)./((rho_p).ˆ2+... (eta_p).ˆ2).ˆ(3/2); % this is the curvature of the L-curve

89 90 91 92 93 94

ind_t = find(kappa_val==max(kappa_val)); % find the point of maximum % curvature ind = find(mu_val>(1.*mu_val(ind_t)),1); % find the reconstruction this % optimized mu corresponds to. In some cases is may be advisable to search % for a value of mu higher than the optimized value

95 96

opt_soln = Tik_all(ind,:);

97 98

rate_distn = opt_soln;

99 100

figure(2223);semilogx(rateax,opt_soln);

4.A.2 1 2 3 4 5 6 7 8 9

% % % % % % % % %

MaxEnt rates.m

MaxEnt_rates.m - Paul Stevenson 2017 An example script for the implementation of Maximum Entropy regularization of the inverse Laplace transform of input data This script demonstrates one possible method - it has not been optimized for speed. The optimization is done using the function newtonraphson.m from the Matlab file exchange. Other optimization methods may be useful

10 11

%% Sanitize workspace

12 13 14 15

close all; clear all; clc;

16 17

%% Define global variables

18 19 20

global prior_dist eta Iexpt rates times var_expt % these variables are % passed between functions

21 22

%% Generate rate axis

23 24 25 26 27

low_l = 1e0; % slowest rate - timescale of 1s hi_l = 1e9; % fastest rate - timescale of 1ns npts = 30; % number of log spaced points between these rates - more points % slow down the execution of the script

4.A Example Matlab Scripts

127

28 29 30

x = linspace(log(low_l),log(hi_l),npts); rateax = exp(x); % this is the rate axis for the distribution

31 32

%% Load the experimental data

33 34

load(’dmpc_T19’); % load in example data

35 36

timeax = logspace(-9,-1,50); % generate time axis for model data

37 38 39

decay_trace = (exp(-(timeax/2e-3).ˆ0.65) - exp(-timeax./1e-6)) ... + 3e-2.*randn(size(timeax)); % generate example data with noise

40 41

%% Determine the value of sigma

42 43

var_expt = (3e-2); % experimental noise estimate

44 45

%% Determine the rate distribution

46 47 48

ChiSqOld = 1000; % initializing variables ChiSq = 100;

49 50 51

eta = var_expt; prior_dist = zeros(numel(rateax).*2,1)+1e-4;

52 53 54 55

Iexpt = decay_trace; rates = rateax; times = timeax;

56 57

eta_inc = 1;

58 59

count = 0;

60 61

while ChiSq>1.5 % loop to minimize ChiSq

62 63

count = count+1 % display the number of interactions

64 65 66 67 68 69

if exist(’xprev’,’var’) x0 = xprev; else x0 = prior_dist’; end

70 71 72

options=optimset(’MaxFunEvals’,100,’MaxIter’,100,’Display’,’Iter’);

128

Chapter 4. Processing and Analysis in Transient Spectroscopy [xout, ˜, F, exitflag, output, jacob] = ... newtonraphson(@finddQ_master, x0, options); xout = abs(xout); % sanitize xout - must always be positive

73 74 75 76

ChiSqOld = ChiSq; ChiSq = findChiSq_neg(xout,decay_trace,rateax,timeax,1); deltaChiSq = ChiSqOld - ChiSq; % used to determine if eta needs to be % increased

77 78 79 80 81

xprev = xout;

82 83

if deltaChiSq50% yield by rinsing the windows with a buffer containing DDM and NaCl, and concentrating the run-off with the spin-concentrator.

A final consideration with sample preparation for FTIR is minimizing the light scattered from the sample. FTIR measurements are performed in the time domain - if some particle drifts into the beam, this will introduce a non-periodic signal into the FTIR signal. This will result in distortions in the frequency-domain representation of the data. One trivial cause of scatter is the windows themselves; this is “trivial” because experiments should be done using windows which have been carefully cleaned and examined. Some more challenging sources of scattering to address are bubbles and sediment in the sample. These effects can be minimized by centrifuging the sample at low speed for10 min prior to sample preparation. This has the effect of concentrating any particulates in the solution to the bottom of the sample, and bursting any bubbles in the sample. Bubbles are a particular concern because the presence of a detergent enables them to persist for a very long time. When pipetting the KcsA solution on to the CaF2 windows, it is advisable to draw up more liquid than required, and only partially depress the pipettor. The final stage in the pipettor is designed to forcibly dispense all liquid in the pipette tip; this often produces a large number of bubbles and is best avoided.

In addition to the spectra of the protein under NaCl and KCl salt conditions, bufferonly spectra were collected to subtract the solvent background. All spectra were collected against a background of dry air at 20 ◦ C. Difference spectra presented were calculated as (KcsA+K+ , buffer subtracted) - (KcsA+Na+ , buffer subtracted)

5.2 Methods

5.2.3

137

Spectral Modeling

To facilitate the interpretation of the experimental spectra, the Local Amide Hamiltonian (LAH) picture discussed previously in §2 is combined with MD simulations. A spectroscopic map is used to translate the positions and local environments of the amide carbonyls in the simulation into a time-dependent Hamiltonian [75]. This yields the IR spectra when combined with transition-dipole trajectories, also derived from the MD simulations.

Molecular Dynamics Simulation MD simulations were performed using the GROMACS software package, and the CHARMM36 forcefield with the TIP3P water model. The construction of a membrane protein system for simulation is not straightforward, either in principle or in practice, and so bears further discussion. The protein structure was based on the crystal structure of the WT KcsA, with the Nterminal chain truncated (labeled here ∆(1-21)WT, PDB accession code 3EFF) [104]. This is the most complete crystal structure of the WT structure at the time of writing; unlike earlier crystal structures, it includes the C-terminal helices. Unfortunately, the positions of ions bound in the selectivity filter are not resolved in this crystal structure. To obtain these coordinates, we use the crystal structure of the smaller KcsA construct (where the Nand C-terminal chains are deleted, PDB accession code 1K4C) which does resolve the ion positions [99] and align the smaller structure to the larger structure using Pymol. From this, we may generate structures with different ion configurations. KcsA has four binding sites which may be occupied by K+ or water. X-ray diffraction studies have shown that at least two of the binding sites are occupies simultaneously, giving two configurations for K+ binding we must consider - KWKW and WKWK. Positions for Na+ occupancy were taken from literature reports [105]. With the protein structure created, the next step is to build an environment for it. One challenge here is the lack of forcefields for lipids and detergents; what forcefields are available are often less mature than their protein and water counterparts. Additionally, the

138

Chapter 5. Infrared Spectroscopic Insight into Ion Binding in KcsA

Figure 5.2: a) Workflow summarizing the simulation procedure. A structural ensemble is generated from MD simulations. From these structures, Hamiltonians can be constructed from the local environment of the carbonyls, and their spatial relation to one another. These Hamiltonians can be used to calculate the eigenstates of the system, and from that IR spectra and Doorway Modes. b) Schematic description of the Doorway Mode transformation.

5.2 Methods

139

aggregation of lipids and detergents into structural assemblies is challenging, particularly in the presence of a protein. The lipid 1,2-dimyristoyl-sn-phosphatidylcholine (DMPC) has a forcefield parameterization which has been demonstrated to reproduce NMR order parameter data and X-ray scattering headgroup areas [106, 107]. DMPC has a preference for forming planar bilayers; this structure is more straightforward to incorporate a protein into than a detergent micelle. The final lipid-protein structure was created starting from the end of a 500ns simulation of a DMPC bilayer (taken from Reference [106]). KcsA was embedded in this structure using the inbuilt GROMACS tool g membed [108]. For simplicity, the system is aligned such that the membrane normal and C4 rotation axis of KcsA lie along the z-axis. The system is then hydrated using the inbuilt functions in GROMACS. With the system constructed, all that remains is to run the simulation. The system was equilibrated at 293 K for 1 ns using an NVT ensemble, then for a further 1ns using an NPT ensemble. A 7 ns production simulation was run in an NPT ensemble, where the ions were subjected to an additional harmonic potential with a force constant of 4 kB T˚ A−2 at 293 K. This allows the ion to experience fluctuations without being able to drift from the specified position. For analysis, the last 5 ns of the trajectory was used. A total of 3 systems were constructed - KWKW, WKWK and Na+ ion configurations. Temperature was maintained with the Nos´e-Hoover thermostat, and pressure was held constant (where appropriate) with the semi-isotropic Parrinello-Rahman barostat.

Amide I Spectral Simulations A time-dependent amide I Hamiltonian was constructed from the MD simulations in the following way. First, protein structures were taken from the simulation at 1 ps intervals. For each structure a site energy for each amide carbonyl was determined using a spectroscopic map to convert the local electrostatic environment to a frequency [75]. This forms the diagonal elements of Hamiltonian. The off-diagonal elements (coupling between oscillators) are determined by a pairwise consideration of the distance and angles between oscillators,

140


using the formula given in §2. Diagonalization of the Hamiltonian only yields the energies of the eigenstates of the system. To calculate the FTIR spectrum, the dipole moment trajectory must be calculated from the structures. The same transformation matrix which transforms the Hamiltonian from the site basis to the eigenstate basis is used to calculate the eigenstate dipoles. The magnitude squared of these dipoles gives the intensity associated with each eigenstate. The computational cost of diagonalizing the Hamiltonian scales as O(n3 ), where n is the number of oscillators in the system. To improve the efficiency of the calculations, the Hamiltonian was initially broken into five regions - the four long C-terminal helices, and the central pore region (this is the region which was initially crystallized in Reference [99]). Spectra were calculated individually and then summed. In some cases, the Hamiltonian was further subdivided, as discussed in §5.3.2.

Doorway Mode Analysis The protein system constructed earlier has 668 residues. For each frame our of trajectory, 668 vibrational modes are calculated. Each spectrum is calculated from 5000 frames - over 3 million vibrational modes are contained within the 100 cm−1 spectral range. Developing a molecular assignment of regions of the spectrum is challenging. One route to developing a molecular interpretation of the different spectral regions is to apply a Doorway Mode Analysis (DMA) to the spectrum [77, 109]. In DMA, a frequency range of interest is defined (typically 1 nm) it is unlikely that there is any significant direct coupling, rather this mode reflects the fact that all four helices have vibrations within this frequency range. This confirms our earlier assignments of Peak 2 as arising from helices. The question of why the helix mode varies with the identity of the ion bound is discussed later in §5.4.2. The remaining three modes have a greater contribution from the selectivity filter. Peak 3 appears to be a combination of the 1674 cm−1 and 1680 cm−1 , which we label as the S3 and S2 modes (this refers to the binding site which has the most contribution to each mode). Similarly, the 1630 cm−1 mode is localized around the binding site of the Na+ ion, though with more contribution from the helix than the S2 and S3 modes. This is because the T75 residue which forms the binding site for Na+ acts as a hydrogen-bond donor to the helix adjacent; the T75 and helix carbonyls are aligned in a near-optimal geometry to interact with each other. These doorway modes are consistent with our earlier heuristic

5.4 Vibrational Assignments and Interpretation

151

analysis.

5.4.2

Interpretation

Armed with this molecular interpretation of the difference spectra, the next task is to relate the experimental spectra to the large body of literature on KcsA, with the goal of gaining some new insight into the dynamics of the pore region. The extent of the delocalization of the vibrational modes provides some insight into the dynamics of the protein. The extent of the delocalization depends on two competing factors, the coupling between oscillators (J) and the difference in site energies in the system (∆ω). Increasing J and decreasing ∆ω leads to increasing delocalization. In the case of the peak 2 (localized primarily on the α-helix; the near-uniform electrostatic environment of the α-helix minimizes ∆ω while also locking the oscillators into a geometry which gives rise to significant coupling between oscillators, resulting in a vibrational mode which is spread over the length of the pore helix. The delocalized nature of this mode is the key to why peak 2 is observed in the difference spectrum. Slight changes in the structure of α-helices, such as a slight bending of the helix from the helical axis, are predicted to change the intensity and line widths of IR peaks [111, 112]. The increase in intensity at 1660 cm−1 in the K+ bound data suggests that the α helices deviate further from an ideal α-helix geometry than the Na+ bound structure. Simulations have suggested that changes in the surrounding pore helices are thought to play a role in whether or not KcsA is in a conductive state [36]. The change in frequency and intensity of peak 2 is consistent with this idea; structural changes in the selectivity filter and the surrounding helical environment are correlated. The remaining three modes we identify are more local in character, because the electrostatic environments that arise from binding K+ and Na+ strongly, and selectively, influences the frequency of the coordinating carbonyls. These conclusions are borne out in the model. For the B mode, the proximity of Na+ lowers the frequency of the four T75 carbonyls by

152


Figure 5.7: Doorway modes for the four vibrations identified in the text. Color indicates the relative phase of the vibrational motion, intensity of color indicates the magnitude of the contribution of that site to the mode.


153

>25 cm−1 relative to the V76 carbonyls. Considering the weaker coupling between T75 and V76 carbonyls (J < 8 cm−1 ), binding Na+ in plane with the T75 residues largely decouples these vibrations from the rest of the extended filter, creating a site-specific infrared transition. For the S3 mode, the similar frequencies of the K+ -binding V76 and G77 carbonyls (∆ω ≈ 10 cm−1 ) and their strong coupling (J > 8 cm−1 ) lead to a vibrational transition that is largely localized on the eight carbonyls defining the S3 site. Similar factors seem to be responsible for the S2 mode which strongly involves the T75 and V76 residues. The confined nature of these vibrational modes offers insight into the electrostatic environment of the selectivity filter. The strong coupling between adjacent residues suggests an environment where the oscillators are held firmly in place. The significant electrostatic frequency shift oscillators adjacent to either monovalent cation experience reveals that, perhaps unsurprisingly, the ion acts as a significant electrostatic perturbation in the system. Together, these insights paint a picture of the selectivity filter as an environment which is conformationally-slaved to the electrostatics perturbation of the ion. The agreement between the simulated spectra and the experimental results also suggest that, contrary to the suggestion in Reference [32], the main occupancy state of the channel is the alternating ion-water configuration, rather than the all-ion configuration. This was also one of the primary conclusions of another IR study of KcsA published later [33]. Neither the experiments in this chapter, nor those in Reference [33] address the nature of the selectivity filter during conduction, however. More experiments are required before the suggestion in Reference [32] that all-ion configurations are responsible for ion transport can be definitively ruled out.

5.5

Summary and Outlook

Different configurations of ions in the selectivity filter of KcsA give rise to distinct spectral signatures. These signatures arise from a combination of frequency shifts from the electrostatic interaction with the nearby ions, and a change in configuration of the selectivity

154


filter. These spectral changes can be reproduced by a combination of MD simulations and spectral modeling. This demonstrates the utility of IR spectroscopy as a probe of both structure and electrostatic environment. The assignment of these modes opens up the possibility for time-resolved IR experiments, such as 2D IR or relaxation experiments, to characterize the fluctuations within the selectivity filter. Studies of these fluctuations, and what environmental changes modulate these fluctuations, have the potential to significantly enrich the current understanding of how ions are conducted through potassium channels. The vital interaction between the carbonyl groups and the ions during conduction makes IR spectroscopy uniquely well-suited as a probe.

5.1 Protein Expression

5.A

155

Protein Expression

It is useful to highlight some of the design principles guiding efficient protein expression, and clarify what is meant by terms such as “Induction”. Native protein expression levels are carefully regulated in all living cells, and play a significant role in allowing living organisms to adapt to changes in their surroundings. One example of this is upregulation (increase in expression levels) and downregulation (decrease in expression levels) of insulin in response to changes in glucose levels; an inability to maintain appropriate levels of insulin is the cause of Diabetes Mellitus. Returning to the E. coli system, it is clear that high levels of expression of a non-native protein will adversely affect the viability of the cells. Therefore, to maximize the total yield of the desired protein, it is desirable to allow the bacteria to grow naturally to an optimum size, before triggering (or inducing) the bacteria to start producing our protein. This induction step can be realized with careful design of the plasmid containing the DNA for the non-native protein. Specifically, preceding the region of DNA which promotes expression, a sequence known as the lac operon is inserted. This region binds tightly to a protein known as the lac repressor, which inhibits the transcription of this plasmid. On addition of certain types of sugar, the lac repressor undergoes a conformational shift which decreases the affinity for the lac operon, freeing the plasmid for transcription. After induction, this is an optimum period of time to allow the cells to produce the desired protein, before harvesting the protein from the cells. It is important to remember that the principle behind protein overexpression is to hijack the cellular machinery of an organism to perform some otherwise-challenging chemical syntheses; however, this cellular machinery also includes an extensive waste disposal system. In particular, proteases can significantly reduce the yield of intact protein, so it is important to optimize the delay between induction and harvesting, and to also include a protease inhibitor cocktail in any extraction medium.

156


Chapter 6

Multidimensional Spectroscopy of Gramicidin D The contents of this chapter are modified from the following publication: • P. Stevenson and A. Tokmakoff, Distinguishing gD conformers through two-dimensional infrared spectroscopy of vibrational excitons, J. Chem. Phys. 142, 212424 (2015) [114]

6.1

Gramicidin D

Gramicidin D is a naturally occurring linear pentadecapeptide which can form transmembrane pores, both in vivo and in vitro [23,115]. The rare combination of biological function and tractable system size in gD have resulted in an extensive body of literature devoted to studying the structure and mechanism of function of gD [116–119]. Despite the small size of the system, the history of gD is steeped in controversy; long after gD was crystallized in organic solvents, the structure of the ion conducting state in membranes remained contentious despite a great number of solid-state NMR studies [120–122]. The primary reason for this debate is the existence of many subtly different conformers of gD, which are predominantly helical but differ in pitch and dimer contacts [115]. These conformers are close enough in 157

158

Chapter 6. Multidimensional Spectroscopy of Gramicidin D

Figure 6.1: Left, Helical Dimer conformation of gD dimer. Right, Dimeric Helix conformation of gD dimer. Amide carbonyl groups are shown by spheres. energy that the nature of the solubilizing lipid or detergent, or even the solvent the sample was initially prepared in, can determine the structure of the sample. Gramicidin D is also rather unique in that the amino acids have alternating l-d stereochemistry [123]. This results in an unusual hydrogen-bonding pattern. Rather than all the amide groups pointing in the same direction along the helical axis, as with α- and 310 helices, the groups alternate direction, giving rise to a hydrogen-bonding pattern similar to β-sheets. For this reason, these structures are usually referred to as β-helices, though it is necessary to distinguish this usage from the parallel β-helices commonly found in antifreeze proteins. Three variants of gD occur naturally and are referred to A, B, and C. They differ only in a single point mutation at residue 11, which is a tryptophan in gramicidin A, a phenylalanine in B, and a tyrosine in C. These mutations are known not to alter the structure or function of gD [25]. Two distinct dimeric species of gD are known to exist in lipid bilayers. They are the Helical Dimer (HD), a head-to-head dimer composed of two β-helices which have parallel

6.2 Sample Preparation

159

β-sheet contacts in the monomer unit and antiparallel β-sheet contacts at the interface [25], and the Dimeric Helix (DH), a double-stranded helix with antiparallel β-sheet hydrogenbonding contacts [124], as shown in Figure 6.1. Though these structures can be distinguished with Circular Dichroism (CD) spectroscopy [125], it has also been observed that differences in the amide I IR spectra can distinguish between the different conformations [126] (shown in Figure 6.2). Changes in the oligomerization state and subtle structural deformations of the HD structure are thought to control the conductive properties of the channel. The ability to map changes in the IR spectrum on to structural changes in the protein offers a new route studying the functional dynamics of gD, while also providing a new test for the predictive power of amide I spectral modeling. Though the HD and DH structures are known to have different IR absorption spectra, this provides only a limited test for a spectral model. Consider the case where an absorption spectrum shows two clear peaks. Do these peaks arise from two different conformations of the protein, or does one conformation of the protein give rise to two IR-active vibrational modes? 2D IR provides a more stringent test of any model of the vibrations of a system. The lineshapes in the 2D IR spectrum report on the femtosecond dynamics of the system, while crosspeaks reports on coupling between modes. A model which satisfies the constraints of the 2D IR spectrum can be used to interpret data from novel, unknown structures with confidence.

6.2 6.2.1

Sample Preparation Conformational Control of Gramicidin D

Gramicidin D (VGALAVVVWLXLWLW where underlined residues are d-isomers and X = W, F or Y for gramicidins A, B, and C, respectively) was purchased from Sigma-Aldrich as a natural product and used without further purification. Labile protons in the protein were exchanged for deuterons by heating gD in d1-ethanol (EtOD, Cambridge Isotopes) at 313 K for 24 h. The solution was then lyophilized to provide a white powder which served

160


Figure 6.2: FTIR spectra of gD in the HD (left, SDS detergent) and DH (right, TX100 detergent) conformations. Regions mentioned in the main text are highlighted. Insets show the CD spectra confirming the secondary structure type of each sample. as the starting material for further sample preparation steps. The conformational preference of gD in detergent micelles can be controlled by judicious choice of detergent identity [25, 125]. This is in contrast to the case for gD intercalated into large lipid structures (such as multilamellar vesicles) where the solvent history of gD plays a significant role [127]. The HD conformation can be prepared by using the anionic detergent Sodium Dodecyl Sulfate (SDS) to form micelles around gD [25]. This is accomplished by dissolving gD in d3trifluoroethanol (d-TFE, Cambridge Isotopes) to form a concentrated solution (50 mM). To 100 µL of this solution, 500 µL of 1M SDS (Fisher Scientific) in D2 O was added dropwise while the solution was agitated. To this solution, a further 500 µL of D2 O was added, before the solution was lyophilized and rehydrated with buffer. During this protocol, it is important to ensure the solution does not become turbid. Gramicidin D is insoluble in water; it will form large aggregates which are visible as white particulate in the solution. These aggregates are kinetically stable, even in the presence of significant concentrations of detergent. The DH conformation is prepared by using the non-ionic detergent Triton X-100 (TX) as the solubilizing agent [125]. Gramicidin D was dissolved in EtOD at a concentration of 20 mM. To 300 µL of this solution, 300 µL of TX was added. This solution was vortexed for


161

several minutes to ensure mixing; TX is approximately 100 times more viscous that water at room temperature. Gramicidin D dissolves completely in this ethanol/TX solution. This solution was then diluted to a total volume of 2.5 mL, then lyophilized overnight. Both the SDS and TX solubilized samples have approximately 100 detergent molecules per gD molecule. Empirically, it was found that using a smaller ratio of detergent to protein resulted in aggregation and sedimentation of gD, resulting in a solution with a detergent-toprotein ratio of approximately 100:1, but lower protein concentration.Samples containing higher detergent concentrations than those used here proved unsuitable for preparation of high-quality optical samples. In particular, the formation, and persistence, of bubbles in detergent-containing samples is a significant problem for 2D IR. This results in light from the incident mid-IR pulses to be scattered into the detector, causing significant spectral distortions and, in some cases, saturating the detector response. The lyophilized samples of gD in SDS or TX were rehydrated in a deuterated buffer (50 mM sodium phosphate, pD 7.5) to an approximate monomer concentration of 5mM. This concentration yielded samples with an optical density (OD) of ≈0.1 from the protein amide I peak. Samples with high OD can result in distortions in the 2D IR spectrum, as the generated signal field can be reabsorbed by the sample. The rehydrated samples were heated for 6 h at 323 K. This is necessary for the TX samples to ensure the HD form is entirely in the left-handed DH structure; several kinetically stable conformations are initially formed which convert to the left-handed DH structure on heating. Literature reports suggest the sample should be exclusively composed of the lefthanded DH conformation after 1 h heating at 323 K. Circular dichroism spectroscopy was used to confirm the structure of gD in micelles from each preparation. Samples were diluted to give an A280 of 0.5 in a 1 mm-pathlength, strain-free quartz cuvette (total gD concentration 500 µM). Data were collected at 1 nm resolution on a Jasco J-1500 CD spectrometer and averaged over 5 data sets. Controlling the conformation of gD through sample preparation requires different approaches depending on whether a lipid vesicle/bilayer or detergent micelle environment is

162


the final goal. Detergent micelles require procedures like those outlined above; there is generally one stable conformation of gD in each type of micelle. Lipid envionments, however, are different. Many different conformations of gD can be prepared in the same lipid environment by changing the identity of the solvent the gD was originally dissolved in. These differences give us some insight into the differences between vesicle and micelle environments; in the detergent environment, the protein is conformationally free enough to find the minimum-energy conformation (or a set of similar-energy minimum conformations). In lipid environments, however, the protein is much more restricted in its ability to explore conformational space; a single dominant conformation of gD can be produced in lipid vesicles, but typically only after heating the sample for extended periods of time. This suggests the gD molecule inserts into the lipid environment in a conformation determined by the solvent, and is then kinetically trapped in this (or a related) conformation.

6.2.2

Optical Sample Preparation

Careful sample preparation is a key part of any spectroscopic experiment, though the mundane details of this are rarely discussed. It is of some utility to discuss, at least in brief, some useful protocols for preparing high-quality samples. These steps must always begin with a careful consideration of the possible sources of scatter and/or impurities in the sample. Using a microcentrifuge to spin the sample for >5 min at >10,000 G is often helpful for removing scatter. Care must be taken, though, when deciding which part of the centrifuges sample to extract from; some impurities may sediment out at the bottom of the sample container, while others (such as lipid aggregates) are actually less dense that water, and will rise to the top of the sample. With large enough sample volume, extracting the purified sample directly from the center of the container with a syringe may be desirable. (This is of particular use when preparing vesicle samples, such as §8.)


163

Window Preparation Additional steps to minimize the scatter from a sample focus on the preparation of a small volume of sample sandwiched between two CaF2 windows separated by a PTFE spacer. Though it may seem trivial, using clean windows is essential; visible inspection of the windows is typically insufficient to assess their cleanliness. An effective washing procedure for preparing lipid and/or detergent containing samples is to first pre-rinse the windows with a solvent based on the last sample used on the windows (e.g. if an aqueous sample was used water is suitable, or if an organic sample was used chloroform may be appropriate), then rinse and dry the windows using water, methanol, water, acetone, and water again. Empirically, windows prepared in this way result in samples with fewer “unexplained” aggregation or scatter events. If a 50 µm thick spacer is used, a sample volume of 25 µL is sufficient to completely fill the sample cell; however, drawing more liquid than necessary (≈ 40 µL) with the pipette and not fully depressing the push button results in fewer bubbles in the sample. The last stage of liquid ejection from the pipette is more forceful, to ensure all the liquid is fully ejected; this is undesireable in this case. Finally, sandwiching the sample by first allowing the two windows to touch at one point, and closing the sample in a “hinging” motion has the advantage that any bubbles which are produced are pushed to the edges of the sample, out of the region probed by the mid-IR beams.

Pump-probe Spectra Well-prepared samples will typically have scatter which is virtually undetectable (> 1000 : 1 ratio) relative to the 2D IR signal. Even in the case where more significant scatter is present, only the component of the scatter which oscillates in t1 is of concern; the time-domain implementation of 2D IR necessarily incorporates a Fourier-filter into the data analysis. However, acquiring a pump-probe measurement (to correct the 2D IR spectrum for stagemotion errors) can be troublesome. Since only a small fraction of the light routed into the boxcar interferometer is used for the pump-probe measurement, even a small amount of

164


scatter can be comparable to the signal levels. This problem is compounded by the need for accurate lineshapes in the pump-probe spectrum if the timing errors in the stage are to be corrected for. One remedy is to have the probe stage oscillate back-and-forth with a period which is not a sub-harmonic of the laser repetition rate (in practice, 13 Hz was found to be adequate) and an amplitude equivalent to a time delay of ≈1 cycle of the center frequency of the range of interest [128]. This has the effect of integrating over a short window of the pump-probe signal and scatter, removing the oscillating scatter signal. This window must be sufficiently short such that no significant changes to the lineshapes of the pump-probe spectrum are expected.

6.3 6.3.1

Spectroscopy of Gramicidin D FTIR Spectra

Buffer-subtracted FTIR spectra for the HD and DH conformations of gD are shown in Figure 6.2. The DH spectrum contains an intense peak centered at 1610 cm−1 . This feature is from the ring-mode vibrations of the aromatic groups in the TX detergent. The oscillator strength of this mode is actually rather weak, but the high concentration of TX relative to gD results in comparable intensity to the protein peaks. The low oscillator strength of this mode, however, means it appears as a very weak peak in the 2D IR spectra. Though not the primary focus of this chapter, the TX detergent shows interesting vibrational dynamics, which are presented in §6.A. Though both conformations are helical in nature, neither of the spectra resemble that of an α-helix (a broad Gaussian peak centered at 1650-1660 cm−1 ). Instead, both spectra resemble the FTIR spectrum of a β-sheet. The HD and DH spectra feature an intense peak centered at 1620-1630 cm−1 , and lower intensity peaks at frequencies >1660 cm−1 . To assist with the discussion of the gD spectra, we define three regions in each spectrum - region 1 covers the intense low-frequency peak at 1630 cm−1 , region 2 covers the range 1630 cm−1 to

6.3 Spectroscopy of Gramicidin D

165

1670 cm−1 and region 3 covers the range 1670 cm−1 and above. The most significant spectral differences between the two conformations are found in region 2; the DH conformation shows a smooth slope in this region, while the HD conformation has a distinct peak. The spectrum of the DH conformation shows a remarkable similarity to that of poly-l-lysine at high pH, where extended β-sheet structures are formed [129].

6.3.2

2D IR Spectra

The features noted in Section 6.3.1 are much more pronounced in the 2D IR spectra. The spectral differences between the HD and DH conformations are much less subtle, as shown in the 2D IR spectra in Figure 6.3. Both conformations show a clear crosspeak between regions 1 and 3; this is in stark contrast to the 2D IR spectra of other helical structures, emphasizing how sensitive amide I spectroscopy is to the spatial arrangements of the amide carbonyls [130, 131]. Differences of the amplitudes of various features in the FTIR and 2D IR spectra are informative; each signal has different scaling with dipole strength (µ2 and µ4 , respectively). This demonstrates that the peak in region 2 of the HD structure, which is almost absent in the 2D IR spectra, arises from a large number of weak (low µ) vibrational modes. This same effect is responsible for the dramatic reduction in intensity of the TX-100 mode in the DH spectrum; this mode is the most intense peak in the FTIR spectrum, but is less intense than the weak high-frequency peak in the 2D IR spectrum.

HD Conformation The lowest and highest frequency peaks observed in the FTIR spectra are clearly resolvable in the 2D spectra, but the peak found in region 2 is not visible along the diagonal. However, the effect of this peak is seen in the 2D spectra in the form of an upper-left quadrant crosspeak between the region 2 peak and the intense region I peak (ω1 , ω3 = 1625 cm−1 , 1660 cm−1 ). The 2D spectrum collected with ZZYY polarization reveals a crosspeak to the high-frequency peak in region 3, evidenced by the vertical elongation of the regions 1 and

166


Figure 6.3: 2D IR spectra of gD. a) ZZZZ polarization of HD gD in SDS micelles. b) ZZZZ polarization of DH gD in TX micelles. c) ZZYY polarization of HD gD in SDS micelles. b) ZZYY polarization of DH gD in TX micelles.

6.3 Spectroscopy of Gramicidin D

167

2 crosspeaks. The effect of the crosspeaks may also be noted in the lower-right quadrant, but the effect is less pronounced due to spectral congestion, with crosspeaks appearing as tails rather than distinct features. The 2D IR spectrum shows only limited change in crosspeak shape and intensity with changing polarization conditions. The only discernable change is an increase in crosspeak intensity to the blue side of the intense peak (growth of intensity at ω1 , ω3 = 1635 cm−1 , 1660 cm−1 switching from ZZZZ to ZZYY). This suggests that the peaks highlighted previously have transition dipoles which lie approximately parallel to each other.

DH Conformation Similar to the helical dimer spectra, crosspeaks are observed between the peak in region 1 and the high-frequency peak in region 3. Notable, however, is the crosspeak between the continuum of region 2 and the high-frequency peak in region 3. This feature, though less distinct than the 1-3 crosspeaks, seems to be centered on 1645 cm−1 , suggesting the continuum feature may not be simply a distribution of conformations, as has been suggested in the literature for poly-l-lysine, but rather a distribution of vibrational modes arising from a single conformation. Unlike the HD spectra, the crosspeaks observed (ω1 , ω3 = 1625 cm−1 , 1675 cm−1 ) do show a significant dependence on the polarization scheme used. Significant enhancement in the relative intensity of the crosspeak is observed switching from ZZZZ to ZZYY, a feature shared by β-sheet spectra. Crosspeaks between perpendicular dipoles are expected to show this polarization dependence.

Transition Dipole Angles One advantage of 2D IR is the ability to resolve angles between the transition dipoles of vibrational modes. The excitation of the sample creates an anisotropic distribution, even if the sample is initially isotropic (as is the case with micelles), so the signal will depend on

168


Figure 6.4: Heatmap of the depolarization ratio spectra of the HD (top) and DH (bottom) conformations of gD. the angle between transition dipoles contributing to a particular peak and the polarization of the pulses. Practically, this means that by comparing the relative amplitudes of the crosspeaks in parallel and perpendicular polarizations, it is possible to determine the angle between transition dipoles. If we define the depolarization ratio as the intensity ratio of a crosspeak between polarizations ρ = IZZY Y /IZZZZ then the transition dipole angle is obtained from [132]:

cos2 θ =

7 − 6ρ . 1 + 12ρ

(6.1)

This expression only constrains θ between 0 and π/2, since the cos2 θ term precludes the unambiguous determination of an angle with respect to its supplement. Figure 6.4 shows a 2D heatmap of the depolarization ratio and corresponding θ calculated with Eq. 6.1 for the upper left quadrant of the HD and DH 2D IR spectra. A comparison of the heatmap with the ZZYY crosspeak intensity contours shows a rather surprising result. The angle heatmap for the DH structure lines up closely with the crosspeak intensity and serves to highlight the underlying structure, such as the regions 2 and 3 crosspeaks discussed earlier. From the depolarization ratio, an angle θ ≈ 60

◦

between

the transition dipoles for the 1630 cm−1 and 1670 cm−1 modes is predicted. In the case of

6.4 Bloch Model

169

HD, the angle heatmap maximum does not correspond to the crosspeak intensity maximum but is shifted by approximately 10 cm−1 along ω1 to higher frequency. This observation indicates that there is substructure within region 1 for the HD spectrum, presumably as a result of overlapping transitions with a large dipole projection angle between them.

6.4

Bloch Model

From the 2D IR spectra, it is clear that the helical dimer and dimeric helix conformations are distinguishable with IR spectroscopy, not simply from peak frequencies, but from the number, strength and coupling, and dipole angles between different absorptions. This is perhaps an unexpected result; both structures are dimeric, forming helical pore-like structures of similar lengths (≈20 ˚ A for HD and ≈25 ˚ A for DH) stabilized by β-sheet hydrogen-bonding patterns. This raises the exciting possibility that we may have enough distinct spectral markers to distinguish subtle conformational changes within these different conformations. However, detailed assignment of the molecular origins of the amide I vibrational modes is necessary for this approach to succeed. Calculations of IR spectra fall into one of two broad categories; either atomistic detail (usually from MD simulations) is used in conjunction with a spectral “map” to generate a time-dependent amide I Hamiltonian and transition dipole trajectory, from which a spectrum may be calculated, or an idealized model of the system is constructed (often constructed such that amide oscillators have a simple geometric relation to one another) which can be solved analytically. The atomistic approach has the advantage of generality; if MD simulations can be run on the system, 2D IR spectra can in principle be calculated. (In practice, the rate-limiting step is diagonalization of the two-quantum Hamiltonian, which limits the size of the system which can be studied.) The analytical approach has the advantage that, where it is possible to construct an appropriate model, it is often possible to gain significant insight into the molecular origins of the vibrational modes. Atomistic approaches to spectral simulations

170


Figure 6.5: a) Helical Dimer structure (PDB accession code 1JNO), and schematic representation of the model helix.(b) Dimeric Helix structure (PDB accession code 1ALX). The model helix is unwound, and grey bars show the connections broken by the projection onto two dimensions. Couplings referred to in the main text are shown on the schematic. are effective at demonstrating whether a particular structure is or is not consistent with the experimental spectral data, but explaining why a structure gives rise to a particular spectrum requires great effort. The small size of gD and simple secondary structure makes it an attractive candidate for an analytical model. In the case of periodic structures, it is possible to develop analytical models based on a Bloch-wave approach. This has been used to a great effect with amide I spectroscopy of βsheets to explain the observed peak frequencies and intensities in terms of their vibrational modes. Two idealized helices consisting of CO groups aligned along the helical axis (z) may be used as the idealized structures for the HD and DH conformations for a Bloch-wave model. In this model, the helical structure does not represent the connectivity of the CO groups via the polypeptide backbone but rather the spatial arrangement of the individual amide I transition dipole moments for each peptide linkage. By defining a screw operation (a rotation, φ, followed by a displacement ∆) which connects adjacent sites (following the work of Loxsom on optical rotation [133]), we may reduce the problem to a pseudo-1D

6.4 Bloch Model

171

problem. This is possible because of the property {Cz (φ)|∆}r = {Cz (rφ)|r∆} which relates any two sites in a regular helix. Now, the parameters which define our helix are the distance between rungs, the number of residues per turn (γ), and the number of residues in the helix (N ). It is possible to take one of the two equivalent approaches to building this idealized model. In the standard approach, a unit cell consisting of an adjacent pair of antiparallel CO groups which are replicated around the helix is defined. Alternatively, the repeating unit may be defined as a single oscillator at each site, where the “Up-Down” alternating nature of the CO groups is introduced afterwards in the transition dipole when calculating intensities. The former case results in a helix exciton band split into sub-manifolds, corresponding to the replication of the symmetric and asymmetric modes of the unit cell. The other approach has the advantage of simplifying the interaction between different “rungs” of the helix, since we discuss coupling between individual CO groups. Because a simple interpretation of coupling between oscillators is important below, we choose to define our repeating unit as a single amide unit. The next step is to generate a simplified Hamiltonian for the system. The strongest couplings in the system will be between hydrogen-bonding partners and covalently linked groups. A schematic summary of the model is given in Figure 6.5. The symmetry of our system requires only two coupling constants to be specified; J1 , the coupling between covalently linked amide groups (i.e. nearest neighbors), and J2 , the coupling between hydrogen-bonded partners on adjacent rungs. The matrix elements of our Hamiltonian are then given by:

Hij = ω0 δi,j + J1 δi,j±1 + J2 δi,j±γ

(6.2)

where ω0 is the site energy of the isolated amide I group and i, j take on values from 1 to N . For periodic systems, H1,N = J1 , while for finite systems H1,N = 0. With these expressions, it is possible to write down the vibrational wavefunctions and

172


energies for a system with periodic boundary conditions (PBC) [134, 135]:

Ψpm

N −1 1 X =√ exp (ikm an)|ni N n=0

Em = ω0 + J1 cos (km a) + J2 cos (γkm a)

(6.3a)

(6.3b)

where |ni is the nth amide oscillator and a is the displacement unit of the screw operation. The exciton wavevector is

km =

2πm Na

(6.4)

and the states are indexed such that:

−

N −1 N 10 ns. However, the lack of anisotropy decay also indicates that the individual detergent molecules are static on the picosecond timescale - no appreciable “wobbling” occurs. The lack of anisotropy decay is also consistent with the idea of vibrationally isolation - proteins

182


Figure 6.8: Left, structure of Triton X-100, where n=9-10. Right FTIR spectrum of Triton X-100 in D2 O with and without Gramicidin D present. carbonyls, which should be similarly static, display a < 1 ps anisotropy decay because of the large number of coupled states in the system. At present, the TX ring modes are an interesting system to demonstrate some principles of 2D IR - such as the existence of all-coherence pathways. (These pathways are rarely discussed in the 2D IR literature, but are vitally important. If they were truly neglected, the expressions for determining the angle between vibrational modes from the depolarization ratio would be significantly different.) However, the Tyrosine sidechain shows very similar spectroscopic behavior to the TX ring modes - perhaps unsurprising given their chemical similarity. This raises the intriguing possibility of using the Tyrosine ring mode as a vibrational probe; at present, the aromatic ring modes in this group seem to be insensitive to the local environment (hence the long coherence times). However, if some environmental factor is found to alter the sidechain vibration, this may make an excellent probe precisely because of its insensitivity to other factors (such as solvent exposure).

6.1 Vibrational Dynamics of Triton X-100

183

Figure 6.9: Top row - pump-probe measurements for ZZZZ and ZZYY polarizations of Gramicidin D / TX100 samples. Middle row - pump-probe trace and residual oscillations of the TX100 1515 cm−1 mode. Bottom left, Fourier Transform of oscillations - note the peak at 100 cm−1 . Bottom right, pump-probe anisotropy of the 1515 cm−1 mode

184


Chapter 7

High-Amplitude Electric Field Fluctuations in Lipid Membranes The contents of this chapter are modified from the following publication: • P. Stevenson and A. Tokmakoff. Ultrafast Fluctuations of High Amplitude Electric Fields in Lipid Membranes, J. Am. Chem. Soc., 139, 4743-4752 (2017) [137]

7.1

Introduction

The function of integral membrane proteins is influenced by their dynamic interactions with the solvating lipid environment [138–141]. For example, the membrane fluidity of the lipid environment surrounding ATPase modulates the enzymatic activity of this protein [142]. Similarly, it has been found that the presence of cholesterol in lipid bilayers, and the effect it has on lateral diffusion of lipids, modulates the time scales of the rhodopsin photocycle [143]. In these examples, lipid bilayers provide a bath of fluctuating forces, which in turn influence protein conformation and function. These forces originate in a hierarchy of coupled dynamics that operate on varying length scales [144] and over many decades in time, from ultrafast short-range headgroup fluctuations and lipid chain isomeraizations to slower diffusive motions and capillary waves that modulate membrane thickness and curvature. 185

186

Chapter 7. High-Amplitude Electric Field Fluctuations in Lipid Membranes

The forces that act on membrane proteins arise, directly or indirectly, from the influence of electric fields. These forces originate in the motion of charged groups and induced dipoles within this heterogeneous and anisotropic electrostatic environment and are of particular importance in understanding the mechanisms of voltage-gated processes [145]. The functional role of different tiers of dynamics in biological systems has been considered for decades in the context of protein conformational dynamics [8, 10, 146, 147]. When fast and slow dynamics are coupled to each other [148], sub-nanosecond fluctuations, though much faster than typical functional timescales, may prove critical for determining the outcome of biological processes. Enzymatic catalysis is a notable (and controversial) example in which the specific nature of fast conformational fluctuations may dictate enzymatic catalysis rates [149, 150]. Specifically, the picosecond dynamics of the P450cam enzyme with different substrates correlated well with the functional behavior of the system [151]. The regioselectivity of the enzyme was found to decrease as the picosecond dynamics became faster [152]. Conduction of ions through ion channels is another case where fluctuations of the protein are functionally essential. Similarly, simulations have demonstrated that a microsecond conformational change between two helical structures can be regulated by the picosecond fluctuations in hydrogen bond strength in the system [8]. A more general argument for the importance of sub-nanosecond membrane fluctuations is their role in the microscopic friction experienced by nanoscopic objects within the membrane. Diffusive motion and activated kinetics in a particular region of space will depend on the magnitude and correlation times of the fluctuating electric fields generated by lipids, protein, and water. In this way, the kinetics of short-range membrane-associated processes are intimately tied to the fast dynamics of the membrane. Further support for this perspective requires experimental studies of fast membrane dynamics. Although small angle x-ray or neutron scattering and NMR/ESR experiments can access fast time-scales, experimental observation of electric field fluctuations in membrane proteins is a challenging proposition. Infrared (IR) absorption spectroscopy can be used to interrogate the electrostatic environments in the membrane [153], as the frequencies of func-

7.1 Introduction

187

Figure 7.1: Gramicidin D in a lipid bilayer, with ester carbonyls highlighted as small spheres, and two Tryptophan sidechains shown as large spheres. The structure of DMPC is shown in the inset. Bottom right shows the frequency shift from a test charge of +1e as a function of distance using the spectral map discussed in the text.

188


tional groups in lipid membranes (such as ester carbonyls) are sensitive to the local electric field. Vibrational Stark spectroscopy experiments have determined quantitative values for the shifts in carbonyl vibrational frequencies in response to the applied electric field, finding values of η ≈ 1 cm−1 /(MV/cm) [154, 155]. Quantitative relationships between electric field and vibrational frequency have been put to extensive use in IR spectroscopy of amide I vibrations in proteins, where experimental spectra can be linked directly to structure and dynamics in MD simulations. However, it is only with the advent of ultrafast time-resolved spectroscopies that direct experimental insight into the fluctuations of these environments has been possible. A combined experimental and simulation approach offers the potential to understand both the timescales and magnitudes of fluctuating electric fields in membranes. Femtosecond two-dimensional infrared spectroscopy (2D IR) can be used distinguish different electrostatic environments and their dynamics. 2D IR tracks time-dependent changes in the vibrational frequency and amplitude of bond vibrations, allowing one to target local regions of large molecules without the need to tag the sample with an exogenous label. Tracking the evolution of these intrinsic labels on an ultrafast timescale can reveal intricate details about the nature of the interactions within the system, and has been utilized to study a diverse range of problems, spanning liquid water dynamics, membrane protein structure and protein-ligand interactions [156–158]. The relaxation of carbonyl lineshapes in 2D IR characterizes the frequency-frequency correlation function (FFCF) for vibrational frequency fluctuations about their equilibrium value, Cωω = hδω(t)δω(0)i [159]. Since vibrational frequency shifts are proportional to changes in local electric fields, δω ≈ M δE, this experiment is effectively a measurement of the time-correlation function for the fluctuations of electric fields about their average: CEE = hδE(t)δE(0)i , or equivalently, the spectrum of ultrafast fluctuating electrical forces. With these force correlations, the microscopic friction in the system can be determined through the fluctuation-dissipation theorem. Here, 2D IR spectroscopy is combined with MD simulations to resolve the amplitude and timescales of electrostatic fluctuations proximate to the ester carbonyls in a model

7.2 Materials and Methods

189

phospholipid bilayer. IR and 2D IR spectra are used to identify hydrogen bonding configurations to the lipid ester carbonyl and the time-dependence of 2D IR lineshapes quantify the fluctuations in vibrational frequency and electric field experienced by these environments. Additionally, the modification of these fluctuating fields by membrane proteins in the bilayer is studied by introducing the ion channel Gramicidin D (shown in Figure 7.1) into the membrane. Gramicidin D (gD) acts as a non-specific defect in the membrane, but also offers the ability to look at specific hydrogen bonds between sidechains and lipid ester groups. MD simulations provide the molecular basis for interpreting the observed fluctuating fields.

7.2

Materials and Methods

1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC, Anatrace Inc) was dissolved in chloroform (Sigma-Aldrich) at a concentration of 10 mM. For samples containing Gramicidin D (Sigma-Aldrich), an appropriate amount was added to give a final concentration of 0.5 mM. Before addition, gD was first dissolved in EtOD (Cambridge Isotopes) to exchange labile NH groups for ND groups, then lyophilized overnight to remove solvent. This yields a 20:1 lipid-to-protein ratio, which minimizes the fraction of DMPC molecules not in contact with gD without forming the hexagonal lipid phase known to occur at high gD density [160]. The lipid-chloroform solution (100 µL) was dried on a CaF2 window. This sample was then placed under vacuum for 24 hours to remove any residual solvent. A small amount (≈500 nL) of D2 O (Cambridge Isotopes) was pipetted on to the sample, to a water-to-lipid ratio of 10:1. It is necessary to use D2 O rather than H2 O because of the absorption of the H2 O bending mode at 1650 cm−1 . This sample was heated to above the phase transition (Tm =24 ◦ C, sample held at 40 ◦ C) and the procedure outlined in Reference [161] was used to prepare aligned bilayer stacks. Briefly, the lipid bilayers were sandwiched between two CaF2 windows at 40 ◦ C, and a mechanical shearing force was applied by slowly rotating the upper window. Aligned layers of bilayers were necessary to minimize signal from scattered light in the sample. Using the extinction coefficients for bulk D2 O [162], and assuming the

190


ester group extinction coefficient is similar to that of the amide I vibration [163], the FTIR spectra are consistent with a 10:1 water-lipid ratio.

7.2.1

IR Spectroscopy

FTIR spectra were collected on a FTIR spectrometer (Tensor-27, Bruker), with 128 averages per spectrum against a background of dry air. The low water content of the samples removes the need for further background subtraction. For FTIR and 2D IR measurements, the sample was placed in a temperature-controlled jacket at 40 ◦ C, which is 16 ◦ C above the gel-fluid phase transition of DMPC [164]. This temperature was chosen ensure both samples are completely in the fluid phase, since gD has been observed to broaden the gel-fluid coexistence range of DMPC [164]. The instrument and methods used for 2D IR data acquisition have been described in §3. All data was collected with perpendicular (ZZYY) polarization geometry in a temperature-controlled cell held at 40 ◦ C. 2D IR data was acquired for waiting times from 0.15-10 ps.

7.2.2

Molecular Dynamics and Spectral Simulations

Since gD is a mixture of Gramicidins A, B and C, the dominant species of the mixture, Gramicidin A, was selected for simulation. (For ease of comparison with experiment, the species being studied is referred to as gD, even though the simulation is not performed for the mixture). Coordinates were taken from the solution NMR structure (PDB 1JNO) [25]. Initial bilayer assemblies were generated using the CHARMM-GUI, with 64 DMPC lipid molecules per leaflet [164–166]. A six-stage equilibration procedure outlined in Ref [165] was used, for a total equilibration time of 375 ps. Production MD simulations were run for 10 ns in an NPT ensemble with the Nos´e-Hoover thermostat with separate temperature coupling groups for the protein, lipid and solvent. The last 5 ns of this run were used for analysis. All stages were run at 313 K, using the CHARMM36 forcefield.

7.4 Results and Discussion

7.3

Results and Discussion

7.4

FTIR and 2D IR Spectroscopy of Lipids

191

Although FTIR and 2D IR spectra of the DMPC bilayers contain a wealth of information, it is contained within a congested, spectral region (spanning > 80 cm−1 ). To analyze the dynamical information in the spectra, it is imperative to understand and assign the underlying peak structure. The DMPC ester carbonyl stretch resonance in the FTIR spectrum observed at 1700-1750 cm−1 has an asymmetric lineshape, and inspection of the second derivative of this feature reveals two broad, overlapping peaks centered at 1728 cm−1 and 1743 cm−1 (Figure 7.2), which are labeled Peak A and Peak B for convenience. The presence of two peaks demonstrates that there are two distinct electrostatic environments experienced by the lipid ester carbonyls. These peaks have been the subject of extensive investigation [167, 168], in particular with the help of site specific

13 C=O

isotope

labeling of the lipid ester headgroups [169, 170]. These studies revealed that a single label results in a ≈44 cm−1 isotope-shift to both peaks, but reached different conclusions regarding the origin of this shift. Blume et al. argued that these peaks resulted from variation in the hydrogen bonding configurations of their ester carbonyls with solvating water [169], while others suggested that the different lipid chains (sn-1 and sn-2) are spectroscopically distinct. The 15 cm−1 splitting between Peaks A and B is consistent with the 16 cm−1 frequency shift observed in the amide I carbonyl stretch in proteins upon formation of a hydrogen bond to water [171]. The IR spectra reflect differences in hydrogen bond configuration to the ester groups from water that penetrates the headgroup region of the lipid bilayer. Carbonyl vibrations act as a sensitive reporter of local electric fields or electrostatic potential. Current evidence indicates that that amide [75], ester [172], and ketone [155] carbonyls are all well described through a linear relationship between their vibrational frequency and the electric field that the surroundings projects along the C=O bond axis. In addition to explaining the sensitivity of the C=O vibrational frequency to hydrogen bonding,

192


this finding is the key to developing interpretive computational spectroscopy tools that relate IR spectra to local molecular structure using MD simulations. Here, it is possible to leverage the extensive literature on amide I vibrational spectroscopy to extend our computational analysis with MD simulations to lipid ester carbonyl stretches. The 2D IR spectrum of the ester carbonyl region at early waiting times (T=150 fs) in Figure 7.3 shows a doublet of peaks, arising from the ground state (ν=0) to excited state (ν=1) transition in red, and the ν = 1 → ν=2 transition in blue. The 2D IR spectrum is elongated along the diagonal, which is normally associated with an inhomogeneous distribution of vibrational frequencies. However, in this case, this diagonal width comes from the overlap between the two peaks identified in the FTIR. We fit the diagonal slice to a sum of two Gaussians, and compared to Figure 7.2: FTIR spectra of the Amide and Ester carbonyls stretches of DMPC with (red) and without Gramicidin D (black). The enlarged region shows the redshift on addition of Gramicidin D. The second derivatives of the spectra are shown in dotted lines, highlighting the underlying two-peak structure. The amide I region spans ≈1600 cm−1 to ≈1700 cm−1 , the ester region spans ≈1700 cm−1 to ≈1770 cm−1 .

the anti-diagonal widths to reveal two peaks with an ellipticity ratio of ≈1.01.25 and a homogeneous linewidth of Γ ≈13 cm−1 (FWHM). This demonstrates the system is better described as two distinct types of ester carbonyl

populations, each with a largely homogeneous lineshape. Here “homogeneous” means that the range of environments available to each population is sampled rapidly compared to the earliest waiting time of 150 fs in our experiment. Configurational dynamics of water’s hy-

7.4 FTIR and 2D IR Spectroscopy of Lipids

193

drogen bonds and closely related electric field fluctuations are one of the few motions that can occur this rapidly. Addition of gD results in an increase in the intensity of Peak A, which is readily apparent in the 2D IR difference spectra comparing bilayers with and without gD (Figure 7.3) This demonstrates that there is an increase in the number of hydrogen-bonded esters. However, the source of these additional hydrogen bonds is not clear from the spectra alone. It is possible that gD disrupts packing in the membrane, allowing greater solvent penetration. Simulations have shown, however, that Tryptophan sidechains of gD form hydrogen bonds to the lipid ester carbonyls [173], which may provide an alternative explanation. The change in intensity of peak A on addition of gD is notably more distinct in the 2D IR spectra than FTIR, which suggests the transition dipole moments of the two transitions are different; Both FTIR and 2D IR signals scale linearly with concentration, but as µ2 and µ4 , respectively, with dipole strength.

7.4.1

Dynamics Revealed by 2D IR Spectroscopy

The DMPC bilayer dynamics are encoded in 2D IR spectra primarily in time-dependent intensities and lineshapes. The 2D IR lineshapes in Figure 7.4 have contours that transition from diagonally elongated to round as a function of waiting time. These lineshape changes may be quantified by the rate of change in the center-line slope (CLS), which measures the loss of correlation between excitation and detection frequencies. In the case of a single peak, the CLS decay is directly proportional to the FFCF, and thereby the correlation function for fluctuating electric fields. Determining the FFCF is key to understanding the fluctuations of forces in the lipid membrane. Bilayers with and without gD show an exponential decay of the CLS. The gD containing sample decays significantly slower than the pure DMPC bilayer (6.2 ps vs. 2.9 ps), indicating gD is a significant spectroscopic perturbation. The presence of multiple peaks complicates the molecular interpretation of this relaxation process [71]. In order to separate the electric field fluctuations encoded by Peak A and Peak B, it is necessary to distinguish it from other processes which may influence the

194


Figure 7.3: 2D IR spectra of the DMPC bilayer ester carbonyls at T=150 fs with (upper right) and without (upper left) Gramicidin D. Lower right panel shows the difference spectrum (each spectrum normalized to diagonal area before subtraction), highlighting the increase in intensity of the lower frequency peak.


195

2D IR lineshape. These factors include different vibrational lifetimes (T1 ), crosspeaks, vibrational energy transfer (VET) between oscillators, and oscillators which exchange during the waiting time. Of these effects, crosspeaks from coupling between oscillators are not a significant contributor, since they are not prominent in our early waiting time spectra, although other studies have reported weak crosspeaks [174]. Integration of intensities over Peaks A and B as a function of waiting time reveals the vibrational lifetime of the ester carbonyl populations (Table 7.1). The lifetime of the Peak B is relatively insensitive to the presence of gD, changing from T1 =1.7 ps to T1 =1.8 ps on addition of gD to the DMPC bilayer. The lifetime of Peak A, however, lengthens by 30% in the presence of gD from T1 =1.0 ps to T1 =1.3 ps. The lifetime of the Peak A in the pure DMPC bilayer is similar to the vibrational lifetime of small ester groups in free aqueous solution [175]. Although there is no simple relationship between T1 and electric field fluctuations, the increase in the lifetime of Peak A on addition of gD provides additional insight into the simultaneous increase in Peak A intensity on addition of gD. If the incorporation of gD into the bilayer simply resulted in a greater population of water exposed ester sites, we would not expect a change in T1 unless this was accompanied by a significant change in the water’s dynamics. It is possible that gD indirectly slows the H-bond configurational fluctuations of water associated with lipid esters. Also possible is a change in the relaxation mechanism of Peak A, present only when gD is in the sample. This is suggestive of protein sidechain-bound esters, as discussed earlier. Since both processes may contribute, the decay is discussed in terms of effective timescale, though it may not relate directly to a single physical process. Two final dynamical processes which may contribute to 2D IR spectra are VET between esters and the making and breaking of hydrogen bonds. Both effects will appear as a crosspeak which grows relative to the diagonal peak. Indications of this are seen in spectra with T> 3 ps, though it is not possible from the spectrum alone to determine which process is responsible. Replacing 5% of the DMPC lipids with sphingomyelin introduces a spectrally distinct amide carbonyl stretch vibration (≈ 1660 cm−1 ) to the system, which can be used

196


Figure 7.4: Selected delays from the waiting time series of the pure DMPC bilayer. The Center-Line derived from vertical slices is highlighted by the white line, showing the change in slope from near diagonal at T=150 fs to near-horizontal at 5 ps. The CLS decays for the DMPC sample with and without Gramicidin D are shown in blue and red, respectively, bottom right. The full waiting time series for both samples are presented in §7.A.

7.4 FTIR and 2D IR Spectroscopy of Lipids Sample

DMPC

DMPC

197

Peak A Vibrational

Peak A Vibrational

CLS Decay

Lifetime (T1 )

Lifetime (T1 )

(τCLS )

1 ps

1.7 ps

2.9 ps

(0.9, 1.2 ps)

(1.4, 2.1 ps)

(2.4, 3.7 ps)

1.3 ps

1.8 ps

6.2 ps

(1.2, 1.6 ps)

(1.4, 2.4 ps)

(4.9, 8.4 ps)

Table 7.1: Experimentally obtained timescales for the DMPC bilayers. The vibrational lifetime measures the decay of intensity in Peaks A and B, obtained by integrating ω1 and ω3 over the windows 1710-1730cm−1 and 1735-1770cm−1 , respectively. Confidence intervals (95%) are given in parenthesis. to test for VET. Figure 7.5 shows the 2D IR spectra of the sphingomyelin-doped bilayers. At T=150 fs, a pair of amide resonances at 1648 and 1665 cm−1 also indicate distinct carbonyl environments for sphingomyelin with and without a hydrogen bond to the amide group. A crosspeak is clearly visible for the longest waiting times, demonstrating that VET occurs between different lipid groups, although considerably slower than the vibrational frequency fluctuations. This is in contrast to intermolecular crosspeaks between amide and ester at short T observed for peptides associated with or folded within lipids [157,176]. The existence of VET pathways does not exclude hydrogen bond switching as a factor. However, experiments on small esters in methanol have estimated hydrogen bond switching timescales to be ≈18 ps, which is significantly greater than the time window of our experiments [177].

7.4.2

A Heuristic Spectral Model for Lipid Bilayers

To extract quantitative dynamical information about both Peaks A and B, the timeevolution of the experimental 2D IR spectra is modeled using an established model for two fluctuating coupled oscillators [59]. For each oscillator, the transition frequency, the dipole strength, the T1 , and the FFCF must be determined. Additionally, the oscillators may transfer energy between each other. Although this results in many parameters, several can be constrained to other experimental data. The parameters used in the model are summarized in Table 7.2, and a fuller discussion of the origins of all parameters is given in

198


the §7.C.4. The primary objective of the model is to determine the FFCF parameters, which are key to understanding the electric field dynamics of the system. The primary metrics for successful reproduction of the spectra are that the transition of the lineshape from diagonally elongated to rounded, but asymmetric in ω1 is captured; the time-dependent frequency shifts are reproduced; and the effect of adding gD on the CLS decay can be explained. The vibrational relaxation dynamics of N-methylacetamide (NMA), a well-characterized single amide I oscillator, are used as a starting point. A biexponential decay of the FFCF with a 56 fs timescale and a second timescale of ≈1 ps were observed in NMA in D2 O [178]. These processes are commonly interpreted as arising from water librational motions and hydrogen bond reorganization, respectively. A biexponential form was assumed for all FFCFs, and the timescales varied to optimize agreement between experiment and simulation. To incorporate the effects of VET and/or chemical exchange, a crosspeak with time-dependent intensity was added, as outlined in Reference [59]. The model and experimental data are compared in Figure 7.6. To obtain satisfactory agreement with experiment, two changes to the biexponential form of the NMA FFCF were needed (Table 7.2). It was necessary to weight the amplitude of the fast term of the FFCF more heavily than the slower decay, and decrease the timescale for the long relaxation process to 4 ps for Peak B, which is the long timescale of the FFCF of NMA in CDCl3 [178]. This gave good agreement between the simulated and experimental CLS decay (simulated 2.6 ps, experimental 2.9 ps). This modification of the model is consistent with the earlier discussion of the two peaks in the experimental spectra in terms of hydrogen bonded and not hydrogen bonded ester groups. The model may be extended to include the influence of gD by adding a third oscillator into the system. By simply shifting the peak frequency of the third oscillator to lie between the two initial peaks (justified by the simulations in §7.4.3), and adjusting the vibrational lifetime and dipole strength to match experiment, it is possible to replicate the slowdown of the CLS decay observed experimentally. These parameters are summarized in Table 7.2.


199

Experimental values for the vibrational Stark shift of carbonyls range from 0.5-1.5 cm−1 /(MV/cm),

which

combined with the FFCF amplitudes provides an estimate the standard deviation of the electric field fluctuations to be ≈10 MV/cm. These results are in excellent agreement with the vibrational Stark shift studies using unnatural amino acids as local probes [153, 179]. Typical values of electric fields across biological lipid membranes are < 1 MV/cm, which is significantly less than the fields determined here (though field focusing effects may reduce this disparity [180]), highlighting the importance of understanding the dynamics of these local electric field fluctuations. Figure 7.5: 2D IR spectra of a 20:1 Though this model is phenomenoDMPC:Sphingomyelin bilayer at early and late waiting times. Crosspeak is highlighted in the bot- logical, it reveals a number of details tom panel. about the DMPC ester vibrations. Oscillators in both environments are dominated by fast ( 5 ps. Similar to the experimentally derived FFCF, all of the simulated correlation functions have a bimodal decay with a prominent fast decay on sub-100 fs time-scales followed by a longer, smaller amplitude decay. For both the hydrogen bonded and non-hydrogen bonded cases the FFCF has decayed to < 1/3 of the initial value within 150 fs, the earliest waiting time which can be confidently measured in these experiments. The rapid decay of the FFCF is consistent with the largely-homogeneous lineshape for each peak. Both FFCF decays involving a hydrogen bond (from solvent or protein) show oscillations. The solvent-bound FFCF shows a beat at ≈ 150 fs (187 cm−1 ), consistent with FFCFs in other hydrogen bonded systems where it has been attributed to stretching motions of the O-O hydrogen bond with water [184]. The protein-bound FFCF shows a lower-frequency beat (≈ 300 fs, 93 cm−1 ), consistent with hydrogen bond configurational changes involving a more massive partner. The zero-time value of our FFCF is proportional to the variance of the electric field fluctuations for each peak: Cωω = c20 hδE 2 i. The standard deviation of the electric field fluctuations in the MD simulations varied from 12.9 MV/cm for the non-hydrogen bonded peak to 23.2 MV/cm for the hydrogen-bonded peak, values in good agreement with previous studies on the range of electrostatic environments in lipid membranes and with our heuristic model estimates. The rapid decay of the hydrogen bonded peak is consistent with the behavior of solvent exposed oscillators. Similar behavior was found in the CN− ion in aqueous solution [185], where the rapid FFCF decay was attributed to a rapidly varying electric field on the oscillator from fluctuations of the hydrogen bond. Though the large electric fields produced by


205

hydrogen bonding are very localized, the ubiquity of solute-water hydrogen bonding, even in membrane proteins [186], suggests the water-induced electrostatic fluctuations should not be neglected. A more surprising result, perhaps, is the rapid decay of the non-hydrogen bonded oscillator, which occurs on the same timescale as the fast component of the hydrogen bonded FFCF. In the absence of hydrogen bonding to solvent, the electrostatic fluctuations in the lipid bilayer will be most strongly influenced by the zwitterionic headgroup. The phosphate group, in particular, is found within only a few angstroms of the ester carbonyl groups.

Fluctua-

tions of this distance between this negatively charged group and the ester carbonyl provides a plausible mechanism for the rapid decay of the non-hydrogen bonded FFCF. There is only limited correlation in our model between the carbonyl frequency of the ester and its distance from the phosFigure 7.7: a) Probability distributions for the three ester environments, as a function of frequency. b) FFCF decays for the three ester environments in the system. Dotted lines are the values from the simulation, solid lines are fits to the functions given in Table 3. The grey region highlights the length of the instrument response function.

phate group located on the same molecule (¯ r = 0.168). However, the lipid is conformationally flexible enough that the charged group closest to the ester is not necessarily the charged group on that lipid molecule.

Thus charged groups from neighboring lipids may be responsible for the fluctuations experienced by the ester group; in this case, it is not clear what the appropriate coordinate to

206

Chapter 7. High-Amplitude Electric Field Fluctuations in Lipid Membranes Carbonyl Group

a1

τ1 / fs

a2

τ2 / fs

a3

τ3 / fs

ω / cm−1

Free Ester

0.79

70

0.21

972

/

/

/

Water-Ester

0.48

56

0.18

380

0.34

67

186

Trp-Ester

0.27

218

/

/

0.73

87

93

Table 7.3: FFCF parameters determiend from MD simulations. Fit to the functional form C(t) = a1 exp(−t/τ1 ) + a2 exp(−t/τ2 ) + (1 − a1 − a2 )exp(−t/τ3 )cos(2πωt) consider is. Whatever the origin of these fluctuations, it is clear that the lipids (and any intercalated membrane protein) are subject to a rapidly varying electrostatic environment even in the absence of local, specific interactions such as hydrogen bonds.

7.4.6

Influence of Gramicidin D

Tryptophan (Trp) has been observed in previous simulations to donate a hydrogen bond to the surrounding esters [173]. This is also observed in the simulations presented here, where 10% of the hydrogen bonds to lipid ester groups within 5˚ Aof the protein are from the Trp sidechains. As seen in Figure 7.7, the FFCF associated with carbonyl groups hydrogen bonded to the Trp sidechain yields a decay qualitatively similar to the solvent-carbonyl hydrogen bonds. The frequency of the beat, however, is much lower (93 cm−1 ), which is consistent with a hindered rotational motion involving a more massive partner. This shift in frequency may provide a route to studying protein-lipid interaction in the future. The frequency trajectories of Trp hydrogen bonded esters span a narrower range than their solvent hydrogen bonded counterparts, and are ≈ 5 cm−1 blueshifted. The fluctuations experienced by the Trp-bound esters are only ≈ 75% of the magnitude of the solvent-bound esters. This is of particular interest, since the hydrogen bonding will not only effect to properties of the carbonyl, but also the properties of the sidechain, such as dipole moment. Thus, these fluctuations between the protein and lipid will induce rapid fluctuations in the properties of the sidechain. These time-dependent changes in sidechain properties may have implications for the functional behavior of gD. Changes to the permanent electric dipole of the sidechains have


207

been shown to result in significant changes in conductance, though no major structural changes occur [27]. This has been attributed to the net dipole from the Trp sidechains playing a key role in the conduction of ions through the channel. Density Functional Theory (DFT) calculations (shown in §7.F) show that the dipole of the indole sidechain in Trpn changes by 0.96 D on the formation of a hydrogen bond. Fluctuations of the hydrogen bond acceptor (such as the ester groups in the lipid) will give rise to a time-dependent dipole. This effect has two implications. First, calculations which utilize the dipole moment of the sidechains should use a weighted average of dipoles arising from different hydrogen bonding configurations. Second the impact of the temporal behavior of the fluctuations should be considered; the mean transit time of an ion through a gD channel is 10-100 ns. Hydrogen bond fluctuations on the ps to ns timescale will present a rugged energy landscape for this ion during transit.

7.4.7

Previous 2D IR Studies of Lipid Bilayer Dynamics

The experiments in this chapter use 2D IR to probe the frequency fluctuations (and from this extract electric field and force fluctuations) of esters in lipid bilayer. However, other time-resolved infrared experiments on the solvation of phospholipids have a bearing on the interpretation of the experiments presented here. Several studies have probed the waterester interactions, finding significant vibrational interactions between ester carbonyls and water on a sub-ps timescale, which is consistent with strongly associated hydrogen-bonding partners [187–189]. Volkov et al concluded, using 2D IR spectroscopy in combination with chain-specific isotope labelling, that the inhomogeneity arises primarily from electrostatic inhomogeneity, not from differences in the local environments of the and sn-1 and sn-2 chains [174]. Significant efforts have also been made to understand the dynamics of membraneassociated water [190,191]. These studies paint a complex picture; the dynamics of membraneassociated waters are significantly different from bulk water. Simulations suggest these studies are primarily sensitive to waters around the charged headgroups, specifically the

208


phosphate groups. Of particular relevance is the slowdown in physical reorientation of water molecules adjacent to lipid membranes [190]. In the context of the results presented here, this would suggest that fluctuations in the system on a sub-ns timescale arise not from “structural” changes, but from subtle variations of strong interactions. Other efforts have studied the dynamics of other regions of the lipid bilayer. Costard et al [192] probed the phosphate headgroup dynamics in reverse micelles directly, and Kel et al [161] used a metal carbonyl probe to study the dynamics of the lipid alkyl chain. Lipid bilayers play host to a wide range of timescales. The phosphate groups, initially inhomogeneous, experience solvent fluctuations on a 300 fs timescale, while the W(CO)6 probe in the alkyl interior of the bilayer experiences fluctuations on a timescale of 10 ps or greater. Ester carbonyls can now be added to this picture. These carbonyls appear to experience fluctuations on a < 100 fs timescale, irrespective of whether the carbonyl groups participates in a hydrogen bond. One interesting implication, is that a transmembrane protein will experience a bath of fluctuating forces with depth-dependent timescales and magnitudes. This is consistent with the simulations of Mukherjee et al [157] who found a distinct positional dependence of the electrostatic fluctuations experienced by the CD3ζ protein, and found that a significant fraction of the electrostatic inhomogeneity experienced by amide carbonyls in the protein originated from interactions with the lipids.

7.4.8

Ultrafast Fluctuations as a Source of Spatially Inhomogeneous Friction

As discussed earlier, fluctuating forces in the system provide a source of friction in the system. The force-force correlation function is related to the friction coefficient by the R∞ fluctuation-dissipation theorem, 2kB T γm = −∞ hF (0)F (t)i. The friction coefficient depends then on both the magnitude and timescales of the fluctuations. This has significant implications for the friction term in lipid membranes. It is possible to calculate values for friction coefficient arising from the electrostatic fluctuations explored here. Values tabulated in Table 7.4 are the friction coefficient experi-

7.4 FTIR and 2D IR Spectroscopy of Lipids mass (amu)

friction γ (ps−1 ))

1

95.7

174 (Arg sidechain)

0.55

2245 (TM helix in VSD)

0.044

209

Table 7.4: Friction coefficients for various masses calculated using the experimentallydetermined FFCF. enced by a group of mass m, with a net charge of +1e. Masses corresponding to relevant to voltage gated processes (e.g. arginine side chain, 174 amu, transmembrane helix in VoltageSensing Domain (VSD) protein, 2245 amu) are used to highlight the importance of this friction term. These numbers, arising solely from the sub-nanosecond electrostatic fluctuations measure here, are the same order of magnitude as the friction coefficient used in protein folding simulations [193, 194]. The importance of internal friction in protein folding is well-established [12, 195]. The comparable magnitudes of friction coefficients in protein folding and determined experimentally in this work suggests the electrostatic friction may have a significant functional role. The forces probed in this study, originating from electric fields will decay in amplitude as r−2 , where r is the distance from the charged group. Thus the magnitude of these forces (and so the friction coefficient arising from them) should be significantly smaller in the hydrophobic region of the lipid bilayer, far from the charged headgroups. This is consistent with the work of Shrestha et al, showing a spatial dependence to the electric field experienced by a vibrational probe in lipid bilayers [153]. A second consideration is the existence of other sources of fluctuating forces. The study by Kel et al on the fluctuations in the hydrophobic interior of lipid membrane revealed a timescale of 10 ps [161]. The distance of this region from the charged headgroups makes it unlikely that these are responsible for the fluctuations probed in this study; a more likely interpretation is that a second set of forces is responsible (for example, structural changes from acyl chain rotational isomerization are known to occur on a picosecond timescale). This suggests it is very likely that the friction coefficient in lipid membranes is spatially

210


inhomogeneous, and that this inhomogeneity arises from differences in force fluctuations. Position dependent friction has been suggested before in the context of membrane protein function. The model developed by Sigg et al to explain the gating kinetics of the Shaker K+ channel (which involves the motion of charged sidechain groups across the membrane) described the process as diffusion across a 1D potential of mean force [145]. To accurately model the experimental µs kinetics, a position dependent friction was required. Our results, in conjunction with previous 2D IR studies of lipids outlined earlier, provide a conceptual basis for this. The heterogeneity of fluctuation timescales and magnitudes suggests this position-dependent friction is a result of the spatial heterogeneity of the dynamics in the membrane. Diffusion across some lower-dimensional potential is a powerful tool for describing the kinetics and mechanism of voltage gated processes. However, these results suggest special care must be taken in their construction and interaction.

7.5

Summary and Outlook

The results presented in this chapter demonstrate that there are significant electrostatic fluctuations on the femtosecond and picosecond timescales. These fluctuations are related to extremely high electric field amplitudes (on the order of 10 MV/cm). The fluctuations shape the underlying potential energy landscape for functional processes; in the language of diffusion across a potential energy surface, these fluctuations act as a source of friction. The results presented here may provide some of the first experimental evidence for the position dependence of the friction coefficient in biological systems, which has been predicted to have a significant impact on the functional behaviour of voltage-gated ion channels. A natural extension of the work presented here is to explore how these fluctuations are modulated by small molecules and different lipid compositions. Preliminary studies in this vein are explored later in §9. An orthogonal line of research would be to probe the position dependence of the fluctuations in the system; this would require the use of non-


211

native probes, such as the series of membrane spanning peptides used by Shrestha et al to study the position dependence of the electrostatic inhomogeneity. What is clear is that understanding the sub-nanosecond dynamics of biological systems is key to fully understanding their functional behaviors. Efforts in this direction have been hindered by a lack of available experimental data; 2D IR promises to offer new insight in this area.

212

7.A


2D IR Waiting Time Series

The full waiting time series for DMPC-containing samples ranges from a waiting time of 150 fs to 10 ps. The earliest waiting time is determined by the pulse length used (100 fs) because of non-resonant response from the sample windows, and ambiguities in the interpretation of signals when pulses overlap significantly. At the longest waiting time (≈ 8× T1 ), the signal has decayed by a factor of 3 × 103 from the initial value.

Figure 7.8: 2D IR waiting time series for pure DMPC bilayers

7.B 2D IR Waiting Time Series

213

One of the main limiting factors here is preparation of samples of sufficient optical quality that scattered light is negligible even at these low signal levels. The effect of this scattered light can be observed in the late time spectra in the form of lineshape distortions giving the induced absorption a diagonally elongated component.

Figure 7.9: 2D IR waiting time series for 20:1 DMPC:gD bilayers

214


7.B

Vibrational Energy Transfer in Lipid Bilayers

To determine whether VET is significant in lipid systems, a sample with 5% of the DMPC lipids replaced by sphingomyelin (SM) was prepared. SM has an amide group, rather than an ester, which shifts the absorption of the carbonyl to ≈ 1650 cm−1 . Growth of a crosspeak between the ester and amide groups can only results from VET. Comparison of the T=150 fs and T=5 ps spectra demonstrates the growth of a crosspeak between the ester and amide carbonyls. Although it would be desirable to use this system to estimate the rate of energy transfer between carbonyls in lipid bilayers, this was not possible in this case because vibrational relaxation to other modes is much more rapid than energy transfer between peaks. A sample with a great enough optical density to quantitatively detect the crosspeak signal at > 3 the vibrational lifetime would cause significant intensity distortions due to re-absorption of the nonlinear signal by the sample, and from saturation of the detector at early times. However, the presence of the crosspeak in Figure 7.5 is enough to demonstrate that VET must be considered in any model of the time-evolution of the spectrum. Energy transfer between adjacent carbonyls in the tri-Alanine peptide has been measured previously. This system is similar in some ways to the ester carbonyls in bilayers. The conformational freedom of such a small peptide means at any point in time, the sample contains a large range of angles between carbonyl groups. The distance between adjacent carbonyls (≈ 3.5˚ A) is similar to the peak in the radial distribution function of ester carbonyls in lipid bilayers determined from MD simulations (≈ 3˚ A). The observed energy transfer rates in Reference [67], 0.07 ps−1 , 0.14 ps−1 are consistent with the data here. VET and chemical exchange are not mutually exclusive; the demonstration here that VET is a factor does not rule out the possibility that chemical exchange also takes place. However, experiments with esters in methanol have found hydrogen bond making and breaking to occur on a 17-18ps timescale [177] – which is significantly longer than the window of observation.

7.C Heuristic Model for Esters

7.C

215

Heuristic Model for Esters

The heuristic model used to calculate the time evolution of the 2D IR spectra of DMPC is based on the model presented in Reference [59]. In this model, the system is described by two coupled oscillators, and the third order response functions are calculated for the Rephasing and Non-Rephasing pathways. Each response function has a lineshape function associated with it (see Reference [61] for the complete expressions) which depends on the pathway and the FFCF of the oscillators. The ester groups are considered to be only weakly coupled to each other, so the eigenstates of the system are well-approximated by the site basis (hydrogen-bonded and non-hydrogen bonded).

7.C.1

Response Functions

For the two-oscillator system presented, there are fourteen vibrational response functions (and their complex conjugates) to be considered (ignoring pathways involving zero- and two-quantum transitions). Each response function takes the form of a product of complex exponentials, and a real-valued decaying dephasing function. The oscillation of the exponential encodes the frequency of the transitions, while the dephasing function will determine the lineshape of the function. Full expressions are given in §2.

7.C.2

Orientational Factors

Unlike Reference [59], where the six-level system describes two oscillators with on the same molecule, here the two oscillators describe two different environments in our system. To determine the angle to use in the orientational factor, we first make the assumption that all our oscillators lie in a plane (which our simulations suggest is a reasonable approximation), and are randomly distributed in this plane. From this, the expectation value of cos2 (θ) is calculated to determine our orientational factors. If we have, P (θ) = 1/2π, then hcos2 (θ)i = 1/2, which corresponds to hθi = 45◦ . This is the angle used in the simulations.

216


7.C.3

Lineshape functions

By analogy between the ester carbonyl groups and NMA, a biexponential functional form for the FFCF is assumed. This is confirmed by the MD simulations. Thus

C(t) =

∆21 exp

t − τ1

+

∆22 exp

t − . τ2

(7.1)

This leads to the Kubo lineshape function for each oscillator:

g(t) =

X i

t t ∆2i τi2 exp − + −1 τi τi

(7.2)

which may then be used with the dephasing function in Reference [61] to calculate the lineshape of each transition. One modification to the expression in Reference [61] is the inclusion of an exponential decay to incorporate the effects of lifetime broadening in the system. This was necessary to achieve satisfactory linewidths with realistic parameters for the FFCF.

7.C.4

Parameter Set

The model described about requires 9 parameters for each peak - the transition frequency (ω), the anharmonicity (δ), the relative dipole strength (µ/µ0 ), FFCF timescales (τi ), and magnitudes (∆i ), VET rate to other groups (τV ET ) and vibrational lifetime (T1 ). The situation is not as dire as it initially seems however, as several of these parameters may be constrained from the experimental data. The transition frequency may be estimated from the FTIR data. FTIR data is used because it is higher resolution; the 2D IR data is collected using a 64 pixel array, which gives ≈ 4 − 6 cm−1 resolution, whereas the FTIR data are collected at 1 cm−1 resolution. The relative dipole strengths may be estimated by comparing the intensities of the peaks in the FTIR and 2D IR spectra – both scale linearly with concentration, but FTIR absorption is the square of the transition dipole moment of a vibration, while 2D IR is the fourth power

7.D Ester Frequency Map

217

of the transition dipole moment. Vibrational lifetimes may be determined by integrating spectral regions of the 2D IR spectra at each waiting time and fitting the decay of this integrated intensity as a function of time. Volkov et al [174] report a value of 20 cm−1 for the anharmonicity for DMPC ester groups, which is adopted here. The average value from Reference [67] for energy transfer between neighboring Amide groups as our estimate for VET between ester groups. No crosspeaks due to coupling between the oscillators were included because they were not apparent in the earliest waiting time spectra, and would not be expected to contribute significantly to this lineshape analysis. The remaining parameters to be determined relate to the FFCF, and are outlined in §7.4.2. Additionally, to replicate the experimental lineshape, it was necessary to include a small amount of inhomogeneity in each peak was necessary (2 cm−1 ). This likely represents the fact that our assumption that the environments of the ester carbonyls can be described by only two oscillators underestimates the inhomogeneity within each environment. However, the small value of this correction suggests this approximation is actually rather good.

7.D

Ester Frequency Map

The simplest functional form for a frequency map is a linear dependence on some electrostatic variable. In previous work, evaluation of the projection of the electric field onto the carbonyl bond axis was found to produce excellent agreement between simulation and experiment. In this case, the site energy of each oscillator is given by ω = ω0 + co Exo , and so has two parameters to optimize - the zero field frequency ω0 and the Stark coefficient co . The metrics for success of the map were chosen to be the center frequency of each peak, and the FWHM of each peak, rather than a direct comparison between the experimental

218


Figure 7.10: Simulated 2D IR waiting time series for pure DMPC bilayers

7.D Ester Frequency Map

Figure 7.11: Simulated 2D IR waiting time series for 20:1 DMPC:gD bilayers

219

220


Figure 7.12: Simulated center frequencies (a) and FWHM (b) for ester peaks as a function of Stark tuning coefficient. Peak A is shown in red, Peak B is shown in black. Dashed lines show experimental values. and simulated spectra because the presence of hydrogen-bonds can result in non-Condon effects which may alter the dipole strength in non-trivial ways. It is found that with this simple map it is not possible to replicate the frequencies and FWHM - only one metric can be optimized. There are several possibilities for this inability to capture the linewidth of the peaks. One is that more sites need to be included to sample the full inhomogeneity of the electrostatic environment. This would require significantly more data points to test, since there are only have four parameters to test with this data (two linewidths and two center frequencies). Anything more complex than a two-site model would be under-determined. More data points could be obtained from chemical modification of the lipid - for example, changing the zwitterionic DMPC to anionic DMPG would change the local electrostatics of the system, and so provide four new metrics to optimize against. Another intriguing possibility is that the fields in the MD simulation themselves are too large. The stark tuning rate from our map is ≈2.5 times smaller than experimental values, which is the factor by which it has been suggested MD simulations overestimate the effect of electrostatic fields. This would result in steeper electrostatic potentials, which would limit the range of conformations available at kB T . One question that follows naturally from this is whether polarizable forcefields would negate this electrostatic overestimation.

7.F Diagonal and Anti-diagonal Linewidths

7.E

221

Diagonal and Anti-diagonal Linewidths

A slice along the diagonal was fit to the sum of two Gaussian function,

S=

X

ai exp

2

− (ω − ωi ) /σ

2

to determine the diagonal linewidths. At the center frequencies of each of these peaks, an anti-diagonal slice was fit to two oppositely-signed Gaussian functions to determine the homogeneous linewidth. Parameters for the fits are given in Table 7.5. The fits reveal that the ellipticity of each peak ranges between 1 and 1.3, demonstrating that each peak is relatively homogeneous. This is consistent with the simulations which suggest there is a rapid < 150 fs component to the FFCF. a1

ω1 / cm−1

σ1 / cm−1

a2

ω2 / cm−1

σ2 / cm−1

Diagonal Slice

1

1728

10.99

0.39

1741

8.23

Upper Anti-diagonal Slice

1

-2

7.73

-0.68

-6

6.85

Lower Anti-diagonal Slice

1

-2

7.99

-0.98

-6

8.9

Table 7.5: Fit parameters for the diagonal and anti-diagonal slices.

7.F

Tryptophan Dipole Moment

To determine the effect of hydrogen bonding on the permanent electric dipole moment of Tryptophan, a series of DFT calculations were performed on the indole sidechain of Tryptophan and a single water molecule. The calculations were performed using Gaussian 09 with the B3LYP hybrid functional using the 6-31(d,p) basis [196]. The structure was first optimized, which resulted in the water oxygen residing ≈3 ˚ Afrom the indole nitrogen. This water was then shifted further away from the indole in 1 ˚ Asteps, along the N-O axis, and the dipole moment at each step calculated. The calculated dipole moment decays exponentially with water-indole separation, and the limiting value of ≈2.4 ˚ Ais in good agreement with

222


Figure 7.13: Fits to diagonal and anti-diagonal slices of the DMPC 2D IR spectrum. Each peak is fit to the sum of two Gaussian functions. Individual components of the fit are shown, and the sum is shown in red. Colors correspond to the colored regions in the contour plot (top left).

7.6 Tryptophan Dipole Moment

223

Figure 7.14: Simulated 2D IR waiting time series for pure DMPC bilayers. Inset shows the structure used in the calculations. dN O is the distance between the indole nitrogen and water oxygen. previous calculations of the isolated indole group. This is shown in Figure 7.14. Over the course of a few ˚ A displacement, the effect of the water becomes negligible.

224


Chapter 8

Conformational Changes in an Ion Channel Driven by a Membrane Phase Transition 8.1

Introduction

Proteins are remarkable molecular machines, capable of performing complex chemical tasks with both speed and precision. Proteins do not operate in isolation, however; their function is tightly coupled to the environment around them. Changes in pH, temperature or ionic strength modulate the function of proteins. This environmental sensitivity is even more pronounced in membrane proteins, where the variations in the membrane environment can significantly alter the functional behavior of proteins. Many examples of proteins responding to changes in the surrounding membrane environment are known. These range from mechanosensitive ion channels (sensitive to changes in lateral pressure) [197], to G-protein-coupled receptors (membrane fluidity) [143], and even the activity of membrane-associated enzymes (membrane thickness) [18]. Though the functional effects are well-established, the mechanism of action is not. Developing this picture is challenging due the many length- and timescales over which membrane-protein interactions 225

Chapter 8. Conformational Changes in an Ion Channel Driven by a Membrane Phase 226 Transition operate. A yet more complicated picture arises when the membrane is viewed not as a static scaffold for the protein, but as a dynamic, fluctuating environment, which will exert timevarying forces on integral membrane proteins. Fluctuations in membranes have been shown to span timescales ranging from ps to >ms, and operate over length scales from 100 over the 6.5 ◦ C range of final temperatures sampled; the variation of τf ast with temperature is weaker, but still significant - rates change by a factor of 10 over the range of temperatures studied. These strong temperature dependences are strongly suggestive of critical slowing (as observed in DMPC vesicles and gramicidin D-containing DMPC vesicles). Though a firm assignment of critical slowing cannot be made on the basis of a sampled temperature range of only 6.5 ◦ C, with no low-temperature points, it is notable that critical slowing has been observed in previous T-jump experiments on cholesterol-containing membranes [212].

274

Chapter 9. Modulation of Multiscale Dynamics in Membranes by Cholesterol

Figure 9.8: Transient response of esters in cholesterol-containing lipids (top left). Representative kinetic trace for the loss of peak B (top right). Temperature-dependence of rates from T-jump experiments (bottom right). Schematic of the changes associated with lipid phase transition (bottom left).

9.5 Sub-Millisecond Dynamics of Cholesterolic Membranes

9.5.1

275

Discussion

Determination of the origins of the two timescales observed is key to understanding the effect cholesterol has on the kinetics of the gel-to-fluid phase transition. The phase diagram of cholesterol-containing membranes has been the subject of many studies, utilizing a variety of techniques. These studies suggest that, under the experimental conditions here, a third phase exists - the liquid-ordered (Lo ) phase, which is characterized by short-range orientational order, but long-range translational disorder [237]. The consensus in the literature is that the gel and Lo phase can coexist, but at the cholesterol concentrations used here, the high-temperature phase should be completely fluid [237]. Multiple initial states (gel and Lo ) which undergo phase transitions provides a plausible explanation for the observation of multiple timescales in the transient data. Even if the coexistence of multiple phases at T < Tm is the origin of the two timescales observed in the experimental data, this does not clarify which timescale corresponds to the gel→fluid transition, and which corresponds to the Lo →fluid transition. Progress may be made, however, by considering whether the gel→fluid transition observed in pure DMPC vesicles would be expected to speed up and slow down in the presence of cholesterol. A simple kinetic picture might argue that the rigid nature of cholesterol (which §9.4 established has effects on even the picosecond dynamics) might introduce a steric barrier to the transition, slowing the transition. However, this Arrhenius-like description does not fully capture the nature of the DMPC phase transition. Rather, the well-established critical nature of the phase transition in DMPC should be considered. The observed slowdown in the phase transition rate as T → Tm arises because the coherence length of fluctuations diverges; at Tm , the coherence length of these fluctuations is “infinite” [242]. This is obviously unphysical in a real, finite, system - though the idea that growth of the coherence length is responsible for the slowdown in rate holds. This is really another way of stating that the transition is cooperative - many lipids undergo the transition simultaneously. Addition of cholesterol, which results in the formation of gel and Lo domains, will limit the maximum coherence length - thus the gel→fluid transition would be expected to proceed faster than

276


Figure 9.9: Proposed mechanism for the change in T-jump rates on the addition of cholesterol. Gel and Lo domains are formed, which further limits the maximum coherence length of fluctuations in the system, resulting in a net increase in rate for the gel→fluid transition.

in pure DMPC vesicles. On this basis, τf ast is assigned as the gel→fluid transition, which leaves τslow as the Lo →fluid transition. This is summarized schematically in Figure 9.9. This formation of small domains in cholesterol-containing mixtures has suggested previously suggested from NMR studies [243] The question remains then, as to why cholesterol slows down the Lo →fluid, relative to the gel→fluid transition. It may be that, on average, the Lo domains are larger than the gel domains, resulting in greater cooperativity. Another explanation may be that the events required for melting (as observed by the ester hydration) are slower in the Lo phase. §8 describes the events leading up to this change in hydration - one of the key events in this process is the conformational change of the acyl tails. These conformational changes are believed to take place on the sub-ns timescale [244]. 2D IR measurements have demonstrated that fluctuations in this acyl region slow down upon the addition of cholesterol [161]. If this trend holds for the conformational isomerization, then this may provide an explanation for the slower rate for the Lo →fluid transition. These factors are not mutually exclusive -

9.6 Action of Cholesterol on Different Length- and Timescales

277

both may be at work. This discussion is consistent with T-jump turbidity experiments. Turbidity measurements are sensitive to structural changes which result in a change in the scattering crosssection of the lipid vesicles, and so are not sensitive to individual gel vs Lo phases. However, on addition of cholesterol, these overall shape changes show an increase in rate, and a decrease in how sharply peaked the critical slowing is around Tm . These observations are consistent with a decrease in the coherence length due to domain formation. Domain formation on the addition of cholesterol has been suggested previously. NMR [236] and simulations [245] have suggested that the extent of these domains is only hundreds of lipids, consistent with the data shown here.

9.6

Action of Cholesterol on Different Length- and Timescales

Cholesterol is one of the simplest naturally-occurring molecules in eukaryotes, yet it is known to introduce complex behavior into membranes. The experiments presented in this chapter reveal yet more complexity - the effect of cholesterol can vary not only with lipid phase, but also with the length-scale being examined. The picosecond electrostatic fluctuations of DMPC lipids becomes significantly slower in the presence of cholesterol, and increase in magnitude. This can be rationalized by considering the condensing (decrease in area-per-headgroup) effect of cholesterol, as well as the molecular rigidity of cholesterol arising from its multiple ring structures. These effects should be considered in light of the discussion in §7 - electrostatic fluctuations give rise to an effective friction which slows the transit of charged groups through the membrane. This friction is proportional to the timescales of the fluctuation, and the variance of the fluctuations. Both of these factors increase in the DMPC bilayers on the addition of cholesterol; the electrostatic friction is estimated to increase by a factor of 8 from the experimentally-determined parameters. This increase in friction would suggest that voltagegating mechanisms which rely on the movement of charged residues through the membrane

278


should become significantly slower in the presence of cholesterol. This is observed in some ion channels, such as the Kv 1.3 channel [246]. However, the very premise of this chapter (that cholesterol has multiple modes of action) complicates the interpretation of these experiments - cholesterol alters many properties of the membrane, such as lateral diffusion, bending rigidity and lipid packing, which may influence kinetics. This is likely responsible observations that in some cases, the presence of cholesterol can speed up the kinetics associated with voltage gating [20]. The effect of cholesterol on the phase transition of DMPC vesicles is more complicated. This single timescale observed in the pure vesicles becomes two when cholesterol is added one faster, one slower. The faster timescale arises from a confinement effect - the presence of two phases (gel and Lo ) restricts the coherence length of the fluctuations at the phase transition. The slower timescale arises from the transition of the cholesterol-specific Ld phase to the fluid phase; the rigid cholesterol acts to slow down the structural rearrangements in the bilayer. This suggests that cholesterol introduces spatial, and dynamic, inhomogeneity into the membrane. Interestingly, NMR experiments have reported observations consistent with rapid exchange of lipids between the gel and Lo phases [236], suggesting that these lipid domains are themselves highly dynamic. For these domains to have an effect on the coherence length of fluctuations in the vesicles, they must be small - the domain size may be estimated to be