12276
J. Phys. Chem. A 2009, 113, 12276–12284
Relativistic Interactions in the Radical Pair Model of Magnetic Field Sense in CRY-1 Protein of Arabidopsis thaliana Artur F. Izmaylov* and John C. Tully Department of Chemistry, Yale UniVersity, New HaVen, Connecticut 06520
Michael J. Frisch Gaussian, Inc., Wallingford, Connecticut 06492 ReceiVed: January 13, 2009; ReVised Manuscript ReceiVed: August 20, 2009
Experimentally, it has been shown that magnetic field sensitivity in living organisms is connected to the presence of blue-light photoreceptor cryptochromes. Cryptochromes transduce a light signal through a chain of chemical reactions involving the formation of intermediate biradicals. It was proposed that an external magnetic field affects the interconversion between singlet and triplet states of biradicals and thus interferes with the signal transduction chain. Theoretical modeling of this process requires an accurate evaluation of all interactions important for singlet-triplet interconversion: electron-electron, spin-orbit, spin-spin, hyperfine, and Zeeman. In the current study we investigate these interactions at the CIS level of theory applied to representative fragments of the CRY-1 protein in the plant Arabidopsis thaliana. We find, in contrast to previous simplified modeling (O. Efimova, O.; Hore, P. J. Biophys. J. 2008, 94, 1565), that the spin-spin interaction is significantly larger than the “exchange” interaction. Thus it is not canceled by the latter but rather dies off with the inter-radical separation. Also, we find that the spin-orbit interaction can play a significant role in singlet-triplet interconversion for short inter-radical distances, and the hyperfine interaction becomes the only coupling interaction for long inter-radical distances. I. Introduction It is quite fascinating that some living organisms can sense even small changes in the surrounding magnetic field: migratory birds can orient themselves based on the inclination of Earth’s magnetic field,1,2 fruitflies (Drosophila melanogaster) can be trained to use a magnetic field to find food,3 and even the plant Arabidopsis thaliana varies its hypocotyl growth rates due to magnetic field changes.4 In the two latter cases it is still unclear why nature supplied these species with a magnetic sensing ability. Experimental studies across different species reveal several common features and similar conditions necessary for magnetic field sensing: the blue part of the light spectrum (470 nm or 2.67 eV) must be available, only the magnetic field inclination and intensity but not the polarity matter, and cryptochrome protein knockout mutants lose magnetic sensing ability. These features provided solid evidence for assigning the key role in this phenomenon to the blue-light photoreceptor proteins cryptochromes.2-5 As implied by their name, cryptochromes’ functional role had been a mystery for some time and was unraveled only two decades ago.6 Currently, it is known that besides magnetic field sensing they regulate growth and development in plants and serve as circadian clocks in plants and animals. Another interesting feature of these proteins is their structural similarity with photolyases, enzymes responsible for repairing DNA damage from ultraviolet light.6 It has also been found that some steps in the cryptochrome photoactivation mechanism are analogous to those in the photolyase repair mechanism.6,7 Recently, the X-ray structure of the photolyase homology domain of the Arabidopsis cryptochrome CRY-1 has become available.8 This protein is deemed to be responsible for magnetic sensitivity in Arabidopsis thaliana.4 Also, as shown by Solov’yov
et al.9 the DNA sequences of the cryptochrome genes in the migratory bird Europian robin and the plant Arabidopsis thaliana are very similar. Thus studying the CRY-1 protein is an important step to understanding magnetic sensing in general, and it will be a primary focus of the current work. The basis for understanding the mechanism of the magnetic field influence is a radical pair model which is quite common for various magnetic field regulated reactions.10-12 In this model, the essential part of the mechanism is the formation of a biradical intermediate that leads to different product channels based upon its spin-symmetry. An external magnetic field in the presence of other magnetic interactions (e.g., hyperfine) can affect the singlet-triplet states ratio and thus alter the yield of a specific product.10 Although some important details of the CRY-1 photoactivation mechanism are still unclear, the main steps are as follows. First, the pterin-based chromophore 5,10-methenyltetrahydrofolate absorbs blue light and transmits energy to the second flavin-based chromophore flavin adenine dinucleotide (FAD). After the indirect excitation, FAD is protonated either from a neighboring aspartic acid (Asp396)7 or from a trapped water molecule13 (see Scheme 1, FADH+ mechanism). The excited FADH+ cation accepts an electron from a neighboring tryptophan (Trp400) and forms the radical pair of the neutral FADH• and Trp400•+ cation. Similar to photolyases, the CRY-1 protein has a chain of three tryptophan fragments near the FAD chromophore: Trp400, Trp377, and Trp324. Therefore, further steps involve electron hops from the next tryptophan fragment (Trp377) to the Trp400•+ cation forming FADH•-Trp377•+ biradical, and from Trp324 to the Trp377•+ cation forming the FADH•-Trp324•+ biradical. The deprotonation of the Trp324•+ cation constitutes the final step in a crucial for magnetic sensing part of the CRY-1
10.1021/jp900357f CCC: $40.75 2009 American Chemical Society Published on Web 10/06/2009
Role of Protein in Magnetic Field Sensing
J. Phys. Chem. A, Vol. 113, No. 44, 2009 12277
SCHEME 1:
Figure 1. Double-well potential with a finite probability of tunneling. E1 and E2 (solid lines) are energies of the first two states, Haa and Hbb (dashed lines) are ground state energies of individual wells as if there were no coupling.
photoactivation mechanism. In order to incorporate the magnetic field influence into the considered mechanism, it is assumed that the back recombination of each biradical is possible only from the singlet state. Hence, if the magnetic field can alter the singlet-triplet ratio it also can affect the rates of the back reactions and facilitate signal transduction by increasing the yield of the final deprotonated product. This description of the photoactivation mechanism closely follows the theoretical interpretation of experimental works by Solov’yov et al.9 and Efimova et al.14 However, experimental results in refs 7, 15, and 16 as well as arguments presented by Kao et al.13 for cryptochrome mechanisms in insects can lead to a different mechanism where electron hops precede the protonation of the FAD fragment (see Scheme 1, FAD•- mechanism). Both mechanisms are plausible, and we are not aware of experimental results that clearly favor one or the other. For instance, a recent study by Bouly et al.16 demonstrated that FADH• is the signaling state, but the study was not detailed enough to attribute the FADH• origin to FAD•- or FADH+. From a speculative perspective in both cases the first process facilitates the second: electron hops are energetically easier with the FADH+ form, and the anionic form FAD•- has a greater proton affinity than the initial FAD fragment. Structural similarities between photolyases and cryptochromes allow us to assign approximate time scale of 10-100 ns for the CRY-1 primary events depicted in Scheme 1; this estimate has been obtained for the Escherichia coli photolyase by Byrdin et al.17 Although electron hops are usually faster events than proton transports, in this time scale both events are probable. In order to postpone choosing between these two paths until we have more detailed knowledge, in the current work we model both pathways. In the system with two unpaired electrons that have equal g-factors, a magnetic field does not couple singlet and triplet states but only lifts the degeneracy of the triplet components. In order for an external magnetic field to influence the singlet-triplet transition, the system must have interactions which can either create a significant difference in g-factors for unpaired electrons or otherwise introduce a singlet-triplet coupling. The former factor is negligible considering the magnitude of the magnetic fields involved (e5 G).4 Thus, the key role is played by electron spins interacting with magnetic fields generated internally in spin-orbit, spin-spin, and nuclear hyperfine interactions.18 It is instructive to split all interactions between singlet and triplet states into two groups: (1) ones which contribute to the singlet-triplet splitting without mixing states of different spin multiplicity (e.g., spin-spin and “exchange”),
and (2) ones which couple states of different spin multiplicity (e.g., spin-orbit and hyperfine). The interplay between these two groups can be illustrated by mapping the singlet-triplet interconversion to single particle dynamics on the one-dimensional potential sketched in Figure 1. This mapping associates the first group of interactions to the difference between Haa and Hbb (energy levels for uncoupled wells), while the coupling elements Hab ) Hba* are related to the second group of interactions. If a particle is initially placed in well a, the wave function of the system at time t can be represented as a linear combination of stationary states of each well with time-dependent amplitudes Ca and Cb. Using a straightforward derivation,19 the individual probabilities to find a particle in each well are
|Ca|2(t) ) 1 -
|Cb|2(t) )
1 - cos2(ωt) r2 + 1
(1)
sin2(ωt) r2 + 1
(2)
where
r)
(Haa - Hbb) 2Hab
(3)
ω)
2Hab 2 √r + 1 p
(4)
Here, the direct analogy of the singlet-triplet ratio is |Ca|2(t)/|Cb|2(t). The magnetic field effect would be equivalent to adding a small energy contribution to the difference Haa - Hbb and corresponding change in the r parameter
r)
(Haa - Hbb + ) 2Hab
(5)
Considering the infinitesimal field limit ( f 0) the amplitude of both probabilities (eqs 1 and 2) is 2
[
Hab 1 2 f 12 2 H r +1 (Haa - Hbb) aa - Hbb
]
(6)
12278
J. Phys. Chem. A, Vol. 113, No. 44, 2009
Izmaylov et al.
Figure 2. Model systems of the flavin adenine dinucleotide (FAD) and tryptophan (Trp) molecules. The FADH model is obtained by adding the hydrogen atom to atom 8N.
Thus, the effect of the Zeeman interaction () becomes smaller for larger zero field splittings Haa - Hbb, and thus the spin-spin and “exchange” interactions should not be more than an order of magnitude larger than the Zeeman interaction if we are to observe a non-negligible magnetic field effect. The spin-orbit and hyperfine interactions associated with the coupling element Hab have a direct impact on the amplitude as well, but their contributions can also be damped by the splitting Haa - Hbb or, in other words, by the spin-spin and “exchange” interactions. This simple model vividly illustrates the importance of an accurate evaluation of both types of interactions to obtain an adequate description of the singlet-triplet interconversion dynamics. In previous work, some of these interactions were considered9,14 but in a simplified manner and not systematically. Also, the spin-orbit interaction that potentially plays a crucial role in the singlet-triplet interconversion was neglected. Therefore, this work extends upon previous theoretical studies with a careful analysis of the main relativistic interactions at the same level of theory which provides grounds for subsequent dynamics simulations. II. Theory A. Models. We model active ingredients of the radical pair mechanism by the fragments depicted in Figure 2. The model of the tryptophan (Trp) fragment has been successfully used in previous studies by Himo et al.20 and Lendzian et al.21 to model hyperfine coupling constants of Trp residues in various biological systems. Trimming the FAD molecule is motivated by its orientation in the crystal structure8 where the flavin ring system is the closest to the tryptophan chain part, and it is natural to expect that it plays the key role in the electron-hopping mechanism. Although our fragments are different from the real FAD and Trp molecules, we will refer to them as FAD and Trp for simplicity. For both mechanisms we consider two- and four-fragment models. The two-fragment model consists of three FAD(H+)-Trpi pairs, where i ) 400, 377, and 324. In the fourfragment model, the entire FAD(H+)-Trp400-Trp377-Trp324 chain is considered. All model fragments are optimized at the B3LYP/6-31G* level of theory and placed in the crystal structure orientation with respect to each other.8 Clearly, the four-fragment model is more computationally demanding, and therefore it has been studied with the moderate 6-31G* basis set. The two-fragment model has also been explored with the larger EPR-II basis set22 which is better suited for calculations of magnetic properties.22,23 All calculations described in this paper are done with a development version of the Gaussian program.24
B. Method. As one of the computationally efficient methods that can describe spin-symmetry, low-lying excited states, and charge transfer (CT) excitations, at least on a qualitative level, the configuration interaction with single excitations (CIS) method has been chosen. Another nice feature of the CIS approach is its straightforward atomic orbital (AO) direct formulation which allows one to treat very large systems.25 The “exchange” interaction can be directly extracted from CIS calculations as an energy difference between singlet and triplet CT states. In order to distinguish the CT solutions from other excited states, we employed the Mulliken charge analysis of excited state densities. There are three main relativistic interactions that depend explicitly on the electron spin variables: spin-spin (SS), spin-orbit (SO), and hyperfine (HF). For all interactions we start from corresponding terms of the Breit-Pauli (BP) Hamiltonian.26,27 In order to estimate individual contributions, we use perturbation theory (PT) to obtain parameters of standard phenomenolgical Hamiltonians for each interaction treated as perturbation. 1. The Spin-Spin Interaction. The SS interaction appears in systems with S g 1, and it lifts the degeneracy corresponding to states with different MS. In our case, it contributes to the splitting of CT triplet components. The SS BP Hamiltonian26 can be written as a sum of the Fermi-contact term 2
ˆ SS(FC) ) - 4πR H 3
∑ δ(rij)sˆ(i)sˆ(j)
(7)
i*j
and the dipole-dipole term 2 ˆ SS(DD) ) R H 2
∑ i*j
[
3(sˆ(i)rij)(sˆ(j)rij) sˆ(i)sˆ(j) rij3 rij5
]
(8)
Here, R ≈ 1/137 is the fine structure constant, indices i and j enumerate electrons, sˆ(i) is the one-electron spin operator, and rij is the interelectron coordinate. The quadratic form of the SS BP Hamiltonian with respect to the spin coordinates allows us to substitute these terms by the phenomenological or spin Hamiltonian within the triplet subspace
ˆ SS ) SˆDSˆ H
(9)
Role of Protein in Magnetic Field Sensing
J. Phys. Chem. A, Vol. 113, No. 44, 2009 12279
Figure 3. Relation between (D,E) parameters of eq 12 and the firstorder spin-spin splitting of a triplet state.
density matrices or using the resolution of the identity technique.28 Although the CIS two-particle density matrix can be written in terms of one-particle densities,31 here we use an even simpler approach of the mean-field approximation (MFA) where the CIS wave function is treated as if it were a single determinant, and hence only a one-particle density matrix of an excited state is sufficient for the D tensor evaluation. The error introduced by this approximation is considered small.29,30,32 The D tensor within the MFA can be obtained by contracting the one-particle spin density of the CT triplet state (P(R-β)) with the appropriate two-electron integrals
where Sˆ is the many-electron spin operator, and D is the traceless tensor
Dkl )
R Dlk ) 〈T+1 | 2 2
∑
[δlkrij2 - 3(rij)l(rij)k] rij5
i*j
ˆ SS ) DXXSX2 + DYYSY2 + DZZSZ2 H
ˆ SS ) D Sˆ02 - 1 Sˆ2 + E[Sˆ+2 + Sˆ-2] H 3
]
where µ, ν, λ, and σ are Gaussian AOs. Another simplification originating from MFA is a cancellation of the SS Fermi-contact component
Dkl(FC) )
(17)
because of the permutational symmetry of the overlap integrals (〈µν|δ(r12)|λσ〉 ) 〈µσ|δ(r12)|λν〉). Using this permutational symmetry eq 16 can be rewritten as27
(12)
with spherical-tensor operators defined as
√2
∑
Pµσ(R-β)Pλν(R-β)]〈µν|δ(r12)|λσ〉
Dkl )
SˆX + iSˆY
4πR2δkl [Pµν(R-β)Pλσ(R-β) 3 µνλσ
(11)
or to the more common spherical-tensor form27-30
(13)
(
R2 [P (R-β)Pλσ(R-β) - Pµσ(R-β)Pλν(R-β)] × 2 µνλσ µν ∂ ∂ ∂ ∂ + + 〈µν|r12-1 |λσ〉 (18) ∂Rk(µ) ∂Rk(ν) ∂Rl(λ) ∂Rl(σ)
∑
)(
)
where R(µ),R(ν),R(λ) and R(σ) are centers of Gaussian AOs µ, ν, λ, and σ, respectively. Thus, calculation of D elements requires contractions similar to those required in direct Hartree-Fock second derivative calculations. 2. The Spin-Orbit Interaction. The SO part of the BP Hamiltonian also contains one- and two-electron terms26
Sˆ0 ) SˆZ
Sˆ- )
(16)
[2sˆz(i)sˆz(j) -
calculated with the high-spin triplet wave function |T+1〉. The spin Hamiltonian can be transformed to the principle axes (X,Y,Z) where
Sˆ+ ) -
∑
Pµσ(R-β)Pλν(R-β)]〈µν|r12-5{δklr122 - 3(r12)k(r12)l}|λσ〉
sˆx(i)sˆx(j) - sˆy(i)sˆy(j)]|T+1〉 (10)
[
R2 [P (R-β)Pλσ(R-β) 2 µνλσ µν
SˆX - iSˆY
√2 ˆ BP ) H ˆ BP1 + H ˆ BP2 H
Parameters D and E in eq 12 can be expressed as
3DZZ D) 2
(14)
1 E ) (DXX - DYY) 2
(15)
and their relation to the SS splitting of the triplet levels is illustrated by Figure 3. Essentially, the diagonalization of the spin Hamiltonian H or D tensor is equivalent to that of the BP Hamiltonian in the triplet subspace, and therefore, parameters D and E are used to assess the first-order PT contribution of the SS interaction. In general, calculating matrix elements with the SS BP operators (eq 10) would require either forming two-particle
2 ˆ BP1 ) R H 2
ˆ BP2 ) - R H 2
2
(19)
Z
∑ r N3 [rˆiN × pˆ(i)]sˆ(i) i,N
(20)
iN
∑ ∑ r13 [rˆij × pˆ(i)]sˆ(i) + r23 [rˆij × pˆ(j)]sˆ(i) i
i*j
ij
ij
(21) where index N enumerates nuclei, ZN is the nuclear charge, pˆ(i) is the one-electron momentum operator, and riN is the electron-nuclear coordinate. In order to eliminate the computational cost related to the two-electron component of the SO BP Hamiltonian, we use an effective one-electron spin-orbit operator
12280
J. Phys. Chem. A, Vol. 113, No. 44, 2009 2 ˆ SO ) R H 2
ZNeff
∑r
3 iN
i,N
[rˆiN × pˆ(i)]sˆ(i)
Izmaylov et al.
(22)
states
where ZNeff is an effective charge of nucleus N and the difference ZNeff - ZN accounts for the missing two-electron SO contribution in eq 22.33 In terms of PT, generally, the SO interaction contributes to the second order, but in the case of small nonrelativistic singlet-triplet splittings it can contribute to the first order as well. In our case of spatially separated radical pairs, the singlet-triplet CIS gaps are quite tiny (see below); therefore, the SO contribution is considered within the second and first orders of quasi-degenerate PT.29,34 The Hermitian SO perturbation matrix up to the second order is
ˆ SO(1-2) |ΦJ〉 ) 〈ΦI |H ˆ SO |ΦJ〉 + 〈ΦI |H 1 1 1 ˆ SO |ΦP〉〈ΦP |H ˆ SO |ΦJ〉 〈ΦI |H + 2 P EI - EP EJ - EP (23)
[
∑
]
where ΦI and ΦJ are CT states of interest (singlet and three components of triplets), ΦP’s are other CIS solutions, and EI, EJ, and EP are the corresponding CIS energies. The secondorder sum over states in eq 23 has been substituted in actual calculations by a sum over the first 60 excited states. Only the triplet states with MS ) 0 are formed in our CIS calculations; ˆ SO operator in eq thus, to evaluate all matrix elements of the H 34 23 we use the Wigner-Eckart theorem for the singlet-triplet couplings and in addition spin shift operators (Sˆ() for the triplet-triplet couplings. In order to present our study in a more concise manner, 16 matrix elements of the SO perturbation matrix (eq 23) are gathered into three groups: (1) the singlet-singlet element that shifts the singlet energy, (2) the singlet-triplet block that creates the coupling, and (3) the triplet-triplet block which acts similarly as the first-order SS contribution. Because of this similarity, the SO and SS corrections to the splitting of the triplet components are considered together by adding the SO triplet-triplet block to the D tensor prior to its diagonalization. It is worth mentioning that because of the difference in spin-symmetry of involved states, the SO triplet-triplet block cannot be directly summed with the D tensor and needs to be transformed first. Appropriate transformations have been described for the general case of arbitrary spin states by Neese et al.35 However, in our case of S ) 1, we take a simpler approach: the SO triplet-triplet block requires only a unitary transformation from the spherical-tensor basis to the Cartesian basis to be summed directly with the D tensor. 3. The Hyperfine Interaction. The HF term originates from interactions of electron spins with spins of magnetic nuclei. The corresponding Hamiltonian can be split into the isotropic (or spherically symmetric)
ˆ HF(iso) ) geβe 8π H 3
TABLE 1: Characteristics of CT States for the FAD•Mechanism Calculated with the CIS Method
∑ gNβNsˆ(i)Iˆ(N)δ(riN)
(24)
Mulliken charges, a.u. Trp377
Trp324
T S T S T S
Two-Fragment Model/6-31G* 6.1307 -1.002 1.002 6.1313 -1.000 1.000 6.7697 -1.006 1.006 6.7697 -0.994 0.994 7.0880 -0.968 7.0880 -1.032
0.968 1.032
T S T S T S
6.1017 6.0816 6.7419 6.7419 6.9173 6.9173
Two-Fragment Model/EPR-II -0.921 0.921 -0.843 0.843 -1.000 1.000 -1.000 1.000 -1.001 -0.999
1.001 0.999
T S T S T S
Four-Fragment Model/6-31G* 6.0889 -1.004 0.994 0.009 6.0905 -0.870 0.861 0.008 6.6997 -1.005 0.058 0.946 6.6997 -1.005 0.057 0.947 7.1688 -1.004 0.003 0.001 7.1688 -1.006 0.003 0.001
0.001 0.001 0.001 0.001 1.000 1.002
spin
energy, eV
FAD
Trp400
parts. Here, Iˆ(N) is the nuclear spin operator, ge, gN, βe, and βN are g-factors and Bohr magnetons for electrons and nuclei. Using the same approach as in the SS case, the HF spin Hamiltonian can be written as
ˆ HF ) SˆAIˆ H
(26)
where Iˆ is many nuclear spin operator and A is the HF coupling tensor. A components for each nuclei can be further split into two parts: the isotropic
AN,kl(iso) ) δkl
8π gegNβeβN 〈Φ| 3 〈Φ|Sz |Φ〉
∑ δ(riN)sˆz(i)|Φ〉 i
(27) and the dipolar
AN,kl(dip) )
gegNβeβN 〈Φ| 〈Φ|Sz |Φ〉
∑ i
[
δkl riN3
-
3(riN)k(riN)l riN5
]
sˆz(i)|Φ〉
(28)
where kl are Cartesian components and |Φ〉 is an electron wave function with nonzero Sz. As in the SS case, contractions (27) and (28) are done with the CT triplet state P(R-β) densities and corresponding one-electron integrals. Usually, the HF interaction is considered within the first order of PT,27 and thus in the current study components of the HF coupling tensor A will be used to assess the HF contribution.
i,N
III. Results and Discussions and the dipolar (or anisotropic)
ˆ HF(dip) ) -geβe H
[
ˆ
3(sˆ(i)riN)(Iˆ(N)riN)
iN
riN5
(N) ∑ gNβN sˆ(i)I r 3 i,N
]
(25)
Tables 1 and 2 present energies and Mulliken charges on each fragment of the lowest CT states for both mechanisms. The CT states have such a pronounced difference in charge distribution compared to other states that we would not expect a change in the assignment of CT states due to using any other standard scheme for extracting partial charges. Two- and four-fragment
Role of Protein in Magnetic Field Sensing
J. Phys. Chem. A, Vol. 113, No. 44, 2009 12281
TABLE 2: Characteristics of CT states for the FADH+ mechanism calculated with the CIS method states
Mulliken charges, a.u. Trp377
Trp324
T S T S T S
Two-Fragment Model/6-31G* 2.8125 -0.003 1.003 2.8123 -0.008 1.008 2.6033 0.003 0.997 2.6033 -0.003 1.003 2.5570 0.001 2.5570 -0.001
0.999 1.001
T S T S T S
2.9027 2.9029 2.6806 2.6806 2.6349 2.6349
Two-Fragment Model/EPR-II -0.001 1.001 -0.011 1.011 0.000 1.000 0.000 1.000 0.000 0.000
1.000 1.000
T S T S T S
Four-Fragment Model/6-31G* 2.6131 -0.008 0.001 0.915 2.6132 -0.009 0.001 0.916 2.7367 -0.009 0.003 0.005 2.7367 -0.009 0.003 0.005 2.7619 -0.008 0.914 0.093 2.7624 0.004 0.903 0.092
0.092 0.091 1.001 1.001 0.001 0.001
spin
energy, eV
FADH
Trp400
models produce qualitatively similar results. On the other hand, the two studied mechanisms have significantly different excitation energies. Energies required for the electron hopping in the FADH+ mechanism (Table 2) correspond exactly to the blue part of the visible light spectrum (470 nm or 2.67 eV), while those in the FAD•- mechanism are beyond the visible range (
