Originally published as Biophys J. BioFAST on July 21, 2006.
doi:10.1529/biophysj.105.079483
Biophysical Journal 91:2778-2797 (2006)
© 2006 The Biophysical Society
How Proteins Trigger Excitation Energy Transfer in the FMO Complex of Green Sulfur Bacteria
Julia Adolphs and
Thomas Renger
Institut für Chemie und Biochemie, Freie Universität Berlin, Berlin, Germany
Correspondence: Address reprint requests to Dr. Thomas Renger, Tel.: 00-49-030-838-54907; E-mail: rth{at}chemie.fu-berlin.de.
 |
ABSTRACT
|
|---|
A simple electrostatic method for the calculation of optical transition energies of pigments in protein environments is presented and applied to the Fenna-Matthews-Olson (FMO) complex of Prosthecochloris aestuarii and Chlorobium tepidum. The method, for the first time, allows us to reach agreement between experimental optical spectra and calculations based on transition energies of pigments that are calculated in large part independently, rather than fitted to the spectra. In this way it becomes possible to understand the molecular mechanism allowing the protein to trigger excitation energy transfer reactions. The relative shift in excitation energies of the seven bacteriochlorophyll-a pigments of the FMO complex of P. aestuarii and C. tepidum are obtained from calculations of electrochromic shifts due to charged amino acids, assuming a standard protonation pattern of the protein, and by taking into account the three different ligand types of the pigments. The calculations provide an explanation of some of the earlier results for the transition energies obtained from fits of optical spectra. In addition, those earlier fits are verified here by using a more advanced theory of optical spectra, a genetic algorithm, and excitonic couplings obtained from electrostatic calculations that take into account the influence of the dielectric protein environment. The two independent calculations of site energies strongly favor one of the two possible orientations of the FMO trimer relative to the photosynthetic membrane, which were identified by electron microscopic studies and linear dichroism experiments. Efficient transfer of excitation energy to the reaction center requires bacteriochlorophylls 3 and 4 to be the linker pigments. The temporal and spatial transfer of excitation energy through the FMO complex is calculated to proceed along two branches, with transfer times that differ by an order of magnitude.
 |
INTRODUCTION
|
|---|
In photosynthetic antennae, light is absorbed by pigments and the excitation energy is transferred to the photosynthetic reaction center. The pigments are held by a protein scaffold at the right distances and orientations for efficient excitation energy transfer. A key question of the structure-function relationships in photosynthetic antennae is: How does the protein environment influence the optical transition energies of the pigments? Besides static and dynamic disorder caused by slow and fast protein dynamics, respectively, there is a specific change of the transition energy of every pigment due to the different protein environment. Different pigment-protein interactions are assumed to influence the optical transition energies of the pigments:- Electrostatic interactions with charged and aromatic amino-acid side chains (1
,2
).
- Differences in Mg2+-ligation (3
).
- Differences in H-bonding (4
7
).
- Differences in the conformation of the pigment; for example, the degree of planarity and the rotation angle of the acetyl group (8
,9
).
A direct calculation of transition energies of pigments in the protein, the so-called site energies, is difficult (2
,9
) because of the complexity of the above interactions. Therefore, the site energies are often treated as parameters that are determined from the fit of optical spectra.
There is one pigment-protein complex, the so-called Fenna-Matthews-Olson (FMO) complex of green sulfur bacteria, which has been extensively used as a model system for larger antenna complexes, starting more than 25 years ago with the pioneering work of Pearlstein (10
) and co-workers. The FMO complex of Prosthecochloris aestuarii was the first pigment-protein complex for which the structure was determined by x-ray crystallography (11
,12
). In the meantime, the resolution of the electron density map has been refined to 1.9 Å (13
) and the structure of the FMO complex of a strongly related bacterium Chlorobium tepidum has been determined as well (14
,15
). The two structures are very similar, but interestingly, the spectra look different. The molecular origin of this difference is unknown; an explanation is suggested in this study. The structure of the FMO complex consists of a trimer, formed by three identical monomers that each bind seven bacteriochlorophyll-a (BChla) molecules, as shown in Fig. 1. The usual numbering of the BChls (original chosen by Fenna and Matthews (11
)) is used throughout this article.

View larger version (21K):
[in this window]
[in a new window]
|
FIGURE 14 Disorder-averaged population dynamics of exciton states calculated in Redfield theory, where the initial population was created assuming that the excitation energy arrives from above the trimer in Fig. 1 as explained in detail in the text. The inset contains on enlarged view on the occupation probabilities of exciton states 46.
|
|
The FMO trimer mediates excitation energy transfer between the chlorosomes, which are the main light-harvesting antennae of green sulfur bacteria and the membrane-embedded type I reaction center (RC) (16
). There is some structural information on the FMO-RC super complex, obtained in an electron microscopic study (17
,18
) and from linear dichroism (LD) spectra measured on FMO trimers and FMO-RC complexes. The latter provide evidence that the symmetry axis of the trimer is parallel to the normal of the membrane, containing the reaction center and the former gave the overall mutual arrangement of the two complexes as sketched in the upper-right corner in Fig. 1. The resolution of the electron microscopic study does not allow us to distinguish between the top and the bottom of the FMO trimer. Hence, either BChls 1 and 6 or BChls 3 and 4 are the linker pigments between the FMO complex and the RC. Since the RC represents an energetic sink, it can be assumed that the low energy pigments serve as linkers.
Different sets of site energies of the seven BChla molecules have been extracted in the past from calculations of optical spectra (19
25
), using different fit algorithms or simply by hand. The different calculations predict different pigments to be responsible for the lowest exciton state. The early exciton calculations by Pearlstein (19
) suggested BChl 7; the calculations by Gülen (21
) found the lowest site energy for BChl 6, and the calculations of Louwe et al. (22
), Vulto et al. (24
), Renger and May (23
), and Wendling et al. (25
) favor BChl 3. Based on the finding that the surface of the FMO trimer is more hydrophobic on the site of BChl 6 it was suggested that the FMO trimer is partially embedded in the membrane (14
), a suggestion that would favor BChl 6 to be the linker between the FMO trimer and the reaction center. However, from the electron microscopy study on FMO-RC complexes, it was seen that the FMO complex is not part of the membrane (17
,18
), leaving again both options for the lowest energy pigment open. In the present study, two independent methods for the calculation of site energies provide compelling evidence that BChls 3 and 4 dominate the lowest exciton state and, for efficient excitation energy transfer to the reaction center, are the linker pigments. The temporal and spatial relaxation of excitons through the FMO complex toward the RC that results from the site energies determined in this study is shown in the lower part of Fig. 1 and will be explained and discussed in detail later.
Until now, the best agreement between calculated and experimental optical spectra of the FMO complex has been obtained by Vulto et al. (22
,24
) for C. tepidum and by Wendling et al. (25
) for P. aestuarii. In both cases BChl 3 was identified as the one with the lowest site energy. The crucial point of those studies (22
,24
,25
) was the use of a smaller dipole strength of BChla (22
) in the calculations of the excitonic couplings than assumed previously. A physical explanation of the rather small couplings is still missing and will be given here. Concerning the calculations of Vulto et al. (24
) it was not clear whether the values for the site energies would change if a better theory for optical spectra were used, since lifetime broadening, vibrational sidebands, and resonance energy transfer narrowing of the optical lines were neglected (24
). It is shown here that this rather harsh approximation does not lead to drastic changes in the relative shifts in optimized site energies. However, the quality of the fit of the spectrum is better if a more advanced theory of optical spectra (26
) is used.
Insights into the dynamics of excitons in the FMO-complex are obtained from time-resolved nonlinear optical spectroscopy (27
32
), which detected subpicosecond as well as picosecond exciton relaxation times. In contrast to pump-probe studies, the newly developed two-dimensional photon-echo technique (31
,32
) allows one to directly detect the coupling between different exciton levels and the energies of exciton states, and to infer from those quantities the spatial and temporal relaxation of excitons. However, such a structural interpretation still requires knowledge of the site energies and excitonic couplings of the pigments. In the recent analysis of the photon-echo data of the FMO complex (31
,32
), the couplings of Vulto et al. (24
) were used, with one exceptionthe excitonic coupling between BChl 5 and 6 was reduced from 94 cm1 to 40 cm1. The site energies by Vulto et al. (24
) were readjusted to fit the linear absorption and two-dimensional photon-echo spectra (31
,32
). In this study, we are aiming at a direct calculation of the excitonic couplings and the site energies.
There is just one study on the FMO complex that tried to obtain the relative shift of site energies directly, from semiempirical quantum chemical calculations (9
). Different factors as nonplanarity of the BChl, rotation of the acetyl group, and charged amino acids (all charged amino acids within 5.5 Å of the pigments were taken into account) were identified, but no optical spectra were calculated. The overall spread of the site energy values (9
) is too large by a factor of 25, judging by the width of the experimental spectra (25
). We find that, even when the site energy differences are downscaled, the agreement between the calculated and experimental spectra is rather poor.
The article is organized in the following way. The theory used for the calculation of optical spectra and exciton relaxation is summarized first. Next, the effect of the dielectric medium of the protein on the excitonic couplings of the pigments is investigated by electrostatic calculations. Afterwards, the site energies of the pigments are obtained by a fit of optical spectra, using a genetic algorithm and, independently, by the calculation of electrochromic shifts due to charged amino acids. The site energies and excitonic couplings are then used to calculate the temporal and spatial relaxation of excitons after a short pulse excitation or assuming that excitation energy arrives from the direction where the chlorosomes are situated in green sulfur bacteria.
 |
THEORY
|
|---|
Hamiltonian
The theory is based on a standard Hamiltonian Hppc for the pigment protein complex, that describes the pigments as coupled two-level systems interacting with vibrational degrees of freedom of the pigments and the protein,
 | (1) |
The exciton part,
 | (2) |
contains the site energies Em of the pigments, defined as the optical transition energies at the equilibrium position of nuclei in the electronic ground state, and the excitation energy transfer couplings Vmn.
The Hamiltonian Hexvib describes the modulation of site energies by the vibrations. A linear dependence of the site energies on the (dimensionless) vibrational coordinate Q
is assumed, which is defined in terms of creation and annihilation operators of vibrational quanta, Q
=
(26
,33
):
 | (3) |
The dimensionless coupling constants g
(m) = g
enter the spectral density J(
) of the exciton vibrational coupling
 | (4) |
which is the key quantity in the expressions for optical spectra and the rate constants for exciton relaxation discussed below. The value J(
) is assumed independent on the site index m; i.e., the same local modulation of site energies by the vibrational dynamics is assumed.
The vibrations are described by an Hamiltonian Hvib of harmonic oscillators
 | (5) |
where Tnucl is the kinetic energy of nuclei.
The coupling between the pigment-protein complex and an external radiation field is described by the semiclassical Hamiltonian Hppc-rad, which reads in rotating wave approximation
 | (6) |
where
is the molecular transition dipole moment of the mth pigment,
the polarization of the field, and h.c. the Hermitian conjugate. The field may be either stationary, i.e., E
(t) = E0, or time-dependent. In the latter case, a Gaussian shape for E
(t) is assumed as
 | (7) |
where the full width at half-maximum (FWHM) of E
(t) is
.
For the calculations of optical spectra and exciton relaxation, the above Hamiltonian Hppc + Hppcrad is expressed in terms of delocalized exciton states |M
, which are given as linear combinations of localized excited states
, where
describes the probability that the mth pigment is excited when the PPC is in the Mth exciton state. The exciton coefficients
and excitation energies
are obtained from the solution of the eigenvalue problem
 | (8) |
with the Hex in Eq. 2. The Hamiltonian Hexvib in Eq. 3, in the basis of delocalized exciton states, becomes
 | (9) |
with the coupling constant
 | (10) |
that contains now diagonal (M = N) as well as off-diagonal (M
N) parts. The former give rise to vibrational sidebands of exciton transitions in optical spectra and the latter lead to relaxation between different exciton states.
The coupling to the radiation field Hppcrad in Eq. 6 now reads
 | (11) |
where the transition dipole moments
of the delocalized exciton states are obtained from the local transition dipole moments
and the exciton coefficients
as
 | (12) |
The above Hamiltonians lead to the expressions for linear optical spectra, exciton state occupation probabilities, created by a short pulse, and rate constants of exciton relaxation, described in the following.
Linear optical spectra
The linear absorption
(
) is obtained from the dipole-dipole correlation function as explained in detail in Renger and Marcus (26
),
 | (13) |
where DM(
) is the lineshape function and
dis denotes an average over static disorder in site energies. A Gaussian distribution function of width (FWHM)
dis is assumed for these energies, and the disorder average is performed by a Monte Carlo method. In the calculation of circular dichroism, the dipole strength
in Eq. 13 is replaced by the rotational strength rM and in the case of linear dichroism by
, where
M is the angle between the symmetry axis of the trimer and the excitonic transition dipole moment
.
The lineshape function DM(
) was obtained using a non-Markovian partial-ordering prescription theory. It is given as (26
)
 | (14) |
where
denotes the real part of the integral. The value DM(
) contains both vibrational sidebands and lifetime broadening due to exciton relaxation. The vibrational sidebands are described by GM(t) and the lifetime broadening by the dephasing time
M (discussed in detail below). Both quantities are related to the spectral density J(
) in Eq. 4. The time-dependent function GM(t) in Eq. 14 is given as
 | (15) |
with
 | (16) |
and
MM being the diagonal part of
 | (17) |
It contains the exciton coefficients
, a correlation radius Rc of protein vibrations (26
) (note that a value of Rc = 5 Å is used, determined from transient spectra of Photosystem II reaction centers in (34
); the stationary spectra calculated here do not depend critically on this value) and the center-to-center distance Rmn between pigments m and n.
The function n(
) in Eq. 16 is the mean number of vibrational quanta with energy 
that are excited at a given temperature T,
 | (18) |
where k is Boltzmann's constant. The dephasing time
M in Eq. 14 is determined by the rate constants kM
N of exciton relaxation, obtained using a Markov approximation for the off-diagonal parts of the exciton vibrational coupling,
 | (19) |
where the rate constant reads
 | (20) |
It resembles the standard Redfield rate constant (e.g., (35
)). The
is the real part of the Fourier-Laplace transform of the pigment's optical energy gap correlation function (26
),
 | (21) |
where J(
) = 0 for
< 0, and
MN =
M
N is the transition frequency between the Mth and the Nth exciton states. The
in Eq. 14 is shifted from the purely electronic transition frequency
M due to the exciton-vibrational coupling
 | (22) |
where
MN is given in Eq. 17 and E
is the local reorganization energy
 | (23) |
The
in Eq. 22 is related to the real part
in Eq. 21 by a Kramers Kronig relation (26
,36
)
 | (24) |
where
denotes the principal part of the integral (details are given in the Supplementary Material).
In the present model, the spectral density J(
) contains both a broad low frequency contribution S0J0(
) by the protein vibrations with Huang-Rhys factor S0 and a single effective high-energy vibrational mode of the pigments with Huang-Rhys factor SH:
 | (25) |
For the normalized low-frequency function J0(
) (with
), we assume that it has the same form as the spectral density, extracted recently (26
) from 1.6 K fluorescence line narrowing spectra of B777-complexes,
 | (26) |
with the extracted parameters s1 = 0.8, s2 = 0.5, 
1 = 0.069 meV, and 
2 = 0.24 meV. The Huang-Rhys factor S0 of the pigment-protein coupling in Eq. 25 was estimated from the temperature dependence of the absorption spectra (28
,37
) of the FMO complexes of P. aestuarii and C. tepidum to be
0.5. A comparison of the function J0(
) with the fluorescence line narrowing spectrum of C. tepidum (37
) is shown at the top of Fig. 2. The experimental spectrum (for the present weak exciton-vibrational coupling) should be similar (26
) to the J(
) in Eq. 25. On the basis of this comparison, the Huang-Rhys factor and the energy of the high-frequency mode were estimated as SH = 0.22 and
H = 180 cm1. By taking into account that the energy of the high-frequency mode is large compared to the thermal energy (at T = 5 K), i.e., 
H >> kT, the function GM(t) in Eq. 15 becomes, using Eqs. 16 and 25,
 | (27) |
By introducing this GM(t) into Eq. 14, using a series expansion for
, the following lineshape function DM(
) is obtained:
 | (28) |
The above DM(
) differs from the one obtained previously (26
) in that it includes now also vibrational satellites of a high-frequency vibrational mode. The dephasing time
M of the Mth exciton transition and the off-diagonal part of the shift in transition energy
in Eq. 22 are determined by the low frequency contribution S0J0(
) of the spectral density.
Relation of the present theory to the earlier theories by Vulto et al. (24
) and Wendling et al. (25
)
The theory used by Vulto et al. (24
) is obtained from Eq. 13 by neglecting any homogeneous broadening, i.e., lifetime broadening and vibrational sidebands, by setting GM(t) = 0 and
in Eq. 14. In this case, the line-shape function DM(
) becomes just a delta-function that peaks at the exciton transition frequency
M, DM(
) =
(
M). To perform the average over disorder analytically, any resonance energy transfer narrowing (38
), i.e., a decrease of the width of the distribution of exciton energies with respect to that of the distribution of local transition energies, is neglected as well. Assuming a Gaussian distribution function
of the same width (FWHM)
for all exciton energies 
M results in
 | (29) |
The mean exciton energies
are those obtained for the mean site energies.
In the theory of Wendling et al. (25
,37
), resonance energy transfer narrowing, lifetime broadening, and vibrational sidebands were taken into account using the following approximations for the function DM(
) in Eq. 14:
- A shift of the transition frequency
M by the exciton vibrational coupling was neglected, i.e.,
in Eq. 22.
- Any dependence of the vibrational sidebands of the exciton transitions on the delocalization of exciton states was neglected by replacing the function GM(t) =
MMG(t) in Eq. 14 by the function G(t), which describes the vibrational sidebands of a monomeric BChla in a protein environment.
- The factor
appearing in Eqs. 19 and 21 for the rate constant of exciton relaxation between the Mth and the Nth state was approximated by a constant
0, which is assumed the same for all transitions M
N, i.e., instead of the
in Fig. 2 (bottom), a scaled J(
) at the top of this figure is used.
The main differences between the present and the two earlier theories are illustrated in Fig. 3. Due to the neglect of vibrational sidebands in the theory of Vulto et al. (24
) (dotted line), there is less intensity in the blue part of the spectrum. The missing lifetime broadening seems to be compensated in part by the neglect of resonance energy transfer narrowing as seen by the width of the two main peaks. In the theory of Wendling et al. (25
,37
), the neglect of the dependence of the vibrational sideband on the exciton states delocalization leads to a stronger vibrational sideband as seen in the blue part of the spectrum and to a change of the relative height of the two main peaks. We note that the differences between the three theories are partly compensated in the fit of optical spectra by different site energies. However, as will be seen below, the overall ranking of optimal site energies is similar in all three approaches.
Excitation by a short laser pulse
The Hamiltonian Hppc in Eqs. 111 is rewritten, using the completeness relation,
, as
where the diagonal part of the exciton-vibrational coupling was used to construct potential energy surfaces (PES) of exciton states
 | (31) |
that are shifted along the coordinate axis by 2g
(M,M) with respect to the PES U0(Q) of the ground state
 | (32) |
The energy 
'M is the transition energy between the minima of the excited state and the ground state PES,
 | (33) |
where 
M is the vertical transition energy (site energy) and E
the local reorganization energy in Eq. 23. It is assumed that the off-diagonal parts of the exciton-vibrational coupling are weak enough that the exciton relaxation during the short excitation pulse can be neglected. In this case, a nonperturbative inclusion of the diagonal part of the exciton vibrational coupling and a second-order perturbation theory with respect to Hppcrad in Eq. 11 yields for the population PM(t) of the Mth exciton state
 | (34) |
where Trvib denotes a trace with respect to the vibrational degrees of freedom, and the vibrationally relaxed initial electronic ground state is described by the equilibrium statistical operator
 | (35) |
and an average over a random orientation of complexes with respect to the polarization of the external field was performed. By changing the integration variable
1
1, and setting t0
, the occupation probability is obtained as
 | (36) |
where
denotes the real part and f(t,
1) is given as
 | (37) |
which for t
becomes the autocorrelation function of the pulse. The function g(
1) contains the average over the vibrational degrees of freedom that is performed using a second-order cumulant expansion, which is exact for harmonic oscillators (33
)
 | (38) |
with the
in Eqs. 15 and 16.
The population of exciton states after the action of the short pulse
is obtained by formally setting t
in Eq. 37. With the E
(t) in Eq. 7, the population
reads
 | (39) |
This result will be used as an initial condition in the calculation of exciton relaxation, discussed in the following.
Exciton relaxation dynamics
Exciton relaxation is described by modified Redfield theory (39
42
). It is assumed that before exciton relaxation between different delocalized states the nuclei relax in the respective excitonic PES. The rate constant kM
N is obtained in second-order perturbation theory in the off-diagonal part of the exciton-vibrational coupling. In contrast to Redfield theory, which neglects nuclear reorganization effects that accompany exciton relaxation, in modified Redfield theory those effects are taken into account by a nonperturbative inclusion of the diagonal part of the exciton-vibrational coupling. Within the present harmonic oscillator model the rate constant is obtained as (41
,42
)
 | (40) |
where the time-dependent functions
 | (41) |
 | (42) |
 | (43) |
are related to the spectral density S0J0(
) via the function
k(t), with k = 0, 1, 2,
 | (44) |
The time-independent part in the integrand in Eq. 40,
MN, is
 | (45) |
using the local reorganization energy E
in Eq. 23 with J(
) = S0J0(
). The coefficients aMN, bMN, cMN, and dMN in the above equations are given by the exciton coefficients and the correlation radius of protein vibrations as (42
)
 | (46) |
 | (47) |
 | (48) |
 | (49) |
If the diagonal part of the exciton vibrational coupling, i.e., the mutual shift of excitonic PES surfaces, is neglected, Eq. 40 reduces to the Redfield result in Eq. 20: The only function in Eq. 40 that does not contain a diagonal part of the coupling is the FMN(t) in Eq. 43 (41
). After setting the remaining functions to zero, Eq. 40 becomes
. If the FMN(t) in Eq. 43 is introduced and the integration over
is carried out, the rate constant becomes
. By noting that (1 + n(
MK)) = n(
KM), the equality of the above rate constant with the Redfield result in Eq. 20 is seen.
Exciton relaxation dynamics is described by the rate equations for the populations PM(t) of exciton states,
 | (50) |
with the initial populations PM(0) =
created by the short pulse (Eq. 39).
Alternatively, exciton states which are formed by the pigments at the top of the trimer in Fig. 1 will be populated at time zero to mimic excitation energy transfer as it occurs in vivo between the chlorosomes and the RC. For the rate constants kM
N the modified Redfield result in Eq. 40 or the Redfield expression in Eq. 20 is used. In matrix form, Eq. 50 reads
, where the Mth element of
is PM(t) and the kinetic matrix A contains in the diagonal the elements
and in the off-diagonal
. The standard solution for
is given as
, where the
i and
are the eigenvalues and (right) eigenvectors of the kinetic matrix A and the constants di are obtained from the initial condition
.
 |
EXCITONIC COUPLINGS IN THE DIELECTRIC ENVIRONMENT OF THE PROTEIN
|
|---|
Although the excitonic couplings, in contrast to the site energies, can be calculated straightforwardly from the structural data, there is an unknown scaling factor that contains the uncertainty about the effect of the dielectric environment on the Coulomb (excitonic) interaction. In point dipole approximation, the excitonic coupling is given as
 | (51) |
where
is a unit vector along the transition dipole moment of the mth BChl, the unit vector
is oriented along the line connecting the centers of BChls m and n, and
is the dipole strength of the Qy transition of BChla in vacuum. The factor f describes the enhancement of the dipole strength and the screening of the Coulomb coupling by the dielectric environment with optical dielectric constant
,
 | (52) |
If the distance between the pigments is large compared to the extension of their ground and excited state wavefunctions, the screening factor g equals 1/
. If two pigments are very close, it can be expected that the factor f becomes distance-dependent, because the two pigments do not form separate cavities with dielectric medium in between, but are instead situated in the same cavity. As pointed out by Knox and van Amerongen (43
), the enhancement of the dipole strength, i.e., the change in quantum mechanical transition probability by the dielectric, can be classically understood by the change of local electric field interacting with the vacuum transition dipole. Recently, Knox and Spring (44
) analyzed the dipole strength of BChla in 15 different solvents and found the empirical relation
between the dipole strength µ(
)2 in the solvent with dielectric constant
and the vacuum value µ(1)2. If this empirical formula is compared with the predictions of two cavity models, the empty cavity model was found to provide the best explanation (44
). Within the latter, the dipole strength of BChla in the solvent is obtained as µ(
)2 = 37.1 (3
/(2
+ 1))2 D2, that contains the vacuum dipole strength
.
In a recent study, Hsu et al. (45
) demonstrated that the excitonic coupling between two molecules in a dielectric environment, resulting from a quantum mechanical treatment with time-dependent density functional theory, can also be obtained classically. The classical calculation just involves the electrostatic coupling between the transition densities of the two molecules taking into account the fast (optical) part of the dielectric response of the medium. Hsu et al. (45
) showed, using a multipole expansion, that the interaction between two transition dipoles in a spherical cavity can be either increased or decreased by the dielectric response of the environment, depending on the geometry of the transition dipoles.
As described in detail below, the present method is based on a numerical solution of the Poisson equation and allows us thereby to also treat nonspherical cavities, as those of the BChls in the FMO-protein. We have tested our numerical procedure by reproducing the spherical cavity results of Hsu et al. (45
), as shown in Fig. 4. In this calculation, the transition dipoles are represented by two transition charges of opposite sign at close distance. The two geometries are depicted in Fig. 4 (left) and the resulting couplings are shown in Fig. 4 (right) in dependence on the dielectric constant.
In the following, we take into account nonspherical cavities of the BChls in the FMO-protein. The dielectric of the FMO-protein is generated by using overlapping spheres for its atoms, with atomic radii taken from the force field used by the program CHARMM (46
), as illustrated in Fig. 5 (right). An optical dielectric constant of
= 2 is assumed for the protein. The BChls are described by empty cavities that contain transition charges
describing the optical transition. The atomic transition monopole (TM) charges qI of Chang (47
) are placed at the position
of the atom I of the mth BChl. The charges were scaled so as to result in a vacuum transition dipole moment square of 37.1 D2 of the empty cavity analysis by Knox and Spring (44
). The excitonic couplings are obtained from the Coulomb interaction of the transition charges of different BChls. The respective vacuum results are shown in column (1
) of Tables 1 and 2. To calculate the excitonic coupling between the BChls in the dielectric environment, the Poisson equation
 | (53) |
is solved for each BChl numerically by a finite difference method using the program MEAD (48
). The value of
equals 2 if
points to a position in the protein and 1 in the case of BChl. From the resulting electrostatic potential
of BChl m, the excitonic coupling with BChl n is obtained as
.

View larger version (63K):
[in this window]
[in a new window]
|
FIGURE 5 Structure of monomeric subunit of the FMO complex of C. tepidum. (Left) BChls are green; positively charged amino acids (Arg, Lys) are red; negatively charged amino acids (Asp, Glu) are blue; and neutral amino acids are gray. (Right) Dielectric volume used in the calculation of excitonic couplings. The protein (gray) forms a cavity around the BChls (green).
|
|
The results for the couplings (column (3
) in Table 1 and 2) are compared to the values obtained for
= 1 in the transition monopole approximation (column (1
)), and point dipole approximation from Eq. 51, using f = 1.0 (column (2
)) and f = 0.8 (column (4
)). A least mean-square fit,
, for the transition monopole couplings in columns (1
) and (3
) results in f = 0.80. The closest couplings between
(
= 2) and the point dipole approximation (column (2
)) were obtained by using a factor of f = 0.81 (P. aestuarii) and 0.82 (C. tepidum), i.e., we obtain nearly the same correction factor as for the difference between
(
= 2) and
(
= 1). It is clear, thereby, that the point dipole approximation works well and that the main effect on the couplings is due to the dielectric environment. The present factor of 0.800.82 is somewhat larger than the factor 0.72 (
= 2) obtained for point dipoles in a spherical cavity.
A closer examination, using extended dipoles instead of transition monopoles, shows that approximately half of the deviations occur because of the more realistic charge distribution and the remaining half due to the more realistic cavity shape and the fact that cavities of neighboring pigments overlap. However, overall, the simple spherical empty cavity model is seen to work surprisingly well, being off by just
10%. In the calculation of optical spectra and exciton relaxation presented below we use the point dipole approximation and a factor f = 0.8 for both P. aestuarii and C. tepidum.
Finally, we note that the Knox/Spring analysis of the dipole strengths was limited to the 00 transition of the BChla Qy transition, i.e., it does not include excitations of intramolecular vibrations. To take into account excitonic couplings that involve intramolecular vibrational transitions would require us to include those vibrational modes via the respective Franck Condon factors explicitly in the theory (see e.g. (49
,50
)), which is beyond the scope of this article. We have performed dimer calculations with and without an explicit vibrational transition to investigate whether, in a simple theory that does not take into account such a transition, the dipole strength of the latter should be included in the calculation of the excitonic coupling. From a comparison of the spectral position of the main peaks in the absorption spectrum obtained in the sophisticated and the simple theory, it is seen that, just taking into account the excitonic coupling that results from the 0-0 transition dipole strength, gives better results in the simple theory.
 |
CALCULATION OF SITE ENERGIES
|
|---|
In the following, two independent methods are used to obtain the site energies of the seven BChla molecules of the monomeric subunits of the FMO complex of P. aestuarii and C. tepidum. First, the site energies are used as parameters that are optimized by a genetic algorithm in the fit of optical spectra. Afterwards, the site energies are obtained directly by electrochromic shift calculations.
Genetic algorithm
An often successful way to treat nonlinear multidimensional optimization problems is to use a genetic algorithm as a fit routine (51
53
). A rough description of the optimization process is given in Fig. 6. In the initial step, N sets of site energies are generated (a so-called population of N chromosomes), one set is provided (start set), and N 1 are randomly created. The corresponding spectra are calculated, and the different sets are ranked according to their fitness value, which is obtained from the reciprocal of the mean-square deviation between the calculated and experimental spectrum. To get better individuals from generation to generation, genetic operations (51
,52
) known as "recombination" and "mutation" (Fig. 7) are performed, taking the ranking into account. In each step, the chromosome with the highest rank is copied to the next cycle (reproduction). The fit is iteratively completed.
The algorithm was tested first by replacing the experimental spectrum by a calculated spectrum, shown as circles in Fig. 8, which was obtained for the site energies given in the last column (True) of Table 3. A simultaneous fit of the absorption (OD), linear dichroism (LD), circular dichroism (CD), and the derivative of the absorption (OD') spectrum was performed. A typical population of 100 chromosomes was used. After 200 iteration steps, the values for the site energies found by the algorithm were within ±2 cm1 of the true values for convenient initial values and within ±8 cm1 for inconvenient initial values. The algorithm converges fast, if all initial values for the site energies are chosen equal at an energy of 12,460 cm1 in the center of the target spectrum. The spectrum resulting for those initial site energies is shown as a dotted line in Fig. 8. The most difficult situation occurred with an inverted order of the initial values with respect to the optimal values. The respective initial spectrum is shown as a dashed line in Fig. 8. If the OD spectrum alone was fitted, the optimal site energies obtained were wrong, even if a population of 500 chromosomes was used.

View larger version (12K):
[in this window]
[in a new window]
|
FIGURE 8 Test of the genetic algorithm. The two Initial spectra are obtained for different initial sets of site energies used in the genetic algorithm as explained in the text. The True spectrum is used in place of the experimental spectrum; it is obtained for a set of known test values of site energies. The Fit-Results curve (solid line) shows the spectra obtained for the two sets of optimal site energies given in Table 3.
|
|
Fit of low-temperature spectra of C. tepidum and P. aestuarii
The genetic algorithm was used for a simultaneous fit of the OD, OD', CD, and LD spectra of C. tepidum and P. aestuarii, measured by Vulto et al. (24
) and Wendling et al. (25
). The temperature of measurement and simulation was 6 K for C. tepidum and 4 K for P. aestuarii. A width of 100 cm1 for the distribution function of site energies gave the best agreement between the calculated and experimental spectra. The excitonic couplings in column (4
) of Table 1 for C. tepidum and of Table 2 for P. aestuarii were used in the optimization. The fit of site energies was performed for the monomeric as well as for the trimeric structure of the FMO complex, to evaluate the effect of the intermonomer couplings. The optimal site energies obtained from the fits are shown in Table 4. As one can see, the maximum deviation of optimal site energies between monomer and trimer calculations is 70 cm1 (BChl 7 of P. aestuarii), and the ranking of the site energies is almost identical. In accordance with some of the earlier results (24
,25
), pigment number 3 has the lowest site energy for C. tepidum and P. aestuarii and the largest deviation in site energies between the two species is obtained for BChl 5. The overall ranking in site energies obtained here, is almost identical to those obtained earlier by Vulto et al. (24
) and Wendling et al. (25
) as seen in Table 4. The spectra obtained for the present optimal site energies for monomers and trimers are compared in Fig. 9 with the experimental data (24
,25
). Although the spectra for the monomers and the trimers are very similar, there are some quantitative differences for the CD spectrum of P. aestuarii that cannot be ignored. Hence, in the following, all calculations of optical spectra are done for the trimeric structure. The quality of the trimer fit is very good for all spectra of C. tepidum and for OD and LD of P. aestuarii, and somewhat weaker for CD of the latter.

View larger version (20K):
[in this window]
[in a new window]
|
FIGURE 9 Low temperature optical spectra calculated for optimal site energies shown in Table 4. The top panels contain absorption (OD), the center panels contain circular dichroism (CD), and the bottom panels contain linear dichroism (LD) spectra. The experimental data of Vulto et al. (24 | |