{'fraction_non_alphanumeric': 0.028619703665193657, 'fraction_numerical': 0.008708084221471277, 'mean_word_length': 4.585220673635308, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimization of energy transport in the Fenna-Matthews-Olson complex via site-varying pigment-protein interactions

Energy transport in photosynthetic systems can be tremendously efficient. In particular we study exciton transport in the Fenna-Mathews-Olsen (FMO) complex found in green sulphur bacteria. The exciton dynamics and energy transfer efficiency is dependent upon the interaction with the system environment. Based upon realistic, site-dependent, models of the system-bath coupling, we show that this interaction is highly optimised in the case of FMO. Furthermore we identify two transport pathways and note that one is dominated by coherent dynamics and the other by classical energy dissipation. In particular we note a strong correlation between energy transport efficiency and coherence for exciton transfer from bacteriochlorophyll (BChl) 8 to BChl 4. The existence of two clear pathways and the role played by BChl 4 also challenges assumptions around the coupling of the FMO complex to the reaction centre.


I. INTRODUCTION
Given the importance of the effect of the environment upon energy transport in photosynthetic systems, it is common practice in theoretical studies to assign to each chromophore, identical spectral densities, usually described by a smooth, simple functional form. While this simplistic approach should still allow the essential physics to be captured, it may nevertheless mask the influence of realistic site variations which could potentially contain interesting physics in itself.
In fact, recent computational chemistry calculations have shown, for instance in the case of the Fenna-Matthews-Olson (FMO) complex of the green sulfur bacteria, that site variations of spectral density can be significant [1][2][3] . Obviously, in biological systems, traits showing some degree of heterogeneity are expected, and in many cases, it is simply random. However, the green sulfur bacteria, due to its ability to survive in extremely low-light conditions, must certainly possess highly-evolved light-harvesting machinery. It is then only natural to suspect that there may be an underlying design principle associated with this heterogeneous system-bath couplings. A previous publication demonstrated that the heterogeneous spectral densities can significantly impact exciton transport dynamics and suggested that this could be a mechanism used to tune energy transport 2 .
In addition, exciton delocalization may also be affected by site-dependent system-bath couplings. Sato a) Electronic mail: coker@bu.edu b) Electronic mail: david.hutchinson@otago.ac.nz and Reynolds 4 found that the length and robustness of quantum coherence in a generic dimer is simultaneously increased under specific ratios of site energy mismatch to site reorganization energy mismatch. Due to the presence of long-lived coherence within the timescale of exciton dynamics and the ability of quantum walks to sample multiple paths simultaneously, there is a large amount of speculation that quantum coherence may remain important to the remarkable efficiency of photosynthetic energy transport. It is generally accepted that the dynamics are modified by coherence and that optimal transport occurs in the intermediate regime between the coherent and incoherent limits [5][6][7][8][9][10] . However, it remains an open question whether or not, in this optimal regime, comparable or even better efficiencies can be achieved by a more classical type of mechanism, e.g. energy funneling. As such, it is entirely plausible that the presence of quantum coherence is to serve a different biological function or that it is simply a by-product of dense chromophore packing, with no specific function of its own. Finding evidence of underlying coherence optimization and a system design which can effectively utilize coherence would certainly provide some confidence regarding its importance in energy transport.
In this work, we present evidence that the sitevarying spectral densities of the FMO complex are indeed optimized for energy transport. This finding has been made possible by the availability of accurate sitedependent spectral densities from recent quantum chemistry/molecular dynamics calculations 1,11 .
The basic theories of energy transport or spectroscopic response assume that the full Hamiltonian can be broken into a system part,Ĥ S , a bath, or environmental component,Ĥ B , and an interaction between these system and bath subsystems,V . The partitioning of the different degrees of freedom between the system and bath subsystems is determined by the resolution capabilities of a given experiment. Thus, for example, high resolution 2D electronic spectroscopy (2DES) can actually resolve transitions between higher frequency intramolecular vibrational and electronic states, so in such experiments the system part is best described by the vibronic eigenstates of the chromophores and the bath becomes the lower frequency environmental, or intermolecular degrees of freedom. For less well resolved experiments, for example, monitoring energy migration from one chromophore to another and not resolving (i.e. averaging over) the intramolecular vibrational structure of the chromophores, the system is best described by their coupled electronic states, while the bath becomes all the intra-and intermolecular vibrations. The study presented here approximately adopts this latter perspective in that it assumes thatĤ S is represented by the eigenstates of the coupled electronic subsystem, while the intermolecular spectral density composed of a continuum of protein environmental modes describes the bath. In these studies we will disregard the influence of the discrete intramolecular component of the spectral density as they are expected to play a fairly minor role in energy transport in chlorophyll systems with small intramolecular Huang-Rhys factors. The possibility of further optimization from vibronic contribution and interplay with the intermolecular component of the spectral density will be explored in a future publication. Preliminary results from these further studies do not materially affect the conclusions in this paper.

II. THEORETICAL MODEL
The total system-bath Hamiltonian is expressed aŝ whereĤ s ,Ĥ b andĤ sb are the system (electronic), bath and system-bath interaction Hamiltonians respectively.
A. The electronic system The electronic system consists of the chromophores in the pigment-protein complex (PPC). Its Hamiltonian governs the coherent part of the evolution and is described by the tight-binding model in the site basis |m : where E m , V mn and N are the site energies, electronic coupling between pigments m and n and number of chromophores respectively. Since our system operates in a low-photon environment, the single exciton manifold approximation is valid here. The kth eigenstate of the system Hamiltonian (also known as the exciton state) can be decomposed in terms of the site basis: In this work, we use the 8-site Hamiltonian for the Prosthecochloris aestuarii (P. aestuarii ) species as presented by Moix et. al 12 .
B. The bath and system-bath coupling The protein environment is commonly modelled as a bath of harmonic oscillators withĤ are the creation (annihilation) operators of excitations of the ith bath mode with frequency ω i,m on pigment m. In this work, the phonon modes on each site are treated as being uncorrelated with each other and coupled linearly to the diagonal part of the system Hamiltonian such that: Here u i,m is the coupling between the electronic transition of the mth site and the ith phonon mode. Physically, this equation reflects the effect of the protein environment dynamically modulating the site energies of the pigments. All information about the system-bath interaction of each pigment m is contained in its corresponding spectral density J m (ω) = π i |u i,m | 2 δ(ω−ω i,m ). The reorganization energy λ m of pigment m is in turn related to J m (ω) by the following integration over frequency ω:

III. NUMERICAL METHODS
The numerical computation of the dynamics was performed using the Modified Redfield Theory (MRT) 13 and its more recent upgrade, the Coherent Modified Redfield Theory (CMRT) [14][15][16][17] . It has been shown that MRT and CMRT are reasonably valid over a broad range of system-bath coupling strengths and provide reasonable agreement with dynamics computed from numerically exact methods 14,16,18 . This, coupled with the convenience and efficiency as compared to numerically exact methods make them suitable numerical methods in this work.
In contrast to MRT, CMRT allows the computation of coherence terms in the density matrix as well as the incorporation of some non-Markovianity. Also while the MRT population vector P (t) can only be in the exciton basis, the CMRT density matrix ρ(t) can be obtained in both the site and exciton basis, with the site basis being the default (Appendix C). The downside of CMRT however, and the reason why both versions of the same method had to be used, is its much longer computation time. This renders it infeasible for the statistical methods and genetic algorithm used in Sections V A and V B respectively, which incidentally require only the population terms. We thus employ MRT in those cases while using CMRT for the rest of the work in this paper. As in Refs. 16 and 17, we solve for the non-Markovian dynamics of CMRT with the non-Markovian Quantum Jump (NMQJ) technique 19,20 .
To justify the use of MRT in place of CMRT, we have checked that the effects of non-Markovianity are not overly significant, or in other words, the general trend of population dynamics (in the exciton basis) is similar with the two versions. For self-consistency, we give a brief outline of MRT, CMRT and NMQJ in Appendices A, B and C respectively. The full derivation and additional details can be found in Refs. 13 and 18 for MRT, Refs. 14-17 for CMRT and Refs. 19 and 20 for NMQJ. All computation is performed for the physiological temperature of T = 300K.

IV. NUMERICAL SETUP
The two most important pigments in the FMO are the initial excitation site and the target site -the former influences the exciton dynamics and the latter determines our assessment of energy transport efficiency. It is now believed that BChl 8 is the most likely linker site between the chlorosome and the rest of the FMO complex 21 . However, BChl 8 is normally lost during sample preparation, and before its recent discovery, 22,23 BChls 1 and 6 were proposed as the possible linker sites. These two pigments are also the usual initial photoexcitation sites in spectroscopic experiments. For these reasons, many theoretical studies of the past and present employ BChls 1 and 6 as the initial excitation site. Since our motivation is to understand the in vivo workings of the FMO, and not for comparison to spectroscopic data, we choose BChl 8 for this work.
In the literature, the linker site to the reaction centre is usually assumed to be BChl 3, which is also the lowest energy pigment with the closest proximity to the reaction centre. There exists some ambiguity, however. Wen et. al. 24 reported that it is the BChl 3 side of the FMO complex which interacts with the reaction centre, but their experiment did not pinpoint the exact pigment(s). Furthermore, it has been mentioned in several publications 7,25-29 that the two lowest energy pigments, BChls 3 and 4, are in the target region close to or in contact with the reaction centre. As such, we consider both BChls 3 and 4 as the possible target sites in this paper.
To aid physical interpretation, we represent the realistic spectral densities in the Drude-Lorentz regularized Ohmic form 30 : where λ m and Ω m are the reorganization energy and cutoff frequency of site m respectively. Here the values of λ m are calculated from the realistic site spectral densities via eq 5. Ω m , which corresponds to the inverse of the bath correlation time, can be obtained by first computing the site-dependent bath correlation function (eq. A6). For simplicity, we assume that the real part of the bath correlation functions can be represented by time decaying exponential functions, from which it is possible to extract the bath correlation time and subsequently Ω m . We have verified that the dynamics with the realistic sitedependent spectral densities can be reliably reproduced with their Drude-Lorentz forms.
The computed values of λ m and Ω m are presented in Table I. It is clear that the FMO complex has a fairly significant range of λ m values where the largest is more than 2.5 times the magnitude of the smallest. It is interesting to note that the mean reorganization energy is also in a sense the most representative since half of the pigments have λ m values close to this value.
The complexity of the problem is now significantly reduced since the spectral density is characterized by only two physically meaningful bath parameters, i.e. λ m and Ω m . It turns out the problem can be further simplified. By comparing the dynamics obtained with site variation in only one of the bath parameters to that obtained with a completely site-independent benchmark, we established that it is primarily the site variation in λ m which modifies the dynamics, with very little contribution from the Ω m site variation. Here the site-independent case is constructed by assigning the mean value of λ m and/or Ω m from Table I. Importantly, we note an improvement in exciton transport to the target sites with the site-varying spectral density. Therefore, all computation will be performed using the Drude-Lorentz spectral density representation, keeping the site-dependent Ω m fixed to the original configuration values in Table I and varying only λ m .

A. Optimality of system-bath coupling configuration
Following the observation of energy transport enhancement with the site-dependent spectral density, a more rigorous assessment of optimality is necessary. To this end, the energy transfer efficiency with the FMO sitedependent λ m configuration must be compared to that from a large sample of random site-dependent λ m configurations. Here the energy transfer efficiency, as obtained using the MRT formulation, is defined as the timeaveraged population at the predominant exciton state of the target site: Here P k is the population in exciton state k, with k = 1(2) for the target site m = 3(4). The participation of an exciton state in each site is denoted by the absolute square of the corresponding expansion coefficient in eq. 3 and ranges from 0 to 1. Here |c 3,1 | 2 = 0.88 and |c 4,2 | 2 = 0.59. We have chosen τ = 1 ps since a large portion of the relevant exciton dynamics and the experimentally observed coherence occur within this timescale. Two types of randomization were performed. In the first set, the efficiencies of all possible site permutations of the original λ m configuration were computed and presented as a histogram in Figure 1.
Notice that the efficiency of the actual FMO λ m configuration falls in an impressive upper percentile range (85th and 78th percentile for the target sites of BChl 3 and 4 respectively). The fact that it is not the most efficient configuration is not a concern here. In fact, it is not surprising in the context of evolution, since biological constraints may be present and traits only need to be sufficiently functional. In contrast, the site-independent case yields an efficiency closer to the mode of the distribution. A similarly high level of optimization is also observed in our second set of randomization, where now the only constraints on the values are that the mean, maximum and minimum are similar to that of the original configuration. For consistency, the same sample size as the first set, i.e. 40320 was used. As shown in Figure 5 of Appendix E, the efficiency corresponding to the original λ m configuration is in the 87th and 81st percentile for the target sites of BChl 3 and 4 respectively. Once again, the efficiency for the site-independent case appears near the mode of the distribution. To provide a visualization of how a random unoptimized λ m configuration could be detrimental to transport efficiency, we present in Figure 2 the comparison of the exciton dynamics with the original λ m configuration to that with the least efficient configuration in Figure 1.
Finally, it is worth mentioning that when we tested BChls 1 and 6 as the initial excitation sites, no such optimization is present. Compared to the site-independent case, there is no appreciable improvement in energy transport in the case of BChl 1, while for BChl 6, there is in fact a deterioration in efficiency. Assuming the distribution of λ m is a result of evolution to optimize energy transport, this observation is further testament that BChl 8 is likely the first pigment to receive the exciton from the chlorosome.

B. Source of optimality
Having established that the actual site-dependent λ m configuration is highly optimized for energy transport to both BChls 3 and 4, we seek to identify the underlying mechanism behind its effectiveness. We note that network connectivity and site energy distribution determine how effective a particular mode of energy transport (dissipative or coherent) would be. Moreover, larger reorganization energies are conducive to energy dissipation while smaller reorganization energies are beneficial for the sustenance of quantum coherence. This suggests an optimized interplay between the FMO λ m distribution and the design of the electronic system. Therefore, to determine if the system design is optimized for a dissipative or coherent mode of transport, we establish whether the best efficiency is obtained with larger or smaller values of λ m . This is achieved through the use of a genetic algorithm (Appendix F). Here, the fitness function is the same measure of efficiency used in the previous section, i.e. the efficiency at the predominant exciton state as evaluated with the MRT formulation (eq. 7). Now, for the target site of BChl 3, the optimal configuration is one where all the sites have the maximum λ m value. Clearly, the optimal mechanism of energy transport to this target site must be a dissipative one. On the other hand, for the target site of BChl 4, the optimal values are minimized for the three lowest energy pigments (BChls 2, 3 and 4) and maximized for the remaining pigments. This suggests that a combined coherent and dissipative mechanism leads to optimal energy transport in this case. We note that the minimization for BChls 3 and 4 could partly be due to increased coherence prolonging linear combinations between two sites. This would lead to increased population in BChl 4 on average compared to if the process was dissipative, promoting only downhill energy flow to BChl 3. However, we fail to see any indication of Rabi oscillations suggesting this process is highly overdamped and negligible here. Next, we performed a second genetic algorithm run, but this time with a constraint on the mean, mimicking the conditions in Section V A. Due to this constraint, it is now possible to assess the significance of each pigment to the energy transport process. This is because the only way now for the algorithm to maximize efficiency is through prioritizing the most important pigments by assigning them the best λ m values, and leaving the remainder to the less influential pigments. Thus, in a dissipative process for example, the more influential a pigment is, the larger the computed optimal λ m would be.
The results from the second genetic algorithm run are tabulated in Table II,  The λ ga m,3 and λ ga m,4 configurations also provide useful optimality benchmarks. Through comparison with the FMO λ m configuration, the source of its high degree of optimality can be identified. With the exception of λ 1 , all the site reorganization energies show some level of optimization for either one or both of the target sites. The two smallest λ m values, λ 4 and λ 5 show a high degree of agreement with the genetic algorithm values for both target sites. λ 2 and λ 6 are optimized for only one of the target sites. Meanwhile, the values of λ 3 and λ 7 are somewhere in between the optimized values for the two target sites, thus achieving a compromise. Lastly, λ 8 , despite being close to neither of the optimized values, show a similar trend to the genetic algorithm results by virtue of its larger than average value. With regards to the mutual optimization for both target sites, three mechanisms of optimization can be identified. The first is where the magnitude of λ m mutually benefits both target sites, for example the larger than average value of λ 8 is advantageous for the dissipative type of energy transport to both BChls 3 and 4. The second mechanism, which is also the most interesting, is where the magnitude of λ m benefits only one of the target sites, while minimizing its negative effects on the other target site. This applies when a particular pigment is involved in different effective modes of transport to each of the target sites. A case in point is the small value of λ 4 , which is advantageous for the partially coherent energy transport to the target site of BChl 4. At the same time, it does not overly impede dissipative energy transport to the target site of BChl 3 since the dominant pathway is not involved. Finally, the third mechanism is one where the magnitude of λ m does not assist energy transport to any of the target sites, but the negative impact is simply minimized. For example, the small value of λ 5 is unfavourable for the dissipative energy transport to both target sites, but the negative effect is mutually minimized since BChl 5 is not on a dominant pathway for any of them.
At this juncture, we must stress that it is not the λ m site variation per se which is responsible for the enhanced performance, since it is really the magnitude of λ m that matters. Rather, it is more likely a case of Nature making the best of an unavoidable situation. Biological constraints in the FMO can lead to site variation in λ m ; for example, pigments located at the protein-solvent interface tend to have larger reorganization energies than those in the interior 2 . The high efficiency of the FMO λ m configuration compared to various other random configurations of similar average (Section V A) is then simply a consequence of the system and system-bath interaction having evolved to complement each other in a very effective manner, e.g. via selection of appropriate dominant pathways and effective transport mechanisms.

C. Relevance of quantum coherence to efficiency
Given the evidence for both dissipative and coherent energy transport pathways, we investigate how energy transport efficiency relates to coherence length in the presence of site-dependent λ m . Since the coherence terms of the density matrix are required, all the computation in this section is performed using the CMRT formulation. For the efficiency η m , we utilize a measure similar to eq. 7 but for the site basis counterpart, i.e: where ρ mn (t) are the elements of the time-dependent density matrix computed using the CMRT formulation, with the diagonal elements ρ mm (t) representing the population at site m. Here m = 3(4) for the target site of BChl 3(4). As usual, τ = 1 ps.
To quantify the degree of exciton delocalization, we use the coherence length defined by 31,32 L ρ (t) ranges from 1 for the case of zero coherence to N for a fully delocalized state, i.e. larger L ρ (t) values indicate a larger degree of delocalization. At t = 0, even though the state is fully coherent, the population is completely localized at a single site, in which case L ρ (0) = 1/N . As a result, L ρ (t) starts from 1/N in all cases, increases then peaks after a certain period of time before decreasing. This reflects the scenario of an initially localized pure state becoming more delocalized before the interaction with the environment gradually destroys the coherence.
With the original λ m configuration, L ρ (t) peaks at the value of about 2 around t = 0.7 ps, corresponding to the fact that within the chosen time scale of τ = 1 ps, the dynamics are still relatively far from thermal equilibrium. Note however, that the state of complete incoherence can never be reached even at equilibrium 33 .
In analogy to eqs. 7 and 8, we define the time-averaged coherence length as: The relationship between efficiency η m and the timeaveraged coherence length L ρ,avg is depicted as a scatter plot in Figure 3. The data points correspond to 200 random site-dependent λ m configuration plus 6 additional relevant data points, namely the original configuration, the site-independent configuration, the optimized solution from the genetic algorithm (i.e. λ ga m,3 and λ ga m,4 ) and the most efficient site permutation of the original configuration (which we shall denote as λ perm m,3 and λ perm m,4 for target sites BChls 3 and 4 respectively). For the target site of BChl 3 (Figure 3), even though there is a positive correlation between L ρ,avg and η m , it is in the weak to moderate regime, with a correlation coefficient of only 0.36. This relatively weak correlation implies that it is unlikely effects due to coherence are significant in the energy transport mechanism to BChl 3, and further confirms the mainly dissipative nature of the energy transport mechanism as demonstrated in the previous section.
We observe a remarkably strong positive correlation between η 4 and L ρ,avg , with a correlation coefficient of 0.90. We also examined the correlation with V 34 set to zero in the system Hamiltonian, in order to address the concern that the remarkable correlation could be predominantly attributed to the strong coupling between BChls 3 and 4. While the positive correlation decreased as expected, it remains strong with a correlation coefficient of  0.76. This shows that even though the strong coupling between the two pigments undoubtedly plays a role, it is not overwhelmingly responsible for the correlation between coherence and η 4 . Figure 3 also clearly showcases the capability of the FMO λ m configuration to accommodate efficient energy transport to both target sites, as we have inferred in previous sections. Even though the FMO λ m configuration does not lead to the best efficiency for either of the target sites, it nevertheless corresponds to relatively high efficiency for both target sites. In contrast, two of the most optimized configurations, λ ga m,3 (λ ga m,4 ) and λ perm m,3 (λ perm m,4 ) only produce superior efficiencies for their respective target site BChl 3(4), but are significantly less remarkable and are even inferior to the FMO configuration for the other target site BChl 4(3). This illustrates the nontriviality of navigating trade-offs to sufficiently accommodate two largely uncorrelated and partially conflicting processes. The highly coherent nature of λ ga m,4 and λ perm m,4 (as reflected by their large L ρ,avg values) while advantageous for the partially coherent energy transport to BChl 4, is unconstructive for the dissipative energy transport to BChl 3 ( Figure 3). Similarly, λ ga m,3 and λ perm m,3 which only corresponds to moderate coherence only leads to average efficiencies at BChl 4 ( Figure 3). Meanwhile, the FMO λ m configuration gives rise to a rather impressive degree of exciton delocalization, where its L ρ,avg value Lastly, we explore the question of independence between the observed optimization of η 3 and η 4 . This is an important point to address since the relatively strong coupling between BChls 3 and 4 raises the possibility that the optimization observed at one target site could simply be a side effect of optimization at the other target site. From the differences in energy transport mechanism and the somewhat different site distribution of λ ga m,3 and λ ga m,4 (Table II), it can be inferred that these two quantities are reasonably uncorrelated. This can be seen explicitly from Figure 4 where we have recast the data, showing the efficiencies are only weakly correlated.

VI. CONCLUDING REMARKS
To summarize, we have proposed a plausible design principle for FMO which involves an effective interplay between the system and system-bath interaction. Energy transport mechanisms and pathways are combined with the fine-tuning of pigment-protein interaction in such a way that the (inevitable) site-varying system-bath interaction strengths are efficiently exploited to optimize energy transport. Interestingly, the optimization is observed not only for the commonly assigned target site of BChl 3, but for BChl 4 as well. More importantly, this optimization for the two target sites is largely uncorrelated, meaning that the observed optimization at BChl 4 is non-trivial.
In fact, it appears that the electronic system itself may be designed to transport energy to the two target sites via two different mechanisms of energy transport. For BChl 3, the optimal mechanism is purely dissipative funnelling while for BChl 4, the optimal mechanism is partially coherent. This implies that the system and system-bath interaction are configured in such a way that both dissipative and coherent processes are reasonably accommodated. Indeed, we observe an optimization of exciton delocalization when the FMO Hamiltonian is paired with the realistic configuration of site-dependent reorganization energies.
This simultaneous optimization of two largely uncorrelated processes of somewhat opposing nature is noteworthy since it suggests the presence of an underlying evolutionary design principle. It requires spatial "engineering" of system-bath interaction such that the magnitude of the reorganization energies either mutually benefits energy transport to both target sites or at least does not negatively impact energy transport to any of the target sites too significantly. Hence, if the simultaneous optimization at the two target sites is not purely accidental, and if coherence does indeed play a role in the FMO complex, then the conventional view of BChl 3 as the sole linker site to the reaction centre is possibly incomplete.
Finally, we would also like to comment briefly on the 180 cm −1 vibronic mode which has been shown to enhance electronic coherence in FMO 34 . We note that this mode is also resonant with the energy difference between BChls 3 and 4. It is possible that this mode may be relevant to the design principle we have presented in this paper. We are exploring this possibility currently. The exciton population dynamics is described by the following rate equation: where K M RT is the Modified Redfield rate matrix and P (t) = [P 1 (t) P 2 (t) . . . P k (t) . . . P M (t)] is a vector of exciton populations at time t, with the total number of exciton states M = N . The matrix element K M RT k,k of K M RT is the population transfer rate from the k th to the kth exciton state and is given by the following time integral: − 2iλ kk ,k k ])e 2(g kk,k k (τ )+iλ kk,k k τ ) .
Here, λ δσ,κµ and g δσ,κµ (t) are the exciton counterparts of the site reorganization energy λ m and site lineshape function g m (t) respectively with the relationships expressed as: where a k,k = c * m,k c m,k . λ k = λ kk,kk is the kth exciton reorganization energy, which when subtracted from the kth exciton transition energy E k (i.e. the kth eigenenergy of H s ), gives the corresponding 0-0 transition energy, i.e. E 0 k = E k − λ k . Meanwhile, the site lineshape function g m (t) is defined by the following double time integral: where C m (t) is the complex bath correlation function which can be written as: Thus, g m (t) can be computed from the spectral density J m (ω) from the following relation: The integrand in eq. A7 contains poles, and thus the lower limit of integration must be changed to an appropriately small finite value for numerical integration. Alternatively, the poles can be removed at the level of C(t) by analytical means, such as in the case of the Drude-Lorentz representation. The derivation of the analytical form of g m (t) for the Drude-Lorentz spectral density is shown in Appendix D.
Lastly, the diagonal elements of K M RT are computed from: Eq. A1 can then be easily solved with P (t) = e K M RT t P (0).
Here the initial population vector consists of the exciton populations corresponding to the initial excitation site m 0 , i.e. P (0) = [c 2 m0,1 c 2 m0,2 . . . c 2 m0,k . . . c 2 m0,M ] . It should be pointed out that the numerical integration in eq. A2 requires a finite cutoff value for the upper integration limit. In general, the stronger the systembath coupling, the more oscillatory the integrand and the longer it takes for the integrand function to taper off. We have taken care to ensure our chosen cutoff value of τ = 1 ps is acceptable by checking that increasing the cutoff has negligible effect on the calculated rates. For the genetic algorithm and statistical evaluation performed in this work, such manual inspection would be impractical. As such, we have applied a "worst case scenario" test by assigning λ m with the largest reorganization energy to each site and confirmed there is no appreciable change in the computed rates when the value of the cutoff is increased to 5 ps −1 .

Appendix B: Coherent Modified Redfield Theory
The formula for the exciton transfer rates for CMRT is similar to eq. A2, except with the upper limit of integration replaced by t to incorporate non-Markovianity. That is, the rates are time-dependent in contrast to the Markovian MRT case. We shall denote this exciton population transfer rate (from exciton k to k) as K CM RT k,k (t). Due to the inclusion of coherences in CMRT, dephasing rates must also be accounted for and is defined by the following relation: It turns out the CMRT master equation can be cast into a generalized Lindblad form: This allows the master equation to be conveniently solved using the Non-Markovian Quantum Jump (NMQJ) approach. The summary of the NMQJ technique is provided in Appendix C while full details and derivation can be found in Refs. 19, 20, and 35. In eq. B2, A kk = |k k | are the jump operators while H e = k E 0 k |k k| is the modified system Hamiltonian that governs the coherent evolution. The jump rates R kk (t) from exciton state k to k are such that: where G k (t) is the kth element of vector G(t) which is linked to the pure dephasing rate L kk (t). This relation is given by G(t) = B −1 D(t) where B is a matrix and D(t) is a time-dependent vector whose elements are defined respectively as: for k < k < M, 0.5 otherwise, The Non-Markovian Quantum Jump (NMQJ) is a tool for the stochastic unravelling of a non-Markovian quantum master equation, and is the non-Markovian generalization of the well-known Monte Carlo Wave Function (MCWF) method 36,37 . It takes advantage of the general definition of the density matrix, where N is the ensemble size and N α is the number of ensemble members in the state |ψ α (t) , and operates on the level of the state vector. In this work, the |ψ α (t) s are the exciton states |k and the time-evolved initial state |ψ 0 (t) . |ψ 0 (t) is a coherent superposition of exciton states in which |ψ 0 (0) = |m 0 , where |m 0 is the initial localized site basis. Due to the construct of the density matrix, the single exciton states are time-independent as global phases are cancelled out. Hence eq. C1 can be rewritten as: Note that ρ(t) is in the site basis due to eq. 3. According to the NMQJ formulation, each ensemble member undergoes continuous time evolution interrupted by discontinuous probabilistic jumps. The propagation of the state vector proceeds in small time steps δt with the deterministic evolution described by where the effective non-Hermitian Hamiltonian is defined asĤ When R kk (t) ≥ 0, an instantaneous positive jump to another state |ψ α (t) may occur, and the state at the next time step would be set to this new state.: A kk |ψ α (t) A kk |ψ α (t) = |ψ α (t + δt) . (C5) Here the probability of the jump occurring through the channel k → k for a given ensemble member state |ψ α is P + α,kk (t) = R kk (t)δt ψ α (t)| A † kk A kk |ψ α (t) . (C6) Unlike purely Markovian dynamics, the non-Markovian transition rate R kk (t) is time-dependent and can become negative. During this time period, the action of the positive jump operator A kk is to bring the target state |ψ α (t) to the source state |ψ α (t) : This implies that the negative jump operator A − kk = |ψ α (t) ψ α (t)|. In other words, a negative jump means the reversal of a previous jump back to a prior state. The probability for a reverse jump is given by: With each jump, the number of ensemble members in the source and target states are updated accordingly for the current time step, i.e. N α (t) − 1 and N α (t) + 1 respectively. The choice between deterministic evolution or jump is determined by a random number 0 < < 1. If is less than or equals to the total jump probabilities of all channels, a jump occurs and vice versa. If a jump is determined, another random number s is generated to randomly select the jump channel.
At t = 0, all the ensemble members are in |ψ 0 (0) , i.e. N 0 (0) = N . This means the propagation of the density matrix starts from a pure state and progresses to a mixed state through positive jumps. At a later time when the rates become negative, negative jumps can undo these positive jumps, and this may include the revival of coherences via a reverse jump to |ψ 0 (t) . For FMO, there are 56 relaxation and 8 dephasing channels. As can be seen from eq. C2, efficient computation and averaging of the constituent density matrices can be achieved by simply updating N α (t) and N α (t) at each time step and performing a one-off time evolution of |ψ 0 (t) .

Appendix D: Analytical form of the Drude-Lorentz lineshape function
Since the Drude-Lorentz spectral density J DL m (w) has simple poles at ω = ±iΩ m , the residue theorem can be conveniently applied to obtain the analytical form of the bath correlation function C m (t) 38 which is: Here ν k = 2πk β are the Matsubara frequencies. In principle, the Matsubara expansion is infinite but in practice, the summation can be truncated at some finite value n r . The number of terms required for convergence is dependent only on temperature, with more terms needed for lower temperatures. Substituting into eq. A5, we arrive at the analytical form of the lineshape function: