Computational study of short-range interactions in bacteriochlorophyll aggregates

Computational and Theoretical Chemistry 998 (2012) 87–97

Contents lists available at SciVerse ScienceDirect

Computational and Theoretical Chemistry journal homepage: www.elsevier.com/locate/comptc

Computational study of short-range interactions in bacteriochlorophyll aggregates J. Alster a, M. Kabelácˇ b,c, R. Tuma d, J. Pšencˇík a, J.V. Burda a,⇑ a

Department of Chemical Physics and Optics, Faculty of Mathematics and Physics, Charles University, Ke Karlovu 3, 121 16 Prague 2, Czech Republic Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, Flemingovo nám. 2, 166 10 Prague 6, Czech Republic c Faculty of Science, University of South Bohemia, Branišovská 31, 370 05 Cˇeské Budeˇjovice, Czech Republic d Astbury Centre for Structural Molecular Biology, Faculty of Biological Sciences, University of Leeds, Woodhouse Lane, LS2 9JT Leeds, United Kingdom b

a r t i c l e

i n f o

Article history: Received 30 May 2012 Received in revised form 25 June 2012 Accepted 2 July 2012 Available online 17 July 2012 Keywords: Chlorosome Bacteriochlorophyll Molecular simulations DFT Ab initio

a b s t r a c t Chlorosomes are light harvesting complexes, interior of which is composed of bacteriochlorophylls that are assembled into lamellar aggregates. In this work, short-range structural parameters of aggregated bacteriochlorophyll-c were explored. The knowledge of bonds involved in aggregation together with other experimental results provides constraints which are satisfied by eight possible structural motifs. Each motif is based on a dimer unit with esterifying alcohols extending to opposite sides of the lamellar layer. The possible models include two of the previously proposed structural models, so-called antiparallel piggy-back dimer [1,2] and a parallel syn-anti model [3]. Structural models of aggregates were built as a single layer consisting of 8  10 monomers and were optimized by Amber program using modified General Amber Force Field. From the optimized models the central tetramer unit was selected for singlepoint energy analyses. Three electronic methods: RI-MP2/SVP, B97-D/TZVP, and PM6-DH2 were chosen for revealing energy relations in the structural models. Calculations disclose that the most preferable structure is the antiparallel chain model corresponding to the antiparallel piggy-back dimer, followed by one of antiparallel sheet models and one of parallel sheet arrangements. All simulations showed that due to the crystal packing a large deformation of monomeric building block occurred. The deformation energies (the difference between the energy of optimized isolated monomer and that of the monomer in assembled structures) are about 15 kcal/mol for all considered computational methods. From energy decomposition, contributions of individual interactions to the overall stabilization energy were estimated. In addition, powder diffraction and optical spectra were calculated for all models and compared with experimental results. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction An extra-membranous light-harvesting antenna of certain phototrophic bacteria (including members of phyla Chlorobi, Chloroflexi and Acidobacteria), so-called chlorosome, is one of the most efficient light-harvesting systems in nature. Structure of the chlorosomes is based on self-organized aggregates of bacteriochlorophyll (BChl) molecules. The arrangement of these aggregates is unique in every chlorosome, which prevents solving their structure using high-resolution crystallography techniques. Nevertheless, recent advances in experimental and computational approaches (see below) revealed many features of the chlorosome structure. A combination of electron cryomicroscopy and small angle X-ray scattering (SAXS) has led to a description of the lamellar nature of BChl aggregates in the chlorosomes [4–7]. X-ray scattering is capable of disclosing key average parameters of the lamellar arrangement in chlorosome ensembles. On the other hand ⇑ Corresponding author. E-mail address: [email protected] (J.V. Burda). 2210-271X/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.comptc.2012.07.001

electron microscopy was instrumental in revealing lamellar arrangement in individual chlorosomes. It demonstrated that curvature is an essential feature of the lamellar layers [5]. The lamellar arrangement ranges from concentric cylinders that were observed in chlorosomes from a triple mutant of Cba. tepidum (bchQRU) to mostly disordered curved lamellae in wild type chlorosomes [5]. Aggregation of chlorosomal BChls (BChl-c, -d or -e) is mediated by their unusual chemical structure, which supports specific short-range interactions between chlorin rings (cf. Fig. 1, for a review see work of Blankenship and Matsuura [8] or Frigaard and Bryant [9]). BChl aggregates are formed by coordination of the central Mg cation of one BChl molecule to a hydroxy substituent at C31 of a neighboring BChl, which can also be H-bonded to the keto group at C131 of a third BChl molecule. Although this scheme is generally accepted, the presence of the H-bond was recently questioned and it was proposed that the keto group at C131 of the second BChl is weakly coordinated to the Mg ion of yet another BChl [10]. In addition, BChl-c, -d and -e lack methoxycarbonyl group at C132, which is common for other chlorophyll molecules and which would cause steric hindrance of the C131 keto and thus prevent aggregation [11].

88

J. Alster et al. / Computational and Theoretical Chemistry 998 (2012) 87–97

Fig. 1. Molecular structure of BChl-c homolog used in this work, 31R-[E,M] BChl-cF. The abbreviation 31R means that the chiral center at C31 is in R configuration, letters in a brackets stand for ethyl and methyl as the substituents at C8 and C12, respectively, and F in the subscript denotes that the molecule is esterified with farnesyl at C173. The red circles show the two stereo centers and the chiral C17 used to define the used nomenclature.

However, there still exist many possibilities how to arrange BChl molecules into aggregates while fulfilling the short-range interactions. Thus, further constraints and modeling are needed to explain the molecular arrangements. Many such models were based on NMR-derived constraints [1–3,12–16]. Complementary constraints are provided by X-ray powder diffraction, which established a dimer as the building block of the aggregates, effectively excluding models based on monomers and models containing more than two molecules in the unit cell. In addition, the esterifying alcohols of the monomers within the dimer should extend to opposite sides of the lamellar layer to enable formation of lamellar structures in aqueous solutions. The motivation for this study is to find plausible structures based on the energy considerations and compare them with the above-mentioned constraints arising from X-ray scattering and self-assembly experiments. Here we show that these requirements are fulfilled by eight possible structural motifs, including two of the previously proposed models: antiparallel piggy-back dimer [1,2] and syn-anti parallel dimer [3]. It is likely that different structural patterns of BChl-c molecules are possible in different species, or even within a single chlorosome, and therefore comprehensive survey of plausible structures is important. A natural way would be to use also constraints arising from the above-mentioned NMR studies, however, this proved difficult since different groups came to different conclusions concerning the proximity of various chlorin ring substituents in the aggregates (compare e.g. [3,15,17]). 2. Model building Structures were built from 31R-[E,M] BChl-cF (see Fig. 1 for explanation of the abbreviation). This homolog was selected

because it represents pigments found in model organisms Chlorobaculum (Cba.) tepidum and Chloroflexus aurantiacus. This homolog differs only by its esterifying alcohol from the single homolog found in Chloroflexus aurantiacus, it is one of the four homologs found in Cba. tepidum (although a minor one) and differs only in one methyl at C20 from 31R-[E,M] BChl-dF found in the triple bchQRU mutant, which exhibits exceptionally well ordered aggregates [5]. BChl-c molecules were assembled into aggregates according to several aggregation motifs (schematically drawn in Fig. 2). A motif is understood here as a pattern of intermolecular interactions and mutual orientations of the BChl molecules. A structure (molecular model) is a concrete realization of the pattern described by the given motif. The goal of this work was to find the most stable structure for each motif. These motifs were constructed according to the following rules: (1) the structure is based on dimers to fulfill constraints from X-ray scattering; (2) the principal interactions between BChl-c molecules is coordination of Mg cations by the hydroxyl groups. In addition the molecules were arranged in a way that distances between hydroxyl and carbonyl groups enabled formation of hydrogen bonds; (3) the structure is periodically extensible to form a layer, with esterifying alcohols alternatingly protruding out of the layer on both sides, in accordance with self-assembly in aqueous solutions, which is driven by the hydrophobic effect [18,19]. These rules/constraints allow to build eight possible motifs (Fig. 2), including several which have not been considered before. To fully describe the relative orientations of all stereocenters, it is necessary to define also the orientation of the OH group (which can rotate about the C3–C31 single bond) and Mg cation with respect to the C17-propionic acid moiety. The C17 chiral carbon is stereochemically pure and thus it provides a suitable reference to define the orientations of the two groups. However, different naming conventions were used in literature to describe these orientations, making a direct comparison difficult [3,14,20,21]. To avoid ambiguity, we adopt a nomenclature according to IUPAC rules: the term ‘‘Syn’’ (abbreviated to S) is used when the hydroxy group at C31 is on the same side of the chlorin ring as the C17-propionic acid moiety and ‘‘Anti’’ (abbreviated to A) when it is on the opposite side. To describe the orientation of the Mg cation with respect to the C17, we adopted nomenclature proposed by Balaban et al. [21]: The configuration with the magnesium cation on the opposite face than the C17propionic acid moiety is designated as a-configuration (abbreviated to a), b-configuration has the magnesium cation on the same face. This nomenclature also agrees with IUPAC nomenclature rules for tetrapyrroles [21]. Specifying the positions of both entities for both molecules forming the basic dimer gives a four-letter code denoting stereochemistry of each motif. Two kinds of structural motifs can be constructed, namely sheet structures and chain structures. Sheet structures are characterized by interconnecting Mg atoms with C31 hydroxyl group throughout the whole stacking BChl-c layer, e.g. from an oxygen of the lower BChl-c to the Mg atom of the BChl-c molecule situated above (cf. lower part of Fig. 2). The chain structures are formed by paired rows of BChl molecules. The Mg cations are coordinated by C31 oxygen atoms within the chain, while only stacking (p–p) interactions occur between the chains (cf. upper part of Figs. 2 and 3). If the mutual orientation of the dipole moments within the dimer is considered (parallel versus antiparallel), the eight plausible motifs can be finally divided into four groups: (a) parallel chains: aAaA, bSbS, (b) antiparallel chains: aAbS (corresponding to antiparallel piggy-back dimer [1,2]) (c) parallel sheets: aSbA (corresponding to a syn-anti model [3]), aSbA-d (‘‘d’’ for displaced), aSbA-dd (‘‘dd’’ for double displaced) and d) antiparallel sheets: aSaS, bAbA (Fig. 2).

J. Alster et al. / Computational and Theoretical Chemistry 998 (2012) 87–97

89

Fig. 2. Structural motifs of BChl-c as they were considered in aggregate models. A shaded rectangle means that the molecule is in a configuration with C17 carbon protruding forwards while a white empty rectangle denotes molecule in an opposite orientation. The side of the molecule facing out of the paper in Fig. 1 is marked with a thick black line. See text for nomenclature. The numbers 1, 2, 3, and 4 in the first parallel chain motif are labels of the monomer positions used in energy decomposition analysis.

3. Computational details

3.1. Structure optimizations

In this study we explore a short-range BChl organization in aggregates by a combination of several types of calculations. The first type is based on molecular mechanics (MM) simulations of larger molecular models compiled from structural motifs using Amber 8 suite of programs [22]. The second type concerns the interaction energy evaluation of a smaller central tetrameric unit of BChl-c performed by quantum mechanics (QM) methods at several computational levels. The third approach was used for predicting diffraction and the last type estimated optical spectra.

In order to be able to construct the structural motifs, the monomeric building block, the BChl-c molecule, was pre-optimized. For this step, a recommended quantum chemical method HF/6-31G(d) was employed with subsequent determination of Merz-Kolman electrostatic potential by Gaussian 03. Partial charges required for the empirical force field were determined by Restrained ElectroStatic Potential (RESP) procedure using the Antechamber program from the Amber package. General Amber Force Field (GAFF [23]) was used for all energy optimizations at the MM level.

90

J. Alster et al. / Computational and Theoretical Chemistry 998 (2012) 87–97

Fig. 3. The top panel illustrates the offset (D) of chain structures. The two possible selections of tetramer are also depicted:‘‘within-chain’’ (dashed green line) and ‘‘betweenchains’’ (dash-dotted blue line). The bottom panel shows the whole optimized molecular structure and illustrates the deformations on the edges. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The central magnesium cation plays an essential role in the aggregate formation. However, the non-bonding parameters for Mg interaction in the chlorin ring are not provided in the standard force-field. Therefore they had to be estimated from the results of QM calculations [24]. The used parameters are in general agreement with DFT calculations at the B3LYP/6-31G(d) level (Fig. S1). The individual models are represented by a single layer of aggregated BChl-c molecules. At first, the influence of slightly shifted and/or rotated monomers was tested for each motif on several structures encompassing a smaller layer of 4  6 BChl-c molecules. The structures were optimized by sander code from Amber package. The most stable conformations were selected and larger (8  10) structures were created and optimized. For chain structures it was necessary to determine the optimal offset between the chains, which is not set by OH group coordination with Mg. Structures with about twenty different offsets were prepared for each chain model (Fig. 3), optimized by molecular mechanics and the resulting offsets and corresponding energies of the optimized structures were compared. The optimal offset, corresponding to approximately half-molecule shift between the two neighboring chains, agrees well with the optimum geometry reported previously for porphyrin molecules from a p–p interaction model [25]. 3.2. Energy calculations From each MM-optimized structure, a central ‘tetramer’ unit of BChl-c was extracted in order to perform single-point (SP) energy calculations. The idea behind this approach was to get a consistently optimized tetramer structure for which interactions of surrounding BChl-c molecules are included (and not just a tetramer exposed to an implicit solvent or vacuum). On the other hand, optimization of structures larger than four BChl-c molecules is not feasible by present QM methods. Therefore the results of MM simulations were used as a source of optimized geometries for QM calculations. Three different QM computational levels were used for SP calculations. Since dispersion energy is of crucial

importance (besides the mostly electrostatic interaction of the C31 oxygen atoms with Mg2+) the comparison of the RI-MP2/def2-SVP [26] (employing Turbomole program 6.3 [27]), and the DFT functional with a dispersion term included B97-D/TZVP levels (as implemented Gaussian 09 program) was conducted. For ‘completeness’ SP energies were estimated by semiempirical PM6-DH2 method using the MOPAC2009 program [28]. At the RI-MP2 level the effect of basis set superposition error was corrected by Boys and Bernardi approach [29] since it is known that it plays an important role for double-zeta bases quality. The stabilization energies of the tetramer structures were calcuP lated according to equation: Estab ¼ Etetramer  i Emonomer . In the case i of chain motifs two different kinds of central tetramers had to be taken into account. The ‘within-chain’ tetramer contains three or four Mg  O–H coordinations in addition to the p–p stacking while the other tetramer unit, designated here as ‘between-chains’, is chosen so that it represents only between-chains stacking interactions (Fig. 3). Therefore the total stabilization energy of chain motifs is determined as an average of both tetramer units. Furthermore, deformation energies were computed as the difference between energies of frozen structure of a given monomer in the tetramer and the optimized structure of an isolated BChl-c. Also, energy decompositions were performed for all the tetramers in order to reveal the importance of individual bonding contributions, namely the role of Mg  O–H coordination, stacking interaction, and H-bond terms. Therefore, pair-pair energies were calculated according to an equation similar to the stabilization enP monomer ergy formula, the only difference was that instead of i Ei summation the energies of two dimers were used. 3.3. Powder diffraction To compare the generated models with experimental data, powder diffraction was calculated for structures corresponding to all motifs. As discussed below, a curvature was applied to the models in order to improve the agreement with experimental data. Since

91

J. Alster et al. / Computational and Theoretical Chemistry 998 (2012) 87–97

Table 1 Geometric parameters of optimized structures for all studied motifs – unit cell parameters (lengths of vectors and angle between them), two largest distances of crystallographic planes and mass density of the resulting aggregate. The density was calculated using the unit cell parameters and the lamellar spacing of 2.1 nm. a (nm)

b (nm)

Gamma (deg)

Distances of crystalographic planes (nm)

Mass density (g/cm3)

bSbS aAbS

0.74(5) 0.760(3) 0.740(8)

1.43(3) 1.438(4) 1.459(5)

81(4) 94(5) 85(1)

1.42 1.43 1.45

0.7 0.76 0.74

1.24 1.15 1.17

aSbA aSbA-d aSbA-dd

0.79(2) 0.85(3) 0.706(5)

1.37(2) 1.34(1) 2.36(2)

80(1) 78(1) 42.0(8)

1.33 1.31 1.57

0.75 0.82 0.59

1.26 1.14 1.13

aSaS bAbA

0.728(8) 0.761(5)

1.466(2) 1.447(8)

90.6(6) 90.8(6)

1.47 1.45

0.73 0.76

1.16 1.14

aAaA

Table 2 Stabilization energies (kcal/mol) obtained by RI-MP2/def2-SVP, B97-D/TZVP, PM6-DH2, and Amber 8 for all the explored structures. The numbers in italics indicate ranking of the three structures with the largest stabilization energies obtained by each method. The two lines for each chain structure contain values for ‘within chain’ (upper) and ‘between chain’ (lower) tetramers, respectively.

aAaA bSbS aAbS aSbA aSbA-d aSbA-dd aSaS bAbA

RI-MP2

B97-D

PM6-DH2

Amber

80.2 80.8 103.3 78.1 90.8 100.6 86.7 97.7

108.7 107.1 135.6 108.8 118.1 129.2 118.9 131.7

88.8 87.2 110.5 81.9 92.0 99.5 97.5 114.9

106.4 106.3 130.4 110.2 122.2 141.9 122.1 130.0

1

2 3

1

3 2

2

3 1

performed. Powder diffraction was computed from the coordinates generated for BChl-c molecules in one layer and bacteriochlorophyllide-c molecules (BChl-c without esterifying alcohol) in the three-layer lamella using Debye scattering equation [30] with corrections for angular dependence of atomic scattering factors [31]. 3.4. Optical spectra

2

It is well known that optical spectrum of chlorophylls exhibits Q and Soret bands. For the reproduction of the monomeric chlorophyll spectrum it is necessary to determine at least 5 excited states. This means that for tetramer structure it is necessary to determine ca (4  5) 20 lowest excited states. So many states for a system of such extent do not allow employing of any high level computational methods. Therefore electronic excitation levels were estimated by semiempirical ZINDO method [32]. Based on our previous experience with calculations of chlorophyll optical spectra [33] we concentrated on the relative difference between tetramer-monomer energy levels rather than on absolute energies of the electron transitions for both the Qy and Soret band. These two ‘spectral shifts’ can be compared with analogous shifts

1 3

Cba. tepidum lamellae are wrapped around with apparent radii of curvature between 2 and 4 nm [5] the planar BChl layers were transformed into a curved lamella consisting of three curved layers (curvature 2, 3, 4 nm), separated by a typical lamellar spacing of 2 nm. Steric clashes resulting from the curvature were partially alleviated by dilation along the arch. No further optimization was

Table 3 Pairing interaction and stabilization energies (kcal/mol). The two lines for each chain structure contain values for ’within chain’ (upper) and ’between chain’ (lower) tetramers, respectively. RI-MP2/def2-SVP

B97-D/TZVP

12–34

13–24

14–23

Stab

12–34

13–24

14–23

Stab

40.5 15.8 47.5 20.5

97.7 49.7 114.2 26.0

73.7 46.3 84.4 34.1

104.1 56.3 121.0 40.5

53.6 29.8 60.5 34.6

117.0 70.9 133.1 53.0

94.7 70.2 103.8 44.7

131.0 86.4 147.0 67.1

aAbS

34.6 15.4

134.1 56.0

126.2 54.4

135.1 47.6

49.0 27.8

157.7 85.4

152.2 77.1

162.0 77.4

aSbA aSbA-d aSbA-dd

24.6 21.9 22.2

70.2 81.3 91.1

63.2 79.9 89.6

143.8 62.8 78.1

38.5 35.1 36.8

92.8 101.5 112.2

87.7 100.9 110.5

175.8 95.3 108.8

aSaS

22.8 23.2

80.7 88.4

72.9 87.0

90.8 86.7

35.3 36.7

104.7 115.6

100.3 114.0

118.1 118.9

12–34

13–24

14–23

Stab

12–34

13–24

14–23

Stab

aAaA

45.4 29.4

86.5 63.9

77.2 55.5

102.5 75.1

53.5 31.4

114.8 68.4

95.9 61.6

132.1 80.7

bSbS

54.2 33.6

86.8 39.4

77.8 39.3

117.8 56.6

63.7 37.0

131.6 43.2

107.2 42.5

151.2 61.3

aAbS

40.1 23.1

127.6 72.1

126.5 60.3

133.8 61.4

44.4 30.0

150.2 80.6

150.3 66.0

159.5 67.6

aSbA aSbA-d aSbA-dd aSaS

31.7 28.7 33.5 31.2 33.7

66.4 80.5 86.2 86.2 101.9

66.9 76.0 80.9 80.7 97.7

143.5 77.5 81.9 92.0 97.5

40.0 39.7 37.2 36.1 38.1

91.8 102.4 125.4 104.6 113.5

88.6 102.4 121.3 103.5 108.3

172.4 88.3 110.2 122.2 122.1

aAaA bSbS

bAbA

PM6-DH2

bAbA

MM/GAFF

92

J. Alster et al. / Computational and Theoretical Chemistry 998 (2012) 87–97

Fig. 4. Central parts of optimized molecular structures for selected motifs: (a) antiparallel chain aAbS, (b) antiparallel sheet bAbA, (c) parallel sheet aSbA-dd. Esterifying alcohol residues and non-polar hydrogens were omitted for clarity.

obtained from experimental spectra of monomer, dimer, and aggregated sample of BChl-c. 4. Results 4.1. Structural parameters Geometry parameters were determined from the 8  10 structural models optimized by MM. The structural parameters were

averaged over all non-rim molecules in each structure and compared with available experimental data. One significant result of MM calculations is a large deformation of the chlorin ring in the aggregate compared to the pre-optimized monomer and a significant displacement of Mg atom out of the plane of the chlorin ring. The displacement is between 0.03 and 0.05 nm and is larger for b-configurations (displacements for individual structures are summarized in Table S1). On the other hand, the arrangement of esterifying alcohols remains almost unaffected.

J. Alster et al. / Computational and Theoretical Chemistry 998 (2012) 87–97

While the Mg  O–H distances in optimized structures always correspond to that of coordination, the H-bond = O  H–O distances were occasionally larger than that required for the H-bond. This is in agreement with much lower energy of the H-bond compared to coordination (see ‘Energy analyses’). Spacing between individual molecules is similar for all structures, but some differences can be found. Particularly, perpendicular distance between chlorin rings in chain structures is smaller within the chain than in-between the chains. On the other hand, the distance is fairly uniform in sheet structures (distances for individual structures are summarized in Table S2). The important geometry parameters of individual structures are summarized in Table 1. Distances between crystallographic planes (Bragg spacing) of 1.17 nm and 0.94 nm (corresponding q values 6.7 and 5.4 nm1) were observed in X-ray scattering of chlorosomes from Cba. tepidum [4] and assigned to the two-dimensional lattice of the layer [6]. Corresponding distances can be calculated for the optimized structures from unit cell parameters. These parameters were determined by averaging translation vectors between equivalent BChl molecules with the Mg atoms taken as the lattice points. The two largest Bragg spacings for each structure and the corresponding unit cell parameters are summarized in Table 1. It can be seen that none of the explored structures exhibits particularly good agreement with the experimental values. This might be due to the intrinsic curvature of the lattice in the chlorosome as discussed below. The density of the chlorosome was determined to be approx. 1.16 g/cm [34]. Considering that proteins account for about 30% of the chlorosome dry weight [8] (protein density is 1.3 g/cm3), it can be estimated that the density of pigments in the chlorosome is about 1.11 g/cm3. This number can be compared with densities estimated for the molecular structures (Table 1). It can be seen that energetically favorable structures (aAbS, bAbA and aSbA-dd, cf. bellow) are among those, which are getting closer to the expected density of pigments in the chlorosome. The lower experimental density can be caused by the disorder present in chlorosomes.

93

bonding interactions contribute to the final stabilization, namely H-bond, Mg  O–H coordination, and stacking. Individual pair–pair energies are listed in Table 3. In the first column (labeled as 12–34, for the monomer numbering see Fig. 2) association energies of two stacked pairs are collected. These values reflect two different situations (a) interaction of two H-bonds + 1/2 chlorin:chlorin stacking in the case of the chain motifs, where pairs of chlorophylls belong to different chains (‘between-chains’), antiparallel ‘within-chain’ tetramers, or in the case of the sheet motifs and (b) interaction of two H-bonds + 1/2 stacking energy + one Mg  O–H coordination in the ‘within-chain’ parallel-chain tetramers. From these values it follows that the Mg  O–H coordination of the ‘withinchain’ structures represents energy of about 20–30 kcal/mol, which can be roughly obtained by subtracting the corresponding values for ‘within-chain’ and ‘between-chains’ (in good agreement by all four methods). Larger values of 12–34 pairing were obtained by MM technique (ca 40 kcal/mol) while MP2 predicts this energy to be only about 20 kcal/mol. The DFT and PM6 methods predict these values in the range of 28–38 kcal/mol. Nevertheless, all methods make the clear energy distinction between sheet and ‘between-chains’ tetramers, e.g. in DFT predicts the 12–34 paring energy above 35 kcal/mol for the sheet structures whereas for chain structures it gives lower values. The latter energies reflect the fact

4.2. Energy analyses The central tetramer unit was selected from the most stable arrangement of each motif and used for further analysis by SP calculations performed at the RI-MP2, B97-D, and PM6 level of theory. The energy characteristics obtained by these methods were compared with the MM energies (Table 2 and 3, and S3). Table 2 has four sections – the first lists stabilization energies at the RI-MP2 level together with labels of the first three most stable structures. The second part summarizes results from the B97-D level followed by PM6 energies. In the last section, results from the original MM are displayed. Comparing the stabilization energies of the considered motifs reveals that the most preferable structure is the antiparallel chain aAbS (at the MP2 and DFT levels) followed by the antiparallel sheet bAbA together with parallel sheet structure aSbA-dd. These structures were found to be the most preferable by all four methods. These motifs are separated by an energy gap of about 10– 20 kcal/mol (depending on the method used) from other motifs. Besides stabilization energies the deformation energies were also determined. Complete set of these energies for all methods is collected in supplementary Table S3. These values clearly demonstrate that large deformation of the monomeric building blocks occurs in the MM optimized monolayer models. The averaged deformation energies are around 15 kcal/mol for all the employed levels, indicating substantial deviations from the geometry of isolated monomer due to the packing effects in aggregates. The energy decomposition gives further insight into the origin of stabilization energy. As mentioned above three important

Fig. 5. (a) Calculated diffraction for a single layer of aAbS aggregate composed of BChl-c molecules. The figure demonstrates that increasing the curvature (i.e. decreasing radius) leads to improved agreement between the observed (vertical red lines) and predicted peaks. The peaks between 4 and 5 nm1 are due to scattering from esterifying alcohols, which remained ordered even after optimization. (b) Calculated diffraction for three layers of aggregates composed of bacteriochlorophyllide c molecules. The vertical red lines denote positions of experimentally observed peaks. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

94

J. Alster et al. / Computational and Theoretical Chemistry 998 (2012) 87–97

Fig. 6a. Calculated optical spectra of monomer (bottom panel) and tetramer structures (aAaA, bSbS, aAbS, aSbA, aSbA-d, aSbA-dd, aSaS, and bAbA from second bottom panel to top) of BChl-c. The expected position of Qy and Soret bands are marked by vertical dashed line.

that perpendicular distance between chains is larger than the analogous distance for ‘within-chain’ or sheet tetramers (cf. discussion of structural parameters and Table 3). For estimation of the Mg  O–H coordination energies also another consideration can be used: the difference between the stabilization energies of the ‘within-chain’ and ‘between-chains’ structures, respectively. In accord with the estimation from 12–34 pair-pair energies the values are about 24 kcal/mol (slightly less for MP6 (ca 19 kcal/mol) and slightly more for MM methods (ca 26 kcal/mol)). In reality these numbers should be a little bit lower since the stacking interaction

is underestimated in ‘between-chains’ due to a slightly larger perpendicular distance, which means that the actual energies should be reduced by ca 4 kcal/mol. In the next column (13–24), association of two H-bonded pairs is calculated. In this case the breaking of all stacking interactions and Mg  O–H coordinations occurs. Depending on the motif there are 4 (antiparallel ‘in-chain’), 3 (parallel ‘in-chain’), 2 (sheet) or 0 (‘between-chains’) Mg  O–H coordinations. This pair-pair energy attains the largest values among the whole energy-decomposition scheme and it can be directly used for determination of H-bonding

J. Alster et al. / Computational and Theoretical Chemistry 998 (2012) 87–97

energy. Since H-bonding is the only interaction, which is not included in the 13–24 energy (both H-bonds survive the separation of 13. . .24 pairs to infinity) the difference between stabilization energy and 13–24 pairing gives the required estimation of H-bond energy: 4–9 kcal/mol/1H-bond (4 kcal/mol at MP2, 8 kcal/mol at DFT and 9 kcal/mol at the PM6 and MM level). Last pairing energy 14–23 is very similar to 13–24 as to the individual values. From Fig. 2 it can be noticed that both H-bonds are broken together with two Mg  O–H coordinations and 2  (1/2) of stacking interaction. It means that only 1/2 of stacking interaction remains included in one of the parts after separation of both parts of tetramer to infinity. In analogous manner to previous 13–24 pair energy, the difference between stabilization energy and 14–23 pair energy reveals 1/2 of stacking interaction of the sheet and antiparallel chain motifs (12 kcal/mol at MP2 and PM6, 19 kcal/mol at DFT and 20 kcal/mol at MM level). Similar values can be obtained for the parallel ‘between-chains’ structures. In the case of parallel ‘in-chain’ 14–23 pair energy, two of three Mg  O–H coordinations are broken so that one still remains after separation of both parts. Based on the energy decomposition it can be concluded that the average energy of H-bond can be estimated to be approx. 8 kcal/mol (for comparison this is higher than that obtained for H-bonding between water molecules), energy of Mg  O–H coordination contributes about 20 kcal/mol. Finally the p–p stacking interaction is about 40 kcal/mol for the completely stacked chlorin ring. 4.3. Diffraction calculations The curvature is an intrinsic property of the BChl aggregates in chlorosomes, however, the modeled structures remained planar. Since it is not clear what is the origin of the curvature, it will be necessary to use molecular dynamics on a larger, multilayer assembly to solve this problem. Such demanding simulations are beyond the scope of this contribution. Here we want to illustrate the importance of curvature for obtaining a better agreement of the modeled structures with experimental results. For that we simulate X-ray scattering from few curved configurations. As the optimized structures do not yield diffraction similar to the experiment, we modified the structures by applying curvature by wrapping them around the x-axis (see Fig. 3 for definition of the x axis), which constitutes the long axis of chlorosomes [6]. Fig. 5a shows the scattering pattern calculated for a single layer of aAbS aggregate (the structure with the best energy score), modified with an increasing curvature (decreasing radius) applied to the optimized structure. The figure demonstrates that the increasing curvature leads to the improvement of the agreement with the experimental results. The diffraction close to the experimentally observed peak at 5.2 nm1 remains at the same position because the corresponding distance (along the x-axis) is not affected by the applied curvature. However, the position of the peak between 6 and 10 nm1 moves to smaller q values as the applied curvature and dilation increase the distance along y-axis. The scattering at 4.5 nm1 (Fig. 5a) is due to the esterifying alcohol tails, as could be proved by comparison of the scattering pattern from the same layer consisting of BChl-c molecules (containing the tails) and bacteriochlorophyllide-c (after the tail removal). The esterifying alcohols remain well ordered in the models after optimization (Fig. 4), probably due to the lack of interaction with adjacent layers, which is a condition unlikely to be fulfilled in the aggregates within chlorosomes. In addition, interdigitation of the alcohols from neighboring layers will lead to further redistribution of the relative distances between the tails towards shorter values. Thus, it is unlikely that this diffraction is observed for chlorosomes. The diffractions presented in Fig. 5a were calculated for an aggregate consisting of a single layer of BChl-c molecules.

95

Therefore, it does not exhibit the experimentally observed lamellar peak corresponding to the distance between the neighboring layers of BChl aggregates. A scattering peak corresponding to this distance was observed around 3 nm1 (spacing 2.1 nm) in chlorosomes from Cba. tepidum. The lamellar peak is nicely reproduced by calculated diffraction from three layers. The three layers used to obtain diffraction patterns presented in Fig. 5b have curvatures with a radius of 2, 3 and 4 nm, respectively, and the distances between chlorin planes in the optimized structures were dilated by 20% in the direction of the arch, to avoid steric clashes at the inner side of the arch. To avoid the presence of the peak caused by the artificially ordered esterifying alcohols, these were removed from the structures. It can be seen that the structures with the lowest energies (aAbS, bAbA, aSbA-dd) exhibit also better agreement with the experimental results. For instance, the position of the peak observed between 6 and 10 nm1 is closest to the experimentally observed 6.7 nm1 for these structures (together with aAaA). These structures also exhibit a peak close to the position of experimentally observed 5.4 nm1, however, this peak is significantly weaker in the calculated diffractions, making the comparison with experimental results difficult. This can be caused by the fact that the relative intensities and widths depend on disorder, which was not taken into account and may differ for individual diffractions. It is interesting to note that the experimentally observed peak at 14 nm1 also appears after formation of the lamellar structure consisting of three layers. This suggests that in order to reproduce the diffraction pattern a larger, multilayer model needs to be considered. 4.4. Comparison of electron spectra Electron spectra were determined for all the explored tetramer structures using semiepmirical ZINDO method for 20 lowest excited states. This corresponds to 4  5 energy levels of BChl-c monomer. The monomer spectrum is presented at the bottom of Fig. 6. Since the error in determination of individual electronic levels is relatively large we decided to examine only the differences between corresponding tetramer-monomer energy levels of the Qy and Soret band. These two ‘spectral shifts’ were compared with analogous shifts in experimental spectra between monomer and dimer or aggregated BChl-c as displayed in Fig. 7. The spectral shift caused by assembling BChl-c into aggregated form is about 80 nm (40 nm for dimer) for the Qy band while for the Soret band it

Fig. 6b. Absorption spectra of BChl-c in various aggregation states: monomeric (in ethanol, blue solid line), dimeric (in aqueous buffer, green dotted line) and aggregated (mixed with lecithin, in aqueous buffer, red dashed line). Maxima of Qy bands are shown.

96

J. Alster et al. / Computational and Theoretical Chemistry 998 (2012) 87–97

amounts to 30 nm. Both expected shifts are marked in Fig. 6 by vertical lines at 450 and 750 nm. Comparing the examined tetrameric structures, the largest shift of the Qy spectral line was obtained for aSbA-dd structure (ca 60 nm). Spectrum of aAaA tetramer displayed the second largest shift of 40 nm. The expected spectral shift of the Soret band is relatively well reproduced by three tetrameric structures aAaA, bAbA, and aSbA-dd. The relatively good agreement between calculated and experimental spectra lends further support to the proposal that the energetically preferred structures represent the realistic configurations that may be found in native chlorosomes.

and JVB is grateful to Kontakt Project No. ME 10149 of Ministry of Education. Authors are grateful to generous access to computational resources from Supercomputation Metacenters in Prague and Brno.

5. Conclusions

[1] T. Nozawa, K. Ohtomo, M. Suzuki, H. Nakagawa, Y. Shikama, H. Konami, Z.Y. Wang, Structures of chlorosomes and aggregated BChl c in chlorobium tepidum from solid state high resolution CP/MAS 13C NMR, Photosynth. Res. 41 (1994) 211–223. [2] A. Egawa, T. Fujiwara, T. Mizoguchi, Y. Kakitani, Y. Koyama, H. Akutsu, Structure of the light-harvesting bacteriochlorophyll c assembly in chlorosomes from chlorobium limicola determined by solid-state NMR, Proc. Natl. Acad. Sci. 104 (2007) 790–795. [3] S. Ganapathy, G.T. Oostergetel, P.K. Wawrzyniak, M. Reus, A.G.M. Chew, F. Buda, E.J. Boekema, D.A. Bryant, A.R. Holzwarth, H.J.M. de Groot, Alternating syn-anti bacteriochlorophylls form concentric helical nanotubes in chlorosomes, Proc. Natl. Acad. Sci. 106 (2009) 8525–8530. [4] J. Psencik, T.P. Ikonen, P. Laurinmäki, M.C. Merckel, S.J. Butcher, R.E. Serimaa, R. Tuma, Lamellar organization of pigments in chlorosomes, the light harvesting complexes of green photosynthetic bacteria, Biophys. J. 87 (2004) 1165–1172. [5] G.T. Oostergetel, M. Reus, A. Gomez Maqueo Chew, D.A. Bryant, E.J. Boekema, A.R. Holzwarth, Long-range organization of bacteriochlorophyll in chlorosomes of chlorobium tepidum investigated by cryo-electron microscopy, FEBS Lett. 581 (2007) 5435–5439. [6] J. Psencik, M. Torkkeli, A. Zupcanova, F. Vacha, R.E. Serimaa, R. Tuma, The lamellar spacing in self-assembling bacteriochlorophyll aggregates is proportional to the length of the esterifying alcohol, Photosynth. Res. 104 (2010) 211–219. [7] G.T. Oostergetel, H. van Amerongen, E.J. Boekema, The chlorosome: a prototype for efficient light harvesting in photosynthesis, Photosynth. Res. 104 (2010) 245–255. [8] R.E. Blankenship, K. Matsuura, Antenna complexes from green photosynthetic bacteria, in: B.R. Green, W.W. Parson (Eds.), Light-Harvesting Antennas in Photosynthesis, Kluwer Academic Publishers, Dordrecht, 2003, pp. 195–217. [9] N.U. Frigaard, D.A. Bryant, Chlorosomes: antenna organelles in photosynthetic green bacteria, in: J.M. Shively (Ed.), Complex Intracellular Structures in Prokaryotes (series: Microbiology Monographs, vol. 2), Springer, Berlin, 2006, pp. 79–114. [10] T. Jochum, C.M. Reddy, A. Eichhofer, G. Buth, J. Szmytkowski, H. Kalt, D. Moss, T.S. Balaban, The supramolecular organization of self-assembling chlorosomal bacteriochlorophyll c, d, or e mimics, Proc. Natl. Acad. Sci. 105 (2008) 12736– 12741. [11] T. Oba, H. Tamiaki, Why do chlorosomal chlorophylls lack the C13(2)methoxycarbonyl moiety? An in vitro model study, Photosynth. Res. 61 (1999) 23–31. [12] A.R. Holzwarth, K. Schaffner, On the structure of bacteriochlorophyll molecular aggregates in the chlorosomes of green bacteria. A molecular modelling study, Photosynth. Res. 41 (1994) 225–233. [13] T. Mizoguchi, K. Hara, H. Nagae, Y. Koyama, Structural transformation among the aggregate forms of bacteriochlorophyll c as determined by electronicabsorption and NMR spectroscopies: dependence on the stereoisomeric configuration and on the bulkiness of the 8-C side chain, Photochem. Photobiol. 71 (2000) 596–609. [14] B.J. van Rossum, D.B. Steensgaard, F.M. Mulder, G.J. Boender, K. Schaffner, A.R. Holzwarth, H.J.M. de Groot, A refined model of the chlorosomal antennae of the green bacterium chlorobium tepidum from proton chemical shift constraints obtained with high-field 2-D and 3-D MAS NMR dipolar correlation spectroscopy, Biochemistry 40 (2001) 1587–1595. [15] Y. Kakitani, H. Nagae, T. Mizoguchi, A. Egawa, K. Akiba, T. Fujiwara, H. Akutsu, Y. Koyama, Assembly of a mixture of isomeric BChl c from chlorobium limicola as determined by intermolecular C-13-C-13 dipolar correlations: coexistence of dimer-based and pseudo-monomer-based stackings, Biochemistry 45 (2006) 7574–7585. [16] Y. Kakitani, Y. Koyama, Y. Shimoikeda, T. Nakai, H. Utsumi, T. Shimizu, H. Nagae, Stacking of bacteriochlorophyll c macrocycles in chlorosome from chlorobium limicola as revealed by intermolecular c-13 magnetic-dipole correlation, X-ray diffraction, and quadrupole coupling in Mg-25 NMR, Biochemistry 48 (2009) 74–86. [17] Z.Y. Wang, M. Umetsu, M. Kobayashi, T. Nozawa, Complete assignment of H1 NMR spectra and structural analysis of intact bacteriochlorophyll c dimer in solution, J. Phys. Chem. B 103 (1999) 3742–3753. [18] P. Klinger, J.B. Arellano, F.E. Vacha, J. Hala, J. Psencik, Effect of carotenoids and monogalactosyl diglyceride on bacteriochlorophyll c aggregates in aqueous buffer: implications for the self-assembly of chlorosomes, Photochem. Photobiol. 80 (2004) 572–578.

There exist eight possible arrangements of BChl molecules in the aggregate, which are in principal agreement with the known bonding patterns and other experimental findings. It is possible to build molecular structures from all of these eight motifs. This suggests that aggregates in chlorosomes may adopt different structural motifs, either in chlorosomes from different species or even present two distinct phases within the same chlorosome. The highest stabilization energies were obtained for the antiparallel piggy back dimer (aAbS), originally proposed on the basis of NMR results [1,2]. Another two structures are plausible, namely bAbA and aSbA-dd. The bAbA is an antiparallel sheet motif, which, to our best knowledge, has not been proposed previously. The two antiparallel models are also compatible with Stark spectroscopy, which revealed an absence of permanent dipole moment difference between the ground and the excited state. Such a result can be easily explained by the antiparallel orientation of the BChl molecules [34]. The aSbA-dd structure is a version of the parallel synanti model proposed by Ganapathy et al. [3]. Although the present approach does not allow to draw a decisive conclusion about the arrangement of aggregates in chlorosomes, it provides basis for further evaluation of the possible models. Despite the fact that individual structures differ by ‘‘only a few’’ kcal/mol, it should be noted that this value represents the energy difference in the tetramer unit, which will be dramatically multiplied by the large number of BChls in chlorosomes. However, it would be also fair to say that for making such statement more persuading we should perform much more extensive averaging not only over our 8  10 model but also in a better sampling of the phase space using e.g. extensive MD simulations. Such simulations exceed the scope of the present study and will be subject of our next explorations. From the energy decomposition the average energy of H-bond can be estimated as 8 kcal/mol, energy of Mg  O–H coordination contributes about 20 kcal/mol and the p–p stacking interaction adds about 40 kcal/mol for a completely stacked chlorin ring. One of the expected goals of the MM optimization was to examine whether the experimentally observed curvature of BChls layers can be reproduced. However, the optimized layers remained flat. Nevertheless, the importance of the curvature became evident from the diffraction calculations. These results show that it is essential to take into account the experimentally observed curvature in order to further improve the agreement between the models and the experimentally observed powder diffraction. However, at present the origin of curvature remains unknown and will be addressed in our future examination. Acknowledgements This work was supported by Czech Science Foundation, project 206/09/0375. MK is grateful to Grant Project No. IAA400550808 (Grant Agency of the Academy of Sciences of the Czech Republic)

Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.comptc.2012. 07.001. References

J. Alster et al. / Computational and Theoretical Chemistry 998 (2012) 87–97 [19] A. Zupcanova, J.B. Arellano, D. Bina, J. Kopecky, J. Psencik, F. Vacha, The length of esterifying alcohol affects the aggregation properties of chlorosomal bacteriochlorophylls, Photochem. Photobiol. 84 (2008) 1187–1194. [20] T. Oba, H. Tamiaki, Which side of the pi-macrocycle plane of (bacterio)chlorophylls is favored for binding of the fifth ligand?, Photosynth Res. 74 (2002) 1–10. [21] T.S. Balaban, P. Fromme, A.R. Holzwarth, N. Krauss, V.I. Prokhorenko, Relevance of the diastereotopic ligation of magnesium atoms of chlorophylls in photosystem I, BBA-Bioenergetics 1556 (2002) 197–207. [22] D.A. Case, T.E. Cheatham, T. Darden, H. Gohlke, R. Luo, K.M. Merz, A. Onufriev, C. Simmerling, B. Wang, R.J. Woods, The Amber biomolecular simulation programs, J. Comput. Chem. 26 (2005) 1668–1688. [23] J.M. Wang, R.M. Wolf, J.W. Caldwell, P.A. Kollman, D.A. Case, Development and testing of a general amber force field, J. Comput. Chem. 25 (2004) 1157–1174. [24] Z. Futera, private communications, unpublished work. [25] C.A. Hunter, J.K.M. Sanders, The nature of Pi–Pi interactions, J. Am. Chem. Soc. 112 (1990) 5525–5534. [26] F. Weigend, M. Haser, H. Patzelt, R. Ahlrichs, RI-MP2: optimized auxiliary basis sets and demonstration of efficiency, Chem. Phys. Lett 294 (1998) 143–152. [27] R. Ahlrichs, M. Bar, M. Haser, H. Horn, C. Kolmel, Electronic-structure calculations on workstation computers – the program system turbomole, Chem. Phys. Lett 162 (1989) 165–169.

97

[28] J.J.P. Stewart, Comparison of the accuracy of semiempirical and some DFT methods for predicting heats of formation, J. Mol. Model. 10 (2004) 6–12. [29] S.F. Boys, F. Bernardi, Calculation of small molecular interactions by differences of separate total energies – some procedures with reduced errors, Mol. Phys. 19 (1970) 553–566. [30] P. Debye, Zerstreuung von röntgenstrahlen, Ann. Physik 351 (1915) 809– 823. [31] D.T. Cromer, J.B. Mann, X-ray scattering factors computed from numerical hartree-fock wave functions, Acta Crystallogr. A-Cryst. Phys. Diffr. Theor. Gen. Crystallogr. A 24 (1968) 321–324. [32] M. Zerner, ZINDO – semiempirical quantum chemical program, in: K.B. Lipkowitz, D.B. Boyd (Eds.), Reviews in Computational Chemistry, VCH, New York, 1991. [33] Z. Vokacova, J.V. Burda, Computational study on spectral properties of the selected pigments from various photosystems: structure-transition energy relationship, J. Phys. Chem. A 111 (2007) 5864–5878. [34] R. Frese, U. Oberheide, I.H.M. van Stokkum, R. van Grondelle, M. Foidl, J. Oelze, H. van Amerongen, The organization of bacteriochlorophyll c in chlorosomes from chloroflexus aurantiacus and the structural role of carotenoids and protein – an absorption, linear dichroism, circular dichroism and stark spectroscopy study, Photosynth. Res. 54 (1997) 115– 126.