Contrasting Modes of Self-Assembly and Hydrogen-Bonding Heterogeneity in Chlorosomes of Chlorobaculum tepidum.

Chlorosome antennae form an interesting class of materials for studying the role of structural motifs and dynamics in nonadiabatic energy transfer. They perform robust and highly quantum-efficient transfer of excitonic energy while allowing for compositional variation and completely lacking the usual regulatory proteins. Here, we first cast the geometrical analysis for ideal tubular scaffolding models into a formal framework, to relate effective helical properties of the assembly structures to established characterization data for various types of chlorosomes. This analysis shows that helicity is uniquely defined for chlorosomes composed of bacteriochlorophyll (BChl) d and that three chiral angles are consistent with the nuclear magnetic resonance (NMR) and electron microscope data for BChl c, including two novel ones that are at variance with current interpretations of optical data based on perfect cylindrical symmetry. We use this information as a starting point for investigating dynamic and static heterogeneity at the molecular level by unconstrained molecular dynamics. We first identify a rotational degree of freedom, along the Mg-OH coordination bond, that alternates along the syn-anti stacks and underlies the (flexible) curvature on a larger scale. Because rotation directly relates to the formation or breaking of interstack hydrogen bonds of the O-H···O=C structural motif along the syn-anti stacks, we analyzed the relative fractions of hydrogen-bonded and the nonbonded regions, forming stripe domains in otherwise spectroscopically homogeneous curved slabs. The ratios 7:3 for BChl c and 9:1 for BChl d for the two distinct structural components agree well with the signal intensities determined by NMR. In addition, rotation with curvature-independent formation of stripe domains offers a viable explanation for the localization and dispersion of exciton states over two fractions, as observed in single chlorosome fluorescence decay studies.

Introduction

The largest light-harvesting antennae found in nature, the chlorosomes, have evolved at an early stage and allow photosynthetic green bacteria to survive under light-limiting conditions and at very low light intensity.[1] In contrast to other photosynthetic light-harvesting antennae, for which a protein scaffold supports the structural order and functionality, the chlorosome antenna function is accomplished by hundreds of thousands bacteriochlorophyll (BChl) c, d, or e molecules that self-assemble into pseudosymmetric syn–anti parallel stacks forming curved sheets and concentric tubes.[2] Although every chlorosome has a unique structure, they all share the common feature of long-range transfer of extended excitons over distances of several hundreds of nanometer.

Whereas high BChl concentrations in solution would result in energy dissipation and loss, owing to a rapid quenching of excitonic energy, nature has managed to avoid this in chlorosomes by a protein-free encoding of distinct, functional BChl packings that display remarkable homogeneity and have gross structural variability at the supramolecular and higher levels. Consequently, there is much interest in obtaining a fundamental understanding of structure–dynamics–function relations in chlorosomes, including how quasi-coherent exciton transfer driven by nuclear motion enables the near-unity yield. This principle can potentially be exploited for more efficient, sustainable, and environmentally clean ways of converting solar energy via the rational design of molecular components for artificial photosynthesis. Currently, it is unclear how the hierarchy of structural and dynamical properties enables an amazing quantum efficiency for transfer within and between tubes, and the dominant mechanisms are strongly debated.[3,4]

Persistent variability in the form of two components that are each well-ordered on a local scale has been concluded from the pronounced splitting of nuclear magnetic resonance (NMR) signals for several mutants.[5] A recent study of single chlorosomes confirms the NMR finding that the variations in the molecular arrangement and the microscopic disorder over the whole chlorosome assembly are limited.[6] The optical linewidths of ∼100 cm–1 that can be extracted for individual chlorosomes, comprising both static and dynamic disorder, confirm the remarkably homogeneous molecular arrangement in two components observed with MAS NMR. The four peaks in polarization-resolved fluorescence-excitation data were attributed to low- and high-energy doublets for separate structure motifs with slightly different characteristics, that is, curvature and molecular arrangement. Their attribution to quantitatively different tubes contrasts with the earlier MAS NMR and transmission electron microscopy (TEM) investigations, which alternatively attributed the two components to transition regions between domains within the tube. The single layer line and other details of the MAS NMR and TEM diffraction responses for different chlorosome species and mutants[5] provide converging evidence that tubes in extended tube structures are quantitatively equivalent.

Our starting point for the current computational study comes from two types of structural information. Information on the short-range structure—an elementary syn–anti dimer—was previously derived from NMR data[7] and later refined by Miloslavina et al. (unpublished results), for details see Figure S7. The arrangement of these elementary dimers, with a medium-range order extending over several tens of monomers, was revealed by a layer line in the power spectra calculated from cryo-EM images, collected from more homogeneous mutant chlorosomes,[7] which can be used for extracting the overall helicity. The helical nature of the chlorosome is conserved upon variation of the BChl composition, and a Bijvoet pair of spots on the equator reveals well-determined intertube spacings around 2 nm, depending on the species.[2]

In addition, several groups have carried out optical measurements to estimate an angle θ for the effective transition dipole moment relative to the symmetry axis of the tubular aggregates from the dichroic ratio.[8−11] The response to excitation is clearly a collective phenomenon, meaning that this angle does not necessarily translate down to the level of individual molecules.[8] However, the ratio of the oscillator strengths along the two principal directions was shown to directly relate to the angle β, see Figure , for monomer transition dipole moments if one assumes a perfect cylindrical lattice model that satisfies special symmetry conditions,[6,12] which provides an additional geometrical constraint. For bchR mutant chlorosomes, the measured ratio translates into β ≈ 56.1° and β ≈ 52.2° for the two structural components, values that are well in line with the orientation of the molecular dipoles in the original bchQRU model by Ganapathy et al.[7] and a structure for the wide type (WT) by Miloslavina et al. (unpublished results). Nevertheless, because assembly structures can notably deviate from the idealized situation that is at the basis of this explicit relation, we will not regard this information as a strict constraint in this work.

Figure 1

Schematic illustration of the angular relations in a tubular structure. The angle between the cylinder axis (O∥) and the induced molecular transition dipole moment (the so-called monomeric transition) is given by β, its angle with the vector of the Bravais lattice by η, and the chiral angle by δ. For δ between 0° and γ, β = |(δ – 90° + η)|. We note that chiral angle with respect to the parallel symmetry axis, δ∥, is obtained as δ∥ = |δ – 90°|. Moreover, we selected a Bravais reference frame such that the rolling vector varies between and . Hence, the chiral angle defined in Günther et al.[6] differs from ours and is equal to |180° – δ|, with δ our chiral angle.

Using unconstrained all-atom molecular dynamics (MD) simulations, we first evaluate the stability of the pseudosymmetric syn–anti motifs, that is, the principle building block of chlorosomes determined from NMR, and then study the propensity of such syn–anti stacks to form concentric chiral tube structures, thereby departing from the ideal tube model. By performing simulations for an initially planar system, containing either BChl c or BChl d, we resolve the molecular origin of curvature and determine the distinct structural variability and global dynamic heterogeneity of the BChls within a curved assembly. Recent computational work based on molecular mechanics and density functional theory (DFT) has not been able to reveal these basic properties[13] because their phase space exploration capability is too limited.[14] As the time scale τ for structural rearrangement scales as τ ≈ ξ, with ξ the correlation length (∼size) and z > 0 a dynamical exponent that depends on the kinetic model, the starting structure for large-scale MD simulations should be chosen with care to avoid kinetic trapping. We start all simulations from a planar Bravais lattice provided by the most recently proposed triclinic unit cells for BChl d.[7] In studies of electronically excited states, the molecular structure is usually generated by projecting regular two-dimensional (2D) Bravais lattices onto ideal tubes.[6,15−21] As a reference to our MD results, we first refine such idealized models by a geometric analysis that links the position of the layer line for the syn–anti stacks to an overall chirality that satisfies all considered experimental constraints.

The nanosecond MD simulations in this work show that defined curvature develops spontaneously for planar systems and confirms that symmetry breaking is an inherent property of parallel stacking of the pseudosymmetric syn–anti BChl dimer building blocks. According to our simulations, a rotational degree of freedom around the Mg–OH coordination bond enables molecules to resolve frustration, at the expense of giving up O–H···O=C interstack hydrogen bonds, and assume adaptive curvature of sheets or tubes. The signal ratios for the two fractions determined from the splitting of NMR signals matches the ratio of hydrogen-bonded and nonbonded domains in the simulations of a single tube, which provides a first sound structural basis for the presence of two well-defined structural components, randomly distributed over the curved helical structure. We note that, in a helical structure, a majority phase is likely to percolate and that the simulated diffraction response for this structure confirms that these diffractions give rise to a single layer line in reciprocal space, in agreement with the cryo-EM diffraction response. We note that the syn–anti parallel stacking leads to polar sheets and cylinders with extended dipoles along the central axis, rendering the lowest exciton states strongly delocalized, twisting around the tubes.[6,9] Because dipoles and transition dipoles of concentric cylinders overlap in space, correlated dynamics will give rise to fluctuations of the cylinder dipoles and time-dependent detuning of energy levels that may give rise to quantum instabilities. Here, we provide the first detailed structural and dynamical information beyond the current, idealized picture of chlorosomes, which is a valuable starting point for future studies of vibrational-excitonic coupling.

Computational Details

The OPLS-AA force field[22,23] was used for all simulations. All initial structures of pigment molecules (BChl c or d) were optimized using DFT calculations (HF/6-31G* or higher level). The equilibrium bond and angle values in the force field were obtained from optimized DFT structures, giving special care to the Mg atom. To maintain an ∼0.21 nm coordination distance between Mg and O, their Lennard-Jones nonbonded interaction parameters were selected based on Roccatano et al.[24] In addition, four Mg–N bonds were added to keep the Mg atom centered within the ring.

Atomic point-charge parameters were obtained by the RESP method[25] at the HF/6-31G* level using Gaussian03 software.[26] Charges obtained at the B3LYP/6-31G* level were also considered. By comparing them, we focused on their effect on packing by extracting lattice parameters for a small simulated system (5 × 5 dimers). We found that lattice parameters vary only very slightly (within 5%) with the chosen method. For larger simulated systems, we expect equivalent variations.

MD simulations with periodic boundary conditions along all dimensions were performed with the Gromacs 5.1.2 software package.[27,28] A particle-mesh Ewald method[29] was used to treat the long-range electrostatic interactions. We performed simulations in a NVT canonical ensemble, maintaining the simulation temperature by a V-rescale[30] thermostat. Additional details are provided in the Supporting Information simulation method section. We simulated the assembly of 172-farnesyl-(R or S)-[8-ethyl,12-ethyl]BChl c or d pigments, which have been considered extensively in the experimental and simulation literature.[5,7,13,31] The ratio between R and S components in the dimer is 1:1 in agreement with Ganapathy et al.[7]

Our initial simulation structures are planar, obtained by prolongating the proposed syn–anti dimer structure periodically. The lattice is defined in terms of two principle directors and , with length a = || and b = || relating to the experimental repeat distances and their mutual angle γ (see Figure S7C). The whole structure, in terms of the number of periodic repeats of dimers along these two directions, is denoted as (A, B), for example, the structure in FigureS7C is for (5, 5).

Results and Discussion

Refinement of the Ideal Tube Model

Since their development in the 1970s, Fourier–Bessel reconstruction methods have become a major tool for reconstructing the three-dimensional structure of naturally occurring helical assemblies such as actin, myosin filaments, microtubules, amyloid fibrils, and bacterial flagella from 2D electron micrographs. The theory and challenges, associated with indexing and assigning helical symmetry for complicated assemblies similar to proteins, can be found in a number of review papers,[32−35] and explicit geometrical relations between the layer line position and the chiral BChl packing on a tubular lattice in terms of (variable) unit cell parameters are provided in the Supporting Information. For these relations to hold exactly, one only has to assume an ideal tube model for the chlorosome; in the remainder, we will perform unconstrained MD simulations to investigate whether this is justified. We note that this idealized representation renders the geometrical description of chlorosomes and inorganic single-walled nanotubes (SWNT) equivalent, apart from the particular Bravais lattice.

In mathematical terms, our analysis focuses on finding appropriate solutions for the underdetermined system of equations that relate the four principal unknowns, that is, the parameters (a, b, and γ) for the dimeric unit cell and the chiral angle δ, to the layer line position (1/d) determined from TEM. Inspired by the NMR and TEM findings, we consider a δ that is invariant under changes of the tube radius R and conserve the unit cell for chlorosomes that feature different layer line positions. Solutions in previous works were based on selecting chiral angles δ from packing considerations that directly related to unit cell parameters, resulting in the well-known packing models: δ = 90° for the WT with d = a = 1.25 nm and δ = 0° for bchQRU with d = b sin(γ) = 0.83 nm.[7] Since their publication, these models have been debated and further refined.[6,36,37] For example, recent optical experiments for single chlorosomes linked 1/d = 1/1.24 nm–1 for a bchR mutant,[6] equal to the WT value (1/1.25 nm–1) within the considered resolution, to two chiral angles δ = (71 ± 2)° for two different structural components. This is equivalent to δ = (109 ± 2)° in our framework (see the caption of Figure ).

On relating the one visible experimental layer line to geometrical properties of the chiral BChl stacks in the Fourier–Bessel framework, our key challenge boils down to selecting the proper helical family on the tube,[34,35] given that each helical family constitutes one layer line in the reciprocal space. In particular, the distance to the equator for that layer line relates to the reciprocal value of the axial repeat d for that helical family. As discussed in more detail in the Supporting Information, we concentrate on calculating d that corresponds to the principal (0,1) and (1,0) helical families, by a geometrical analysis of the Bravais lattice in real space because they, together with the (1,1) helical family, represent the closest packing or most stable spacings in the tube in the considered range of δ. Moreover, it is known that a small tilt of the tube with respect to the plane of imaging can produce a meridian signal[32] and that replacing idealized point scatters by actual vibrating molecules can give rise to an amplification or weakening of the layer line intensity, which explains why only one reflection of these principal helical families is observed in the noisy power spectrum. A power spectrum calculated directly from a 3-tube aggregate with a distinct layer line at 1/d = 1/1.25 nm–1, Figure S6, clearly illustrates this effect. Adopting the nomenclature of Moody et al.,[38] the axial repeat d is the repeat distance along the tube axis between adjacent helices (for a n-start helical family, d = P/n), and the pitch P is defined as the axial advance after one complete helix turn. Figure shows the theoretically predicted axial repeat distance d versus chiral angle δ for three candidate helical families—d and d, corresponding to the assembly along the and direction, respectively, and darmchair, corresponding to an analogue of the armchair conformation in SWNT—rendered for the equilibrium unit cell parameters of the simulated (30, 30) BChl c system: (a, b, γ) = (1.48,0.98,124.3°) at 300 K. For these unit cell parameters, our geometrical analysis predicts a unique solution δ0.83 = 13° (or, with respect to the tube axis, δ∥ = 77°) for bchQRU and three candidate solutions for the WT. One of these candidates, δ1.25 = 50° (δ∥ = 40°), implies that chirality varies only gradually with a change of the molecular side groups because both BChl c and d stacks wrap along the same direction. The two other candidates for the WT, δ1.25 = 80° (armchair, δ∥ = 10°) and δ1.25 = 112° (δ∥ = 22°) (wrapping along ), imply that a change of the side groups gives rise to a topological transition.

Figure 2

Axial repeat d, d, and darmchair, relating to the H(1,0), H(0,1), and H(1,1) helical families, respectively, as a function of the chiral angle δ ∈ [0, γ]. The solid lines were calculated using geometrical analysis; symbols highlight the same relations for discrete equidistant values of δ. For each helical family, the axial repeat d determines the layer line position in the diffraction pattern as 1/d. We have highlighted special cases, for example, δ = 0,90. For the lattice parameters (a, b, γ), we adopt the values determined from our MD simulation (1.48, 0.98, 124.3°), respectively. Between δ = 0 and δ = γ, the cylinder topology changes, meaning that the lattice rolling direction changes from along to along , with a cross-over chiral angle δc given by a cos(δc) = b cos(γ – δc). The predicted chiral angles based on the layer line position are labeled by white circles: 12.7° for bchQRU type, 49.6° or 112.3° for WT chlorosomes and the corresponding predicted β angles are 46.9°, 10.0°, and 52.7°, respectively, see more details in Figure S9.

The predicted chiral angles δ0.83 = 13° for bchQRU and δ1.25 = 112° for the WT are consistent with the earlier proposed values in Ganapathy et al.[7] if one accounts for the different dimensions of the unit cell. In particular, the -axis in our larger unit cell needs to be tilted to match the experimentally determined pitch length; hence, we find δ1.25 = 112°. Following this reasoning, the unprecedented prediction δ1.25 = 50° can be understood in terms of a similar tilt in the opposite direction.

As mentioned in the introduction, a recent study for individual bchR mutant chlorosomes, composed of slightly modified BChl c, determined an average angle β ≈ 55° (with respect to the tube axis) for the induced monomer dipoles, from the ratio of parallel and perpendicular oscillator strengths. Yet, the angle θ estimated from the dichroic ratio by earlier studies is often smaller (θ = 15–25°) or distributed over a broader range, which has been interpreted in terms of structural variation between species[11,39] or in terms of the preparation method.[9,36] Straightforward analysis of the average angle β between the orientation of the BChl head groups and the symmetry axis in our setup, see Figure S9, shows that δ1.25 = 112° agrees well with β ≈ 55°. Nevertheless, MD simulations show that prepacked nested three-tube BChl c assemblies are stable for both tested values of δ1.25, that is, δ1.25 = 50° and 112°, and layer line positions determined by FFT are in excellent agreement with the predicted values.

We may benefit from the fact that we optimized local structure by MD to extract specific unit cell parameters directly from the simulated (30, 30) BChl d system (see Figure S8 for details). We considered the most cylindrical structure (1 ns at 50 K) to allow for easy extraction of unit cell parameters and δ, giving rise to (a, b, γ) = (1.46, 0.93, 122.5°) and δ = 19.1°. Inserting these unit cell parameters in the general expression for d provides the alternative unique solution δ = 19.8° for d = 0.83 nm. Comparing it to the extracted chiral angle δ = 19.1°, we find that they match very well; therefore, we will consider δ0.83 = 19.8° in the remainder. Thus, although our simulated BChl d aggregate is too partial to provide a well-defined layer line by a direct procedure, see discussion in the Supporting Information, this equivalence between simulated and theoretical values leads us to conclude that the longer-range structural properties of the simulated BChl d aggregate match very well with the ones for the bchQRU mutant containing an excess of [E, M]BChl d.

As a comment to earlier work,[6] we find that the experimentally observed sequence of tube radii is fully consistent with two fixed chiral angles, δ0.83 and δ1.25, for nested tubes in the two systems. We considered δ1.25 = 112° in this analysis, but note that it is illustrative for the other two values of δ1.25. Restricting the maximum variation of δ to a tiny |Δδ| < 0.2°, lattice periodicity plays a role in determining which radii are commensurate, but we find that commensurate R very consistently matches the experimental values, see Figure S4 in the Supporting Information. This analysis also shows that a maximum variation of only 1° is sufficient to produce any tube radius R, given that the discrete increment δR that stems from lattice restriction is on the order of the experimental resolution. Our simulations thus show that the tube radius is essentially a scale-free parameter, that is, not restricted by any inherent length-scale in the system. This interesting observation agrees well with the nucleation and growth formation kinetics of chlorosomes[31,40] because it allows for flexibility to adopt any radius to tightly wrap new layers around an existing (and possibly defected) nucleus.

MD Simulation Results

Unconstrained MD simulation allows us to gain fundamental insight in the deviations from the ideal case and the origin of curvature, as well as relate local structure to global properties of the assembly. Our analysis of MD simulation results focusses on the following complementary aspects: (I) large-scale structure and evolution, as well as properties (radius R, pitch P, and chiral/pitch angle δ) of the spontaneously formed tubular structures, including the overall conformation of the farnesyl-tails. (II) Local structural details by an in situ analysis of the curved aggregates, concentrating on stabilizing factors and interactions, and on the relative rotation between syn–anti and anti–syn pairs, which provides quantitative insight into the driving factors for curvature formation. (III) Since local dynamics changes the detailed network of interactions between BChls, in particular, the hydrogen bonding, we also evaluated the dynamics of this rotational degree of freedom and its role in the formation of hydrogen bonded domains. We focused on the angular sampling statistics for an isolated dimer and consequently evaluated the role of these twisting motion in our simulated structures.

Large-Scale Structure: Intrinsic and Flexible Curvature

Our MD simulations provide computational evidence that the proposed local syn–anti stacking model[7] gives rise to a stable aggregate, with initial forces that are well within bounds and a tendency to spontaneously curve on the scale of the aggregate. Figure S12 shows MD simulation snapshots for an initially planar sheet, built from (30, 30) BChl c syn–anti dimer stacks, that wraps quickly as a whole to form a helical tube structure well within 1 ns (full trajectory available in Supporting Information Movie S1). The obtained curved structure remains stable when further simulated in hexane solvents. To determine the effective tube radius at the end of the simulation, we used a least-square algorithm to fit a cylinder to the coordinates of the Mg atoms, see Figure S12a. Helical features of the BChl packing along the tube can be observed both visually and from the characteristic peaks in the calculated 2D Fourier spectrum for the Mg atoms, see Figure S12b. The fast kinetics and accompanying jump in the potential energy (Figure S12c) suggest that the curved structure is thermodynamically much more stable than the planar structure, and that the driving force for this transition is rather large. This argues against a curvature-invariant system, with curvature purely induced by thermal fluctuations, which is usually a very slow process.

The packing parameters obtained by simulation, (a, b, γ) = (1.48, 0.98, 124.3°) for BChl c and (1.46, 0.93, 122.5°) for BChl d, see Table S1, are based on the analysis of the local packing of neighboring dimer units assuming a planar arrangement, which is a valid assumption for large tube radii. Overall, we find that the packing into a curved tubular structure is quite regular, with small static heterogeneity, which is in agreement with the visual observation of regular packing of Mg atoms (Figure S12b). We have additionally applied a procedure introduced by Connolly[41] on part of the simulated structure, using the standard 1.4 Å Connolly radius, to determine its effective volume and consequently the BChl c density as 1.348 g/cm3, which compares well with the density of 1.31 g/cm3 determined by X-ray crystallography for ethyl chlorophyllide a dihydrate.[42]

The tube structure obtained by 50 K NVT simulation remains stable after increasing the temperature to 300 K. The packing parameters show almost no dependence on the simulation temperature, that is, the relative change in the packing parameters is within 2% upon equilibration at 300 K, which is negligible. Yet, globally, the radius of the fitted cylinder can be seen to change significantly from 9.335 to 6.840 nm (Figure b). This variation of the tube radius, for a structure with very similar local packing features, reveals the sensitivity of our assembly structure to temperature, and the flexibility to adjust the assembly properties on a large scale while conserving the packing on a local scale. We further consider this relation by performing a geometrical analysis.

Figure 3

Simulation results obtained for different systems: (a) (30, 30) BChl c system at 50 K, (b) (30, 30) BChl c system at 300 K, which started from the final structure obtained at 50 K, and (c) (40, 30) m-BChl c system at 300 K, i.e., in the absence of the farnesyl tail. Cylinders are fitted to the coordinates of the Mg atoms using the same procedure as in Figure S12. All sizes are scaled equivalently, enabling a proper visual analysis of relative sizes. From left to right, the radius of the fitted cylinder is 9.3, 6.8, and 6.1 nm, respectively. For clarity, side groups of the head part and tails are not shown.

Large-Scale Structure: Helicity

The simulated assemblies for BChl c are defected but clearly tubular and, on a local scale, pigments can be seen to adopt a helical arrangement; therefore, it makes sense to perform a quantitative analysis in terms of large-scale structural properties: the tube radius R and pitch P and the earlier discussed chiral angle δ of the helix, also known as the pitch angle. For structures obtained after 1 ns simulation, simulated at different temperatures, see Figure b, we find that both R and P decrease monotonically with increasing temperature. However, the third parameter, δ, adopts a constant value of δ ≈ 31°. The only exception is T = 50 K, where the chiral angle is significantly lower (Figure c). However, as our simulations are performed sequentially, using the 1 ns result obtained at 50 K as input for a simulation at higher T after (local) equilibration, the deviating chiral angle for 50 K can be explained in terms of the reduced simulation time, if we disregard the alternative explanation of a genuine structural transition in the range of 50–100 K.

Figure 4

(a) Illustration of the principle helical family H(1,0) and its relation to the radius, pitch, and chiral angle. (b,c) Large-scale properties (cylinder radius, pitch, and chiral angle) for temperatures in the range 50–300 K, as determined for the H(1,0) helical family for the (30, 30) BChl c system after 1 ns simulation: (b) Fitted cylinder radii R and pitch P at different simulation temperatures and (c) fitted chiral angles δ at different simulation temperatures. With increased temperature, both the radius and pitch decrease, whereas the chiral angle first increases and then assumes a constant value, suggesting a structural transition between 50 and 100 K. These results show that the chiral angle is the only conserved structural property, with the pitch adapting to a particular radius.

All simulated curved aggregates are only partial tubes after 1 ns simulation, meaning that the periodicity imposed by the Bravais lattice, which may affect helical properties, does not play a role at that stage. In the absence of the farnesyl tail, see Figure , the initially flat m-BChl c sheet curves and closes upon itself within 1 ns to form a closed (complete) tube. The extracted chiral angle δ for this system is somewhat larger (∼40°) than that for the partial tubes, but we lack insight in the particular role of the tail on this property. An axial repeat d = 1.0 nm, close to b = 0.98 nm, was extracted (see Figure S5) from the layer line in the simulated electron microscopy images for the full tube. This value agrees with the geometrical prediction for δ = 40° in Figure and the real-space observation that the tube symmetry axis is almost parallel to .

All partial tubular structures appear stable on the 1 ns time scale of the simulation. Continuing the simulation for the highest considered temperature, T = 300 K, shows that the partial tube does close up after roughly 2 ns. Closing up stabilizes the structure: the tube diameter is constant and the rugged connection zone does not notably change, either via (local) reorganization of dis-/reconnection, in the following 8 ns simulation. Apparently, although a flat structure is unstable, the thermodynamic driving force for structural reorganization away from a nonoptimal δ and/or radius is tiny and/or the rearrangement kinetics is too slow to be captured by MD. On the basis of these arguments, we conclude that the processing history and/or the environment (e.g., solvent or lipids) is vital for obtaining the appropriate chirality by MD, which agrees well with experimental information on this system.[43] We may, however, assume that enhanced kinetics will give rise to a chiral structure featuring the predicted δ1.25 = 50° rather than the other two options, δ1.25 = 80° and δ1.25 = 112°, which require a topological transition. We also note that partial tube structures that are observed in experiments can be interpreted as a lack of material needed for fully closing up. Overall, these findings are consistent with a templated formation of concentric tubes, in which newly forming layers can easily adapt their curvature to an existing nucleus while reproducing the local packing characteristics of that nucleus. Helicity plays a key role in this proposed mechanism, as a helical pigment arrangement is the only option for accommodating a continuously varying tube radius.

Large-Scale Structure: Symmetry Breaking

In principle, there are two out-of-plane directions in which the planar symmetry can be broken. Each of them gives rise to a curved structure with different chirality and with different (syn or anti) tail domains at the concave or convex side of the curved structure. We may distinguish these options in terms of ± × , that is, the cross-product of two principle directors, see Figure S7c. In particular, tails distribute to either side of the tube, depending on the head type (syn or anti) to which they are attached, owing to the mutually parallel orientation of the pigment heads in the elementary syn–anti dimer,[7] showing that local packing plays a key role in the distribution of tails. This tail distribution is in strong contrast to antiparallel head–head dimer stacking which would lead to both tails naturally ending up at one side of the stack; see Pandit et al.’s work[44] for a classification of molecular stacking options for related zinc-chlorins.

Despite limited statistics, only one curving direction (− × ) is consistently found in all simulations. Also, without the (farnesyl) tails, that is, when considering m-BChl c instead, the head stacks curve in the same direction. This preferential direction is equivalent to syn tails always ending up at the concave side for all curved structures (see Figure ). The appearance of a preferential helical chirality is a feature of chiral asymmetry.

Figure 5

Details of the syn–anti packing for the curved structure obtained after 1 ns MD simulation of (30, 30) BChl c system at 50 K. Syn and anti BChl c molecules are shown in red and blue, respectively. (a) Curved structure (repeated from Figure S12a) and fitted cylinder. For clarity, side groups of head part and tails are not shown; (b) syn–anti packing of heads, showing high order of head parts even in curved domain; (c) projection of the pigment assembly along the cylinder axis with and without tails, displaying the spatial extent of the tails in the curved structure. The concave (convex) side of the curved structure is covered by syn (anti) tails, respectively; and (d) details of the tail structure, illustrating the groovelike nature of the tail phase at the concave side of the curved structure, leaving enough free space for other molecules to occupy, as well as the liquid-crystal-like order of the tail phase for this low temperature.

The actual chirality of natural/mutant chlorosomes could be extracted by circular dichroism (CD), which would enable us to validate our finding of chiral asymmetry. In practice, however, CD spectra for such psi-type assemblies show undesired variability.[2,11,45] Theoretically, the magnitude and sign of these spectra were not only shown sensitive to various particular structural properties[46] but also the same three-band (−, +, −) CD signal was calculated for quite different tubular dipole arrangements by related methods based on point-dipole approximations.[8,12,37] Nevertheless, for chlorosomes that most closely match the studied systems, that is, the bchQRU containing primarily BChl d(47) and an artificial one composed of only BChl cs,[11] the same type I spectrum[11] was measured, which may be attributed to our type of dipole arrangements based on the available information.[12]

Distributed between the tubular stacks, tails determine the intertube structure and spacing. Ideal tube models, that is, single tubes obtained from a lattice of dimers by projection,[6,15−21] usually consider tails to be fully extended orthogonal to the tubular domain, as the experimental data do not provide sufficient information on tail conformations in the assembly. We find that tails adhere to head groups belonging to adjacent dimers, see Figure d, which explains why simulations for dimers fail to capture these tail configurations. Consequently, the orthogonal projection of the end-to-end distance (d⊥) is quite small (|d⊥|/|d| ≈ 0.2, see Figure S10), when compared to the fully extended chains (|d⊥|/|d| ≈ 1.0) in current projection models.

To analyze whether our tail configurations directly relate to the interlayer distance dexp = 2.0–2.1 nm determined from Bijvoet pairs for multitubular WT and mutant chlorosomes,[7] we determined the thickness of the simulated tail layer. For our aggregates, which comprise a single tube in vacuo, we find that the interlayer distance matches the experimental value 2.1 nm quite well (see Figure S11). In contrast, fully stretched tails in ideal models would translate into a substantially increased interlayer distance. This finding gives us confidence to conclude that the tail configurations in our MD simulations are also realistic for a multitube system, suggesting that tails from different tubes do not experience notable interpenetration.

From the simulation snapshot shown in Figure d, tails can be seen to bundle together on a larger scale, leaving sufficient space to accommodate small molecules similar to the carotenoids that were extracted from natural chlorosomes.[48] Although less ordered than heads, tails adopt a liquid crystallike order at low temperatures when considered locally (Figure b,d).

Although tails fill the interstitial space between tube layers and determine their spatial separation, they are not a key factor in curvature formation. The latter can be concluded from simulation of a (40,30) m-BChl c sheet at 300 K (Figure c), where large scale curvature forms spontaneously in the absence of farnesyl tails. Thus, we conclude that curvature formation originates from head–head packing.

In Situ Analysis: Stabilizing Factors

Next, we analyze the interaction network between neighboring BChl c by considering the nearest distance distribution of interacting groups in the head domain and their overlap ratio. Two types of pair interactions stand out: (i) the interactions between magnesium and oxygen atoms (Mg···O–H in Figure a) coordinate neighboring pigments along the packing direction . (ii) Hydrogen bonds formed between the hydroxyl group and the carbonyl oxygen, see O–H···O= in Figure a, stabilize the molecular assembly along the direction. Because each pigment molecule contains a donor (hydroxyl group oxygen) and an acceptor (carbonyl group oxygen), which are both capable of forming hydrogen bonds, three states are needed to describe this bonding: “2-hb”, that is, both donor and acceptor form hydrogen bonds, “1-hb”, that is, either donor or acceptor forms a hydrogen bond, and “0-hb”, that is, no hydrogen bond is formed.

Figure 6

Analysis of key interactions in the assembly, as determined from the curved structure obtained by MD simulation of (30, 30) BChl c system at 50 K after 1 ns. (a) Schematic representation of the interacting groups that promote pigment assembly along the packing direction and ; (b) distribution of nearest distances between Mg and O atoms, important for coordination between Mg and OH groups; (c) distribution of nearest distances between OH and O=, important for hydrogen bonds formed between neighboring carbonyl oxygen and hydroxyl group; (d) distribution of nearest distances between stacking layers, important for π–π interactions. The mean value is 0.356 nm, which is a characteristic value for π–π stacking; and (e) head–head stacking overlap ratio along and directions. The overlap ratio calculation is based on the van der Waals geometry of carbon atoms in the head part, ignoring substitutes.

We find that the distribution of Mg···O–H distances is narrow with an average around 0.21 nm (Figure b). In contrast to the Mg and O–H groups, which always coordinate to each other, we find that not all neighboring O= and O–H groups always form hydrogen bonds. We quantify this property by considering the ratio between the actual number of hydrogen bonds between donors and acceptors and the total number of donors and acceptors. For the (30, 30) BChl c system at 50 K, this ratio is approximately 0.7, meaning that 70% of all possible hydrogen bonds is actually formed, see details in Figure S13. We note that, for a system with a similar packing, modification of groups that are capable of hydrogen bonding did not destabilize the assembly.[49] This confirms our finding that the formation of hydrogen bonds along the direction is not vital for the stability of the assembly structure. To conclude this analysis, the ratio extracted from the simulated BChl d structures at 50 K is around 0.9, meaning that both simulated values are consistent with the signal-splitting ratios in the NMR experiments of chlorosome containing an excess of BChl c (7:3) and BChl d (9:1).[5,7]

We attribute the overall stability of the assembly along both and directions to the unique π–π stacking between the head parts along both directions. Earlier static calculations for Zn-chlorins[49] confirm this conclusion. The packing distance between head parts for our system is found to be 0.356 nm (Figure d), as determined from the average line distance, which is a characteristic value for π–π stacking. The really unique feature in this structure is that, along the direction, heads only partially overlap owing to the Mg···O–H coordination, which is compensated by a partial overlap along the direction: the head–head packing overlap ratios along the and directions were determined as 0.39 and 0.2, respectively, see Figure e. In this way, π–π stacking plays a role along both directions, strengthening the interaction framework of the assembled structure.

In Situ Analysis: Curvature-Inducing Factors

Dimers that satisfy the key stabilizing factors of the larger assembly, that is, coordination and π–π overlap, have only one degree of freedom left: to adjust their rotation relative to their neighbors. In particular, one can anticipate that curvature should originate from consistent adjustment of relative rotation angles throughout the whole structure. For this reason, we extracted the relative rotational angle α between adjacent pigments, which is defined as the angle between Mg···O= and Mg···Mg vectors after being projected on one molecule’s rigid ring surface (Figure a), with the aim of providing quantitative molecular insight in the origin of large-scale curvature formation.

Figure 7

Orientational details of head–head stacking, extracted from the curved structure obtained in the MD simulation of (30, 30) BChl c system at 50 K, after 1 ns. (a) Typical intradimer and interdimer stacking configurations associated with a positive and negative relative rotational angles, respectively. Although the syn–anti dimer is the basic structural element; on the level of individual pigment, there is an alternating packing of syn–anti (intradimer) and anti–syn (interdimer) configurations. The relative rotational angle α is defined by the angle between Mg–O= and Mg–Mg vectors (labeled as red and blue arrows) after being projected on the top molecule’s rigid ring surface. The sign is determined by the direction of the cross-product of the two vectors; (b) individual values for the rotational angles along the assembly structure; and (c) distribution of intradimer and interdimer rotation angles, illustrating the asymmetry in packing.

It should be understood that, although a syn–anti dimer is the structural element, packing in a larger structure gives rise to alternating syn–anti (intradimer, αintra) and anti–syn (interdimer, αinter) pairs, see Figure . The sign of α is determined by the cross-product of the two key vectors; for typical pairs, αintra > 0 and αinter < 0. Analyzing all pairs in the curved structure for T = 50 K (1 ns), we find that the signature of alternating positive and negative α is well-maintained throughout the large assembly (Figure b). Analyzing their distributions over the curved structure, see Figure c, shows that a broad range of α values is sampled. This rotational heterogeneity may contribute to the broadening of the optical spectra. Closer examination suggests a slight asymmetry, that is, ⟨αintra⟩ + ⟨αinter⟩ = −2.78°. Improving the statistics, however, by adding snapshots of a tubular structure along the simulation trajectory, does not give a unimodal distribution but clarifies that the sampling of this rotational degree of freedom is broad (see Figure S14) with a mild preference for specific α values.

To allow for more firm conclusions about symmetry breaking, we performed additional analysis. In particular, we introduced a simplified representation with handedness and consider Δβ = βintra – βinter, with the angles βintra and βinter determined as in the inset of Figure . Using geometrical arguments, it is easy to show that one obtains a tubular structure with helical ordering for , where the bar represents the average value. Analysis of simulated curved structures along the simulation trajectory at 50 K indeed suggests that such a nonzero exists, as can be seen from the peak position of the Δβ distribution (Figure ). Increasing the temperature (to T = 300 K) transforms this bimodal distribution into a unimodal Gaussian-like distribution centered around . The unimodal nature of this distribution is clearly the result of increased thermal fluctuations (Figure S15), which helps the system to overcome energy barriers along this rotational coordinate. This finding is consistent with the large-scale analysis of the helical tube structures (see earlier section).

Figure 8

Distribution of difference (βintra – βinter) between neighboring intradimer and interdimer rotation angles in the representation that incorporates handiness, as extracted from the curved structures along a MD trajectory for a (30, 30) BChl c system at 50 K, between 500 ps to 1 ns. The inset shows the intradimer and interdimer configurations and the way the corresponding angles βintra and βinter are determined. When the intradimer (anti–syn) and interdimer (syn–anti) configurations are (anti-)symmetric, the packing-angle difference βintra – βinter is 0°. Nevertheless, a major distribution peak around 12.6° > 0° appears, suggesting consistent asymmetry between the chiral intradimer and interdimer configurations in the curved structure. Such asymmetry is the origin of the transition from the translational symmetry in the initial flat structure to helix symmetry in the simulated curved structure. Moreover, reconstructing structures containing alternating syn and anti representation with a fixed positive βintra – βinter value gives rise to an anticlockwise growth and a curved structure. Such a curved structure is consistent with the curved structure obtained in simulation when syn tails are located in the concave side of the curved structure.

Our detailed analysis illustrates how the pseudosymmetry of the elementary structure develops into a well-defined asymmetry, via symmetry breaking in the head–head packing, which is the basis for large-scale curvature formation at a molecular level. The global flexibility to adopt varying curvature should also be seen to originate from the local head–head packing. Because the distribution of Δβ is rather broad, the assembly is rather flexible in incorporating various curvatures on a larger, global scale.

In Situ Analysis: Hydrogen Bonding

We consider the curved structure at 50 K to analyze the formation of hydrogen bonds using the three-state description introduced earlier. We find that different states are organized into striped domains on the projected 2D structure, see arrows in Figure , meaning that they form helices on the tube, as discussed in the theory section in the Supporting Information. Interestingly, almost all “0-hb” states (no hydrogen bonding) are next to “2-hb” states (fully hydrogen bonded).

Figure 9

Hydrogen bonding states for individual molecules along the curved structure obtained in MD simulation of (30, 30) BChl c system at 50 K, after 1 ns. Depending on the number of hydrogen bond formed on a molecule, three states are defined: “2-hb”(dark blue), “1-hb”(light blue), and “0-hb” (white).

With the appearance of “0-hb” and “1-hb” states in Figure , it is clear that not all residues capable of hydrogen bonding are “saturated”. As mentioned previously, around 70% of these residues do form hydrogen bonds in this curved BChl c structure at 50 K. Correlating the two local properties, rotation and hydrogen bonding, via the absolute value of the relative rotation angle distribution |α|no-hbonding for molecules that are not-hydrogen bonded (the absolute value of the rotation angle is sufficient because we do not need to distinguish between inter- and intradimer pairs) shows that hydrogen bonds are broken when the relative rotations exceed some threshold value, with |α|no-hbonding distribution peaking around 35°. In particular, hydrogen bonds do not interfere with the interaction between neighboring Mg and carbonyl group oxygen atoms (the Mg···O= interaction is strongest for |α| = 0°, see Figure S16 for distribution and analysis details.)

Dynamics: Rotational Sampling Distribution

Although the total number of hydrogen bonds can be seen to converge to a constant value with time, see Figure S13, their distribution over the curved structure varies considerably with time. Hydrogen bonds between molecules in neighboring stacks are formed and broken in a concerted manner to conserve their overall number, but no diffusion of striped “0-hb”, “1-hb”, or “2-hb” domains as a whole was observed, see Supporting Information Movie S2. In combination with the narrow distribution of donor–acceptor spacings, see Figure , these findings clarify the correlation between the relative rotation angle α and hydrogen bonding between BChl stacks and signify the role of individual BChl c rotations in the curved structure.

We further zoomed-in by performing spectral analysis of time traces of the relative rotational angle α for a number of well-chosen dimers, making sure to select combinations of monomers that belong to different states, that is, “0-hb”, “1-hb”, or “2-hb”, at the end of the trajectory. Analysis of individual traces shows that hydrogen bonding does affect this rotational dynamics, by modulating the amplitude of the variations in α. Yet, we find that all calculated spectra in the frequency domain consistently show a band around 125–180 cm–1, see the Supporting Information for details. We tentatively assign this band to the 145 cm–1 low-frequency mode observed by ultrafast spectroscopy.[3,4]

Finally, we analyzed the sampling distribution of the relative rotational angle α, see Figure , for a dimer to obtain a better understanding of the rotational dynamics in the larger aggregate. This reduced dimeric system enabled us to perform a much longer (50 ns) MD simulation, which considerably improves the sampling statistics. Our choice for the m-BChl c dimer means that we neglect the influence of tails and of the matrix of surrounding dimers.

Figure 10

Distribution of the relative rotation angle for m-BChl c dimer, obtained by simulating a syn–anti dimer at 200 K. The whole 50 ns MD simulation trajectory was considered. Three characteristic configurations are selected for potential energy comparison between fully atomistic and quantum modeling, see Table , and are denoted by arrows; the fourth, a transition state, is only considered in the context of MD. The simulation temperature was chosen to provide sufficient phase space sampling while maintaining a chlorosome-like syn–anti packing.

Table 1

Potential Energy Comparison of the Three Stable Characteristic Configurationsa

state	I	II	III	IV
angle α (deg)	–0.2	21.5	57.5	98.1
MD (kJ/mol)b	0	4.0	15.4	9.0
AIMD (kJ/mol)c	0	15.7		36.7

For MD, the potential energy of the transition state is added to give an idea of the energy barrier. For AIMD, the average energy of the reference state is reported, and the reported energies for other configurations are differences.

Reference value 2563.9 (kJ/mol).

Reference value −782.449 (hartree).

The simulated rotational angle distribution for this dimer reveals that several configurations are stable: a higher probability relates to more stable configurations or a lower potential energy via the Boltzmann factor. We identify two dominant rotation angles: αI ≈ 0° and αII ≈ 20°. Indeed, the Mg and O= groups are in closest contact for αI, which is associated with the highest probability density. For αII, the probability is only slightly reduced and the intermediate depression is not very pronounced, suggesting that the energetic barrier between states αI and αII is rather modest. Beyond these two peaks, that is, α > αII, the clear drop of probability density hints at the presence of a more pronounced energetic barrier, which is consistent with the probability distribution data extracted for the large assembly in Figure c.

Four characteristic configurations (labeled by the associated α) were selected for direct potential energy comparison (Figure ): three values for stable states (αI,II,IV) along this angular coordinate and an intermediate or transition state (αIII). We find that αI indeed has the lowest potential energy, followed by αII with an increase in the order of 5 kJ/mol. The change of the potential energy for αIV is still modest (see Table ), whereas the transition state αIII is indeed associated with the highest increase. The reduced head–head overlap that accompanies an increasing α clearly gives rise to destabilization of the structure or an increased potential energy. This is consistent with the larger assembly structure, where the probabilities of these larger α’s are indeed considerably reduced, see Figure . To test the consistency of our MD force fields, we complemented this analysis by a potential energy calculation by ab initio MD (AIMD), using the characteristic configurations from MD as input for AIMD equilibration. AIMD confirms the stability of the three stable MD states (see Supporting Information for details of the procedure). Moreover, the potential energy trends are very similar, see Table , albeit that the potential energy surface for MD is clearly flattened compared to AIMD.

Final Discussion

The heterogeneity in chlorosomes appears because of frustration, as a result of the competition between the stabilizing network of interactions and curvature-driven rotation away from the preferential conformation in a syn–anti dimer motif. The aggregation of BChl in sheets and tubular assemblies are apparently not compatible with the formation of an array of H bonds that are commensurate with the Bravais lattice formed by the BChl. The resulting H-bonding defects form stripes that are randomly distributed over the system but appear pinned to the underlying lattice. In particular, because the rotational mode is soft, H bonds can be switched on and off, and we may think of parallels to the folding funnel in proteins. The random distribution of pinned defects over the pseudosymmetric lattice shows that every chlorosome is unique all the way down to the supramolecular level.

Our computational results provide an alternative explanation for the two sets of superradiant states that were tentatively attributed to two types of tubes in a single chlorosome:[6] the full width at half maximum of the frequency spectrum associated with our rotational mode is 30–40 cm–1, that is, well within the residual broadening of ∼100 cm–1, inferred from fluorescence excitation experiments, that should cover both static and dynamic heterogeneity. In the ground state of the imperfect lattice, several modes are enhanced and contribute to a disordered energy landscape. Because the energy levels of the BChl will depend on its position along the rotational coordinate and the H bonding, it is tempting to conclude that such modes can give rise to nonadiabatic coupling elements for quasi coherent energy transfer by resonant coupling upon crossing of exciton levels. The relatively strong vibrational band centered around 160 cm–1 agrees well with the time scale of 100 fs for rapid exciton transfer in chlorosomes determined by coherent optical spectroscopy.[3]

Conclusions

Detailed information on the structure and heterogeneity in natural chlorosomes provides necessary insight toward the understanding of the unique excitonic properties of these assemblies and enables the design of synthetic mimics with similar efficiency. In this study, we have combined a formal geometrical analysis and unconstrained all-atom MD to obtain unique insight in the direct relation between molecular detail and structure, curvature formation, symmetry breaking, and stabilizing factors, as well as static and dynamic heterogeneity within a curved structure. An intriguing finding is that out-of-plane rotation and hydrogen bonding compete, which identifies a new type of heterogeneity that is consistent with NMR observations.

From a structural perspective, the overarching question is how and why nature selects particular chiral angles for structuring within the considered chlorosomes. Here, the “how” is clarified: by building in asymmetry at the molecular level, which translates into specific chirality at the macroscale. We find that our MD simulations in vacuo properly capture this mechanism for the BChl d system, where the δ = 19.1° extracted from the simulated aggregate agrees very well with the predicted δ0.83 = 19.8° of bchQRU. For BChl c, our theoretical analysis predicts three chiral angles δ1.25 that match the experimental NMR and cryo-EM information. The tube topology for the chiral angles extracted from simulation, δ = 31° with farnesyl tail and δ = 40° without tails, that is, curved along the direction of strongest interactions, agrees with one of these predictions, δ1.25 = 50°. Yet, optical measurements for single chlorosomes[6] were combined with theory to predict an orientation of the monomeric transition dipoles that only agrees with δ1.25 = 112°, that is, curved along the direction, albeit that this prediction is uncertified for structures that do not agree with the special symmetry conditions of the theory.[12] An independent way to ascertain the actual helicity is to numerically calculate the optical spectra for our candidates based on a Frenkel Hamiltonian,[8] which will additionally allow us to evaluate details, such as signals originating from two distinct structural components that are lost in the averaging procedure. Nanosecond MD simulations of preassembled structures (single and multiple tubes) for the alternative δ1.25 confirm their stability.

Overall, clearly many factors secure the molecular organization in natural chlorosomes, such as the formation kinetics, matrix templating, and/or nesting, which agrees well with the experimental finding that spectroscopic properties of reconstituted structures crucially depend on processing conditions. The detailed insight provided in this study, into the molecular origin of large-scale curvature and structure, local packing and dynamics within the confines of the curved aggregates, has a general value and validity and will serve as input and inspiration for more detailed computational investigations of excitonic-vibrational coupling.