Optimal number of pigments in photosynthetic complexes

We study excitation energy transfer in a simple model of photosynthetic complex. The model, described by Lindblad equation, consists of pigments interacting via dipole-dipole interaction. Overlapping of pigments induces an on-site energy disorder, providing a mechanism for blocking the excitation transfer. Based on the average efficiency as well as robustness of random configurations of pigments, we calculate the optimal number of pigments that should be enclosed in a pigment-protein complex of a given size. The results suggest that a large fraction of pigment configurations are efficient as well as robust if the number of pigments is properly chosen. We compare optimal results of the model to the structure of pigment-protein complexes as found in nature, finding good agreement.


Introduction
Photosynthesis is the main natural process for harvesting Sun's energy on Earth, providing a food source for a great variety of organisms ranging from highly evolved photosynthetic systems in higher plants to simpler bacteria [1]. In spite of such diversity, basic underlying principles are shared among majority of light harvesting organisms. The initial stage of photosynthesis usually involves multiple pigment protein complexes (PPC) that consist of a number of pigment molecules (i.e., chromophores) held in place by a protein cage. PPCs are employed to absorb incoming photons as well as to transport the resulting excitation to the reaction center, where the excitation is used to initiate chemical reactions. The absorption of light and in particular channeling of the absorbed energy to the reaction center is known to achieve high efficiency [1].
One of the most studied PPCs is the Fenna-Matthews-Olson (FMO) complex that is found in the photosynthetic apparatus of green sulfur bacteria. It is the first PPC for which the atomic structure has been determined [2]. FMO is composed of a large protein envelope that encloses a tightly packed group of 7 pigment molecules ‡ called bacteriochlorophylls (BChl). In FMO the main role of BChl pigments is actually not to absorb light but instead to transport electronic excitations from the input pigment, which is close to the antenna complex, towards a "sink" pigment channeling excitation to the reaction center. Recently discovered long-lasting quantum coherence in FMO complex [4] gave an additional boost to studies of excitation transport in PPCs. Previously it was namely believed that the transport in PPCs is predominantly of a classical nature. Many aspects of time dependence have been studied [5,6,7], including the functional role of quantum coherences [8,9,10,7], as well as the possibility of transport enhancement via environmental interaction [11,12,13,14].
Also important is the question of structural characteristics of PPCs that enable efficient excitation transfer, for instance, why has a given complex precisely the shape found in nature. Positions and orientations of pigments, known with high precision from crystallographic measurements, on casual inspection show no clear ordering or organization that would enable an easy classification of efficient configurations. It is also not clear, how special efficient configurations are -whether efficient configurations are a result of a long-lasting process of evolutionary improvement, or, can efficient configurations be readily achieved probing few random configurations of pigments. Prior to the availability of crystallographic structural data, the average minimal distance between pigment molecules was estimated based on the comparison between florescence yield of in vivo and in vitro chlorophyll solutions [15]. Recently, the efficiency of random configurations within simplified models of PPC was inspected [16,17,18,19], suggesting that the efficient configurations are relatively probable. This complies with the results obtained for the Photosystem I from plants and cyanobacteria [20,21], where random orientations of pigment molecules were probed, and high efficiencies of excitation energy transfer (EET) were obtained irrespectively of the pigment orientation. Also, the PPC configuration was shown to be robust to the removal of a pigment from the complex [20].
Fundamental, yet still unanswered question that we address in present work is, why has a particular photosynthetic complex exactly the specified number of pigments. In other words, what is the optimal number of pigments for a given size of the complex §, or, equivalently, what is the optimal size of the complex for a given number of pigments. For instance, why do we find exactly 7 pigments in the FMO complex and not more or less. Using a simple model whose parameters are taken from experiments, that is without any fitting parameters, we calculate the optimal number of photosynthetic pigment molecules for different complex sizes and compare these theoretical predictions with the actual number of pigments found in naturally occurring complexes. We find a very good agreement for a variety of PPCs in different organisms. To judge the optimality we use two criteria: (i) the average efficiency of the excitation transfer, where the averaging is performed over random configurations of pigment molecules, and (ii) the robustness of the efficiency to small variations of pigment's locations. A rationale behind these two choices is that "good" PPCs should have high efficiency but at the same time also be robust. A specially "tuned" configuration, that has a very high efficiency which though is very fragile, will obviously not work in a natural environment with its changing § Note that the size of the PPC might be fixed by external factors like for instance the membrane thickness.
conditions. Additionally, from an evolutionary perspective, the efficiency should be stable with respect to different foldings of the protein cage. It is advantageous to have a PPC with such a number of pigments that will results in high average efficiency, i.e., in many close-to-optimal pigment configurations. Thereby, a small change in the environment, be it of a chemical origin or for instance a genetic mutation, will still result in a functional PPC. Robustness of quantum coherence to structural changes in the PPCs has been also found experimentally [22]. The model we use to describe the excitation transfer across the PPC consists of a Lindblad master equation describing a dipole-coupled pigments with an on-site excitonic energies being determined by the distances between disc-like pigments. If two discs come too close, i.e., if they overlap, this effectively rescales their energies, introducing a disorder. We should say that the optimal number of pigments that we predict is quite insensitive to details of the underlying model. Agreement between predictions of our model and naturally-occurring PPCs shows that Nature has optimized PPCs by using just the right number of pigments so that the resulting PPCs are highly efficient and robust at the same time.

The model
To calculate the efficiency of the excitation transfer in the PPC, we need equations of motion describing dynamics of excitation on multiple chromophores, coupled to the environment. It is the ratio of chromophore-chromophore interaction strength to the chromophore-environment coupling that determines the applicability of various models. When interaction between chromophores is small compared to the environmental coupling the Förster theory [23] is applicable, leading to a picture of incoherent hopping of excitation between chromophores. In the opposite limit of strong inter-chromophore interaction and weak coupling to the environment, the excitation dynamics can be described by quantum master equation, either the Redfield equation or, employing a secular approximation, the Lindblad equation [24]. For nonperturbative parameter ranges, more advanced methods have been developed [25,26,27], usually at the expense of higher computational complexity.
Our main optimality criterion, namely the transfer efficiency, is relatively insensitive to details of excitation time-dependence. Thus, we are going to use the simplest description of EET with the Lindblad equation, which retains coherent nature of transport, while still taking environmental interactions into account. We expect that more exact descriptions, which in general enhance excitonic oscillations, lead to similar results due to our averaging procedure. Also, these oscillations appear on the time scales of few 100 fs, which is much shorter than the time scale of excitonic transfer.
In the following, we will introduce the Lindblad equation for the overlapping disc model, used for the description of excitation dynamics in PPCs. Optimality criteria for the efficiency of PPC configurations as used in latter sections will also be presented.

Lindblad master equation
Internal dynamics of the system of N chromophores within a single-excitation manifold is determined by the Hamiltonian of the form [28] where a state |n represents an excitation on the n-th chromophore site, i.e., the electronic state of the n-th chromophore being in the 1st excited state. Because EET is sufficiently fast, events with two excitations being present at the same time are rare and it is sufficient to consider only zero and single-excitation subspace [9]. The coupling V mn is due to dipole-dipole interaction between chromophores of the form where r mn = x m − x n is a vector, connecting the m-th and n-th chromophores, d n is a transition dipole moment between the ground and the 1st excited state of the n-th chromophore.
Because the system of chromophores is coupled to the protein and nuclear degrees of freedom it is described by a reduced density matrix ρ. Decoherence due to environmental interaction, recombination of excitation to the ground state, and transfer of excitation to the sink, are modeled by Lindblad superoperators that augment the von Neumann equation for the time evolution of density matrix, To model the effects of the environment, we have taken a simplified picture of purely dephasing Lindblad superoperators (i.e. Haken-Strobl model), which is believed to capture the basic environmental effects and was used in various previous studies [29,11,13,12]. For longer time, relevant for the efficiency of PPC, it has been shown that the description with the Lindblad equation accounts for the main features of the dynamics [30,31]. Dephasing Lindblad superoperator destroys a phase coherence of any coherent superposition of excitations at different chromophores, and is given by where a site-independent dephasing rate is given by γ and { , } represents the anticommutator. Irreversible transfer of excitation from the N -th chromophore to the sink |s is modeled by Lindblad superoperator where κ denotes the sink rate. The irreversible loss of excitation due to recombination is given by an analogous term with a site-independent recombination rate given by Γ. The ground state of a chromophore system (state without any excitation) is represented as |0 . Note that L sink and L recomb can be equivalently represented by an antihermitian Hamiltonian at the expense of non-conserved density matrix probability [12], avoiding the need for an additional sink and the ground state. Relevant environmental parameters going into Lindblad equation (3) are dephasing strength γ, recombination rate Γ and the sink rate κ. We use standard values inferred from experiments and used before [12], sink rate κ = 1 ps −1 , recombination rate Γ = 1 ns −1 and dephasing rate at room temperature γ = 300 cm −1 . Dephasing rate, being a product of temperature and the derivative of the spectral density, can be estimated by using experimentally determined parameters of the spectral density (as done e.g. in reference [12]). This value approximately agrees with the optimal dephasing rate at which transfer is most efficient [12,11]. We note that the results shown depend very weakly on the actual value of the dephasing rate as long as it is of the same order of magnitude as γ = 300 cm −1 . In Appendix A we show that the values of γ = 150 cm −1 and 600 cm −1 give almost the same optimal PPC size.

Overlapping discs model
Because we want to study the dependence of efficency on the size of PPC, keeping the number of pigments fixed, we have to account for the size-dependence of the Hamiltonian. On-site energies and interaction strengths in equation (1) will be determined from the geometry of pigments configuration. Changing the PPC's size two gross effects are at play. First, as the distance between pigments is reduced, the dipoledipole interaction between pigments gets larger, enhancing the transfer of excitation among them; this effect is already taken into account by the ∼ 1/r 3 dependence of V nm in equation (2). Secondly, as chromophores get even closer together, approaching the distances comparable to the extent of the chromophore electronic orbitals, the Hamiltonian (1) is not sufficient for the description of EET anymore because the effects due to electronic orbital overlap become important. A detailed analysis of processes that take place as chromophores get close together would require advanced quantum chemistry methods and is out of scope of this paper. However, the main effect can be effectively taken into account by appropriately rescaling parameters of equation (1). Because the pigment molecules will be deformed, their excitation energies will also change. As the on-site energies n , being of the order of few eV, are about ∼ 100 times larger than V nm , even a small relative change in n can have a large effect. Effectively, close or even overlapping chromophores will therefore result in widely different values of on-site energies n at different sites, i.e., in a disorder. Thus there are two competing factors that determine the optimal size of PPCs: reducing inter- The values entering the Lindblad equation (3) should be in units of frequency, e.g. s −1 . The conversion from inverse centimeters cm −1 as used traditionally in spectroscopy is given by chromophore distance increases EET, while the overlapping of chromophores introduces a disorder that effectively suppresses EET. On-site energy has an additional contribution due to local environment (e.g. because of pigment-protein interactions), which is however of the order of ∼ 10 2 cm −1 and is usually much smaller that the on-site disorder due to pigment overlap, which is proportional to the unperturbed on-site energy of ∼ 10 4 cm −1 . Therefore, the effects of pigment-protein interaction were neglected when obtaining the results presented in the main text. To verify whether neglecting of on-site disorder due to protein interactions is justifiable, we have also calculated the optimal size for N = 7 pigments with random onsite disorder added to each random sample of chromophores. The results (see Appendix B) show that disorder of such magnitudes indeed has no gross effect on the results.
Each chromophore in our model is represented as a disc -a thin cylinder -of radius r and height a. Each disc is supposed to represent an approximate size of the electronic cloud of the orbitals involved in the EET (highest occupied electronic orbital, lowest unoccupied electronic orbital). We have estimated the height to be a = 1 Å, while the cylinder radius is taken as r = 4 Å. Size of this cylinder in comparison to a BChl pigment molecule can be seen in figure 1a. Radius r = 4 Å is chosen so that it contains 16 closest non-hydrogen atoms to the Mg atom located in the center of the pyridine ring. For given locations and orientations of discs, we then determine if there are any overlaps between discs. If two discs overlap, for an example see figure 1b, we rescale the radius of one of them to the new radius r n so that they do not overlap anymore but instead only touch. After eliminating all overlaps we end up with disc's radii r n , for an example see figure 1c.
Provided the radius of the n-th disc r n is different from the non-overlapping size r, we have to appropriately rescale the on-site energy n . If the effective size of the electronic cloud is reduced from r to r n , the kinetic energy of electron increases by a factor r 2 /r 2 n . We therefore estimate that the energy of excitation on a resized disc will also scale quadratically with its size, giving the on-site energy dependence where (0) is the excitation energy of non-deformed pigments, i.e., the energy difference of two lowest electronic states on a pigment. In FMO (0) is approximately (0) ≈ 12 300 cm −1 . The overall offset of on-site energies is irrelevant for the dynamics in the model, therefore we shift all energies by (0) . Such quadratic on-site energy scaling can be rigorously shown under an assumption that the electronic eigenfunctions of the rescaled pigment are just the rescaled eigenfunctions of the original pigment of radius r.
Let orbitals ψ i be the eigenfunctions of the Hamiltonian is kinetic energy operator and U (x) is a confining potential. The on-site energy of a given chromophore n is the difference between the energy of ground and 1st excited state, n = E 2 − E 1 . Assuming that eigenstates ψ i are just scaled to a smaller volume, ψ * i (x) = λ 3/2 ψ i (x/λ), the scaling of on-site energies from equation (7) is obtained by comparison of eigenvalue equations for the original eigenstate Hψ i = E i ψ i and the scaled eigenstate Dipole strength of the chromophores is similarly scaled linearly with the radius r n of the cylinder, where d is the bare transition dipole moment of the original chromophore of size r, and d n is the scaled dipole strength of the resized disc. This can be justified on the same grounds as the scaling of on-site energies, by inserting the rescaled wavefunction ψ * i (x) into the expression for transition dipole matrix of relevant chromophore transition, There are different possibilities of how to precisely resize discs in order to avoid overlaps. While different procedures lead to different on-site energies, the determined optimal complex size changes by little. Results in the main text were obtained by sequentially inspecting each pair of discs, resizing only the disc having greater radius afterwards, while keeping the other disc intact. We have verified other resizing procedures, for instance, resizing both discs in pair to the same size. Such resizing effectively reduces disorder of on-site energies as even strongly overlapping pigments will have identical on-site energies. Nevertheless, for such resizing procedure, the determined optimal radius of PPCs are within 2 Å of the values obtained by the resizing procedure used throughout the paper, and are thus within the error bounds of the model.
To summarize, in our overlapping disc model the matrix elements of H are calculated for given PPC configuration (positions, as well as disc and dipole orientations) by first resizing all overlapping discs, obtaining rescaled radii r n and then scaling dipole strengths and on-site energies according to equations (7) and (8).

Optimality criteria
We have already introduced equations of motion that govern the dynamics of excitation on chromophores, as well the overlapping disc model that allows us to determine the Hamiltonian for a given configuration of chromophores. What is left are criteria that will enable us to determine whether a given configuration of chromophores is efficient in terms of EET. The efficiency of the PPC complex is characterized by the probability that the excitation, initially localized on the input site, will be funneled to the reaction center trough the output site. For an example of time evolution see Appendix C. The efficiency in the model is not unity because the excitation can be lost. The probability that the excitation will be transported to the reaction center can be expressed as which will be used as our main efficiency criterion. Closely related is the average transfer time, which signifies the speed of transfer of excitation to the reaction center, and is expressed as with smaller transfer times being better.
As an additional viability criterion of PPC, robustness of efficiency to static disorder will be also inspected. Dynamic disorder due to thermal motion is already effectively described by the dephasing Lindblad terms in equation (3). Static disorder due to structural changes of PPC, for instance due to changes in biological environment, like temperature, electric charges, etc., should be treated separately. A given configuration of pigments in PPC is robust to the static disorder if random displacements of pigments from the original locations do not induce large changes in PPC's efficiency η (or equivalently, the average transfer time τ ). To put it on a more quantitative ground, we define the pigment configuration robustness σ η (x) for a given configuration of pigments ¶ with positions x = (x 1 , x 2 , ..., x N ), as a standard deviation of efficiency η when pigment coordinates are varied in the neighborhood of original positions, Probability density w(y) defines the neighborhood of a given configuration, and is localized around the original location of the pigments. The most straight-forward ¶ We omitted the disc and dipole orientations from the definition of pigment configuration robustness to simplify the expressions. However, no qualitative differences are to be expected if orientations are also varied when inspecting the robustness. choice for the distribution w(y) is a product of uncorrelated normal distributions at each pigment location, where σ defines the size of neighborhood in which the robustness is being probed. With given probability distribution, the robustness σ η is a function of original pigment locations x and size of deviations from original locations σ. In the limit of small deviations, σ → 0, the expression can be simplified to where the sum goes over all components of pigment coordinates. The robustness σ η in the limit of small pigment displacements is thus proportional to the amplitude of pigment displacements.

Optimal number of pigments: the case of Fenna-Matthews-Olson complex
In this section, the efficiency and robustness of random configurations within the model will be considered on an example of FMO complex. FMO consists of N = 7 BChl pigments. On-site energy for BChl was chosen to be (0) = 12 300 cm −1 , which is within the range of on-site energies for BChls in FMO as determined in the literature [32,33,34]. The strength of transition dipole moment d = |d| was taken as d 2 = 26 D 2 + (note that published values for d from calculations and experimental data vary considerably [35]). On-site energy (0) and dipole strength d used hold for BChl pigments in general and therefore the results presented are expected to be valid also for other PPCs containing BChls, not just for the FMO complex.
To determine the optimal size of PPC for a given number of chromophores (or equivalently, optimal number of chromophores for a given size), we considered two criteria based on overall behavior of efficiencies and robustness of random configurations. In the following, the motivation for choosing such optimality criteria will be given, and the results for the case of FMO will be presented.

Average efficiency
We shall use the average efficiency η , averaged over random positions and orientations of pigments enclosed in a predefined volume. The reason to use random averaging with uniform distribution is twofold: first, high average EET efficiency under uniform averaging will mean that there are many different configurations that have high efficiency, i.e., high efficiency is globally robust. Choosing averaging over random configurations therefore offers insight in how special efficient configurations of chromophores are within the space of all configurations. Second reason is that we have a priori no knowledge what would be the appropriate measure for possible pigment configurations under say different protein cage foldings due to for instance mutations. A uniform measure represents in this case a "least-information" distribution. * Using configurations sampled according to a uniform distribution over chromophore positions within a ball of radius R and random orientations of dipoles and discs, the average efficiency is calculated. Formally, it can be written as where X contains positions and orientations of chromophore discs and dipoles (apart from the positions of input and output sites which are fixed on the poles of the sphere), and w conf (X) ∝ 1 signifies a uniform distribution of chromophores inside the sphere.
Observing the dependence of the average efficiency η(R) for different number of chromophores and different radii R of the enclosing sphere, we can determine the optimal number of chromophores for a given radius R, or equivalently, the optimal radius R opt for a given number of chromophores.
To obtain a more detailed information about efficiencies of random configuration, we also observed the probability distribution over efficiencies p R,N (η), defined as p R,N (η) = R δ(η(X) − η)w conf (X)dX. For the number of chromophores as found in FMO (N = 7), the probability distribution over efficiencies p R,7 (η) is shown in figure 2. When going from large radii R to smaller, configurations tend to get more efficient, which is expected as the chromophores are closer to each other, thus increasing the dipole coupling. However, as R is reduced even further, overlapping of chromophores gets more probable, causing an on-site energy disorder. This leads to the localization of excitation on chromophores not connected to the sink site. Such configurations have low efficiency. Therefore, as R gets smaller the distribution p R,N becomes bimodal, with lower efficiency mode due to overlapping configurations and high efficiency mode for non-overlapping configurations.
Low efficiency of overlapping configurations therefore leads to a maximum in the average transfer efficiency η(R) at the optimal radius R opt . The average transfer efficiencies at the optimal radius are rather high, e.g. for N = 7 in figure 2 it is * We have to note that we also checked other distributions, for instance a uniform distribution on the surface of a sphere of radius R, and obtained practically the same results. For instance, the difference in the position of the maximum in figure 2 was within our error estimate of 1 Å (seen as an "error" band in figure 3). η(R opt ) ≈ 0.95 with a large fraction of configurations having even larger efficiency than the average. Thus within the model, high efficiency is not due to finely tuned pigment positions and orientations, but occurs for majority of pigment configurations for parameters estimated to be relevant in PPCs. For the FMO case with N = 7, the optimal radius was estimated to R opt ≈ 16 Å, which fits the actual configuration of pigments very well (see figure 3). The average transfer time τ is also minimal at R = R opt (see figure 2). Optimal average transfer time of ∼ 30 ps is so large due to the contribution of very inefficient configurations of chromophores. Looking at the average transfer time of the 5% of most efficient configurations, we get the value of 5 ps, which is comparable to the transfer times as determined using different models of the FMO in references [12,13,14,36].
The estimated optimal radius is quite insensitive to small variations of input parameters, e.g. dipole moment d, chromophore disc radius r or its thickness a, or the scaling of on-site energies and dipole strengths of resized discs. For instance, decreasing disc radius to r = 3.5 Å decreases R opt by ≈ 2 Å, changing disc thickness to a = 0.5 Å or 1.5 Å changes R opt by ≈ ∓2 Å, while changing quadratic energy scaling to a linear or cubic one again changes R opt by ≈ ∓2 Å. Similarly, changing the dephasing rate γ by a factor of 2 changes R opt by ≈ 2 Å, see Appendix A. Details of the disc resizing procedure also change R opt for less than 2 Å, as the extreme case of resizing each overlapping disc pair to the same size reduces R opt by ≈ −2 Å.
The optimal radius of the enclosing sphere was obtained from the average efficiency over all random configurations within a sphere. However, even though the evolutionary drive to more efficient configurations might not be very strong if majority of configurations are already efficient, still some optimization is to be expected. Thus one might argue that the optimal enclosing volume of the natural PPCs should be determined considering only the ensemble of more optimal configurations. We will denote such averages with η p where p specifies a portion of most efficient configurations that should be taken into account when calculating the average (e.g. η 0.05 is the average of η over 5% of most efficient configurations as shown in figure 4). As p is reduced, the overall value of average efficiency η p will increase. The increase will be more pronounced in the region of R < R opt , where the distribution is bimodal. The location of the maximum of η p will be thus moved to smaller values of R, indicating more densely packed chromophore configurations. However, as we will see in next subsection, robustness of such densely packed configurations deteriorates very quickly, supporting our choice of estimator for the R opt .
Note that the overlaps between pigments and protein cage are not considered explicitly in the model. If overlaps with protein cage would be taken into account, R opt would represent the size of a protein cage, whereas in our model without pigmentprotein overlaps, R opt is the size of a sphere that contains all pigment centers. For instance, looking at figure 3b, we can see that the sphere with R opt contains all pigment centers, while parts of few pigments still protrude the bounding sphere. If overlaps of pigments with the protein cage would be taken into account explicitly, R opt would be approximately by a disc radius r = 4 Å larger, i.e. in corresponding figure, the bounding sphere would enclose all pigments completely.

Robustness
Robustness of PPC configurations to static disorder should also be taken into account when determining whether a given configuration of pigments is feasible, as the conditions in which PPCs operate are subject to constant environmental changes. In previous subsection, we have inspected the probability distribution of efficiencies η over random configurations, showing that majority of random configurations achieve relatively high efficiency when the enclosing volume is optimal. In this subsection, we will present analogous analysis of the robustness of random configurations, in particular of those with high η. We shall show that highly efficient configurations in small enclosing R are very fragile.
We have defined robustness of efficiency σ η in equation (11). In simulations we have displaced the pigment positions according to normal distribution with a width of σ = 0.1 Å, which is small enough to quantify the robustness in the neighborhood of specific configuration, while larger than displacements due to thermal vibrations, which are already effectively described by the Lindblad equation. We are specifically focusing on a subset of the most optimal configurations in terms of η. The average robustness of the subset of optimal configuration is denoted as σ η p , where p specifies a fraction of most optimal configurations in terms of EET efficiency η. As an example, we will consider robustness of top 5% of efficient configurations σ η 0.05 . The dependence of the average robustness on the radius of the enclosing sphere R is shown in figure 4. The average efficiency of optimal configurations η 0.05 is also shown in the figure. While the average efficiency of top 5% of optimal configurations continues to rise as the enclosing sphere radius R is reduced, we can see that the average robustness σ η worsens very quickly as the R drops below R opt .
Quick worsening of EET robustness with reducing sphere radius suggests that even if the PPC configurations occurring in nature are indeed optimized in terms of pigment positions and orientations, the excessive stacking of pigments is not favored as it makes PPC configurations very sensitive to any displacements of pigments. The transition from robust to non-robust regime takes place at a radius comparable to R opt at which the average efficiency η has a maximum. This is not surprising as both, the efficiency of configurations and robustness of configurations, are strongly influenced by the overlapping of pigments which gets more pronounced for R R opt .
We have presented results for the robustness σ η of the 5% of most efficient configurations, with pigment displacements σ = 0.1 Å. General characteristics of σ η p however do not quantitatively change for different p (the case of p = 0.15 is also shown in figure 4) or displacements σ. Most importantly, the radius R at which the robustness of configurations drops significantly takes place at approximately the value of R opt . Same behavior of robustness is observed also in the limit of infinitesimal robustness from equation (13) where σ → 0.

Optimal pigment numbers in other PPCs
In previous section we have calculated the optimal size R or the optimal number of pigments for the FMO complex. We also demonstrated that a large portion of chromophore configurations has high efficiency when the enclosing volume is properly chosen (∼ R opt ). Additionally, robustness of configurations to chromophore displacements starts to deteriorate quickly once the enclosing volume is reduced below R opt . Based on these two observations, we argue that the enclosing volume of PPCs occurring in nature should be close to the optimal volume as determined by our simple model. In this section we will present similar results for the PPCs containing chlorophyll (Chl) chromophores.
We compare results of the model to the structure of PPCs from the Photosystem II (PSII) [37], found in cyanobacteria, algae and plants. PSII consists of multiple functional units, which are either part of the outer light-harvesting antenna or the inner core, to which excitations are funneled. In the light-harvesting antenna we will consider the light harvesting complex II (LHCII), while in the core we will focus on the PC43 and PC47 complexes that funnel excitations to the reaction center and thus have a similar role as the FMO complex in bacteria. A monomeric unit of LHCII contains 14 chlorophyll molecules (8 Chl-a and 6 Chl-b), while CP43 and CP47 contain 13 and 16 Chls respectively.
Model parameters for the sink rate, dephasing and recombination are kept the same as in the FMO case, while the transition dipole strength and on-site energies are different for Chl molecules. Transition dipole moment of Chl molecules is chosen as d 2 = 15 D 2 and the on-site energy (0) = 15 300 cm −1 , where values were taken according to reference [38] (we take the average between values for Chl-a and Chlb). For the CP43 and CP47 complexes we have simulated random configurations of 13 and 16 chromophores enclosed into a sphere as the actual chromophore positions are distributed relatively uniformly in all directions. The shape of LHCII is however significantly elongated in one direction. We therefore choose the cylindrical volume, having only one additional parameter that has to be provided, i.e. the ratio between the cylinder radius R c and cylinder height A. Based on positions of the LHCII chromophores we have estimated the ratio of the two to be R c /A = 0.34.
The CP43 and CP47 primarily play a role of an exciton wire, making the model with input site and output site at the opposite sides of the sphere applicable. The optimal radius as predicted by the model is R ≈ 18 Å for CP43 and R ≈ 20 Å for CP47. As the LHCII also has to transport excitations from adjacent complexes, we have also determined the optimal shape of LHCII with input and output sites located at the opposite sides of the enclosing cylinder. With the ratio R c /A fixed, we have varied the height A of the cylinder and determined the optimal height to be A opt ≈ 43 Å.
In addition to acting as an excitation wire, CP43 and CP47 complexes are also directly involved in the absorption of photons, in which case the role of the input site can be taken by any chromophore site. This is even more common scenario in the LHCII complex, whose primary role is the absorption of photons. To verify whether the findings about optimal enclosing volume are also valid when the main purpose of PPC is the absorption of photons, we randomly placed the input site inside the enclosing geometry for each configuration in random ensemble. General characteristics of the distribution over efficiencies p R (η) do not change considerably, however, the distribution is somewhat shifted to the higher efficiencies because in many random configurations the input site is considerably closer to the output site than the diameter of the enclosing volume. This results in the optimal size of enclosing volume being somewhat larger, R opt ≈ 20Å for CP43 and R opt ≈ 22Å for CP47. For the LHCII we have moved the output site to the midpoint on the side between top and bottom of the cylinder, where the actual output site is supposedly located [41]. For such geometry and previously used R c /A = 0.34, we have obtained the optimal height of the enclosing cylinder at A ≈ 47Å.
The optimal enclosing values obtained from the model (averaged between the case for fixed input site and random input site) were compared to the actual configurations of pigments as obtained from spectroscopic data, and are shown in figure 4a-c, showing good agreement. For the spherical geometries, we have centered the sphere of optimal radius R opt to the arithmetic mean of locations of BChl/Chl centers. For the LHCII, where cylinder was used as the enclosing geometry, the cylinder axis was determined such that N i r 2 ⊥i was minimal, where r ⊥i is the distance from the cylinder axis to the position of i-th Chl center. Interestingly, three cylinder axes do not lie in a plane but are instead tilted at an angle 15 • to the plane containing three cylinder centers. It is not known if this plays any functional role.

Conclusion
We have studied the efficiency of excitation energy transfer in protein-pigment complexes for random configurations of pigments. The Hamiltonian part of Lindblad master equation is determined from the geometry of pigment configurations. If pigments are too close, so that they overlap, this introduces a disorder in on-site energies, effectively inhibiting excitation transport. Fixing the enclosing volume in which pigment molecules are located we have calculated the average efficiency over random pigment configurations as well as robustness of efficiency to variations of pigment locations. Doing this we have determined the optimal number of pigments for a given size of the complex. Even though the model is an oversimplification of actual processes that take place in nature, statistical predictions obtained from the model are robust to its variations.
Comparing theoretically predicted optimal number of pigments with several naturally-occurring complexes we find good agreement. This might indicate that PPCs are not optimized just to have the highest possible efficiency -in fact, efficient configurations are quite common -but instead to be robust to variations in pigment locations. Namely, it turns out that configurations optimized for the highest efficiency, that is those with specially tuned positions and dipole orientations, are very sensitive to small perturbations. The number of pigments in nature is therefore chosen in such a way that the probability of having efficient configurations that are at the same time also robust is the highest.
The presented findings could be in principle verified experimentally by modifying the structure of known PPCs and probing the efficiency of excitation transfer. For the FMO complex the structure was already changed by mutation of genes encoding the structure of BChls, as well as by substituting the carbon 12 C atoms with 13 C [22]. Comparison of excitonic spectra revealed no distinctive differences in the dynamics of excitations, complying with the hypothesis that configurations are not highly tuned but instead very robust. An additional intriguing possibility would be also to inspect the characteristics of FMO with mutated protein cage, modifying positions and orientations of pigment molecules. One could also compare our predictions for the optimal sizes (e.g. figure 3a) with other complexes occurring in nature.

Appendix A. Dependence on the dephasing rate γ
In the simplified model used, environmental interaction is described by dephasing rate γ, having the same value for all chromophore sites. In principle, environmental interaction requires more involved equations of motion (e.g. HEOM [7]), taking into account the spectral density of the environmental modes. However, due to crude nature of the model, simplified description is expected to account for main environmental effects influencing the efficiency of exitation transfer (see e.g. [30,31] for more detailed comparison of approaches). The adequacy of simple Lindblad-type description of dynamics for determination of optimal size is also justified due to the high robustness of the results to the actual choice of dephasing value γ, as seen in figure A1a, where R opt only changes for ∼ ±2Å as dephasing rate γ is changed by a factor of 2. Optimal size R opt of the PPC is somewhat smaller as dephasing rate gets stronger, which is expected as larger dephasing rate enables transfer across the sites with greater on-site energy mismatch, getting increasingly common in more compact configurations of chromophores.

Appendix B. Random on-site disorder
To verify whether effects of the local chromophore environment due to e.g. pigmentprotein interaction can affect the findings about the optimal PPCs sizes, we have amended the Hamiltonian in equation (1) with random on-site disorder rand n of the magnitudes as found in naturally occurring PPCs (i.e. on-site energy differences in the order of ∼ 100 cm −1 ). The values of disorder for each realization of random PPC were calculated according to Gaussian distribution with variance σ rand . Results are shown in figure A1b. In the region R < R opt , where average transfer efficiency is strongly affected by disc overlaps, an addition of random on-site energy disorder has no noticeable effect. The effect is more pronounced for R > R opt where overlapping of discs is not the limiting factor of transfer efficiency any more. The estimated optimal size R opt however is not changed considerably by an addition of random on-site disorder.

Appendix C. Time evolution of site populations
To provide some insight into the temporal dynamics of the excitation transport, we present the time evolution of site populations for two different realizations of PPC within the Lindblad model. In figure C1a, time evolution for the FMO Hamiltonian from reference [6] is shown, and in figure C1b, the time evolution of randomly generated PPC of radius R = 18 Å. The values of dephasing, sink rate and recombination rate are the same as used in the main text.