## Abstract

Photosynthetic pigment-protein complexes (PPCs) are a vital component of the light-harvesting machinery of all plants and photosynthesizing bacteria, enabling efficient transport of the energy of absorbed light towards the reaction centre, where chemical energy storage is initiated. PPCs comprise a set of chromophore molecules, typically bacteriochlorophyll species, held in a well-defined arrangement by a protein scaffold; this relatively rigid distribution leads to a viewpoint in which the chromophore subsystem is treated as a network, where chromophores represent vertices and inter-chromophore electronic couplings represent edges. This graph-based view can then be used as a framework within which to interrogate the role of structural and electronic organization in PPCs. Here, we use this network-based viewpoint to compare excitation energy transfer (EET) dynamics in the light-harvesting complex II (LHC-II) system commonly found in higher plants and the Fenna-Matthews-Olson (FMO) complex found in green sulfur bacteria. The results of our simple network-based investigations clearly demonstrate the role of network connectivity and multiple EET pathways on the efficient and robust EET dynamics in these PPCs, and highlight a role for such considerations in the development of new artificial light-harvesting systems.

## 1. Introduction

Photosynthesis is the most important process on our planet, converting sunlight into the chemical energy which powers nearly all biological systems [1–4]. Estimates suggest that, every year, around 10^{14} kg of organic matter are produced globally [5]; this biomass is roughly 20×10^{3} times larger than the estimated mass of the Great Pyramid of Giza in Egypt. The majority of this organic matter is associated with photosynthesis in common green plants, algae and phytoplankton, as well as in more exotic photosynthetic bacteria. Of course, besides providing the energetic basis for production of organic compounds such as sugars and carbohydrates, photosynthesis also results in evolution of oxygen; in short, without photosynthesis, life as we know it would not exist.

All photosynthetic systems operate along the same broad principles (figure 1). Photosynthetic supercomplexes, such as photosystem II found in green plants and algae [8–11] are arranged such that a large collection of light-harvesting antenna complexes, including light-harvesting complex II (LHC-II [12,13]), surround a collection of core antenna complexes, such as CP43 [14]; these core complexes are themselves situated in proximity to the reaction centre (RC). In the first step of the overall photosynthetic process, commonly referred to as light harvesting, the peripheral antenna complexes absorb sunlight, generating an electronically excited state; in a series of electronic energy transfer steps, this electronic excitation energy is subsequently transported to the RC, where it is used to drive the initial charge separation reactions of chemical energy storage. The light-harvesting process and the excitation energy transfer (EET) from the light-harvesting antenna complexes to the RC are facilitated by pigment-protein complexes (PPCs); in green plants and algae, these include LHC-II, CP26 and CP29 [11,13–17], whereas in green sulfur bacteria such as *Chlorobium Tepidum*, the Fenna-Matthews-Olson (FMO) complex is the common PPC structure [18–25]. The known crystal structures of PPCs, including both FMO and LHC-II considered later in this article, show that such complexes are built from a network of chromophore molecules, such as bacteriochlorophyll a (Bchl a), chlorophyll (Chl) or carotenoids, embedded in a protein scaffold which enforces spatial organization on the chromophore system (figure 1). The absorption spectra of these PPCs typically span a broad spectral range from 400 to 900 nm [26], although individual PPCs are often much more finely tuned to absorb sunlight in a smaller window, presumably in a response to evolutionary and environmental pressures and availability of solar radiation. Most importantly, photosynthetic PPCs demonstrate outstanding light-to-charge conversion efficiencies, typically close to 100% [11], while simultaneously allowing EET to proceed rapidly enough (typically tens to hundreds of picoseconds) to avoid competitive excited-state quenching paths [3].

While an enormous amount of research effort has been directed at understanding the core electron transfer chemistry which is induced upon arrival of excitation energy at the RC [3,4,9–11,27–29], our particular interest lies in the light-harvesting process itself, with a focus on the quantum-mechanical dynamics of the energy transportation from peripheral PPCs to the RC. Specifically, *biological* PPCs exhibit many of the properties which would be desirable in *artificial* light-harvesting systems for solar energy generation [11,1]; for example, biological PPCs exhibit highly efficient light-to-charge generation and demonstrate highly directional energy transport towards the RC. Furthermore, biological PPCs often exhibit resistance to photodamage from high-intensity light levels *via* quenching relaxation pathways (typically facilitated and associated with carotenoid chromophores [16]), while the absorption spectra of different PPCs demonstrate that there is a reasonable degree of spectral tunability which is accessible within the chemical domain of such systems. Finally, from a commercial and economical viewpoint, it is notable that biological PPCs achieve these amazing feats of energy transportation using predominantly organic building blocks, without (toxic and expensive) heavy metals.

It is, therefore, not surprising that much computational and experimental research over the last decade or so has focused on developing an atomic-level, or even electronic-level, view of how biological PPCs facilitate efficient and directional EET [3,9,10,16,17,20,30–56]. Much of this work has focused on understanding the role of quantum mechanics in guiding EET in PPCs, despite the initial perception that such effects might not be important [26]. In particular, biological PPCs operate in a regime where thermal fluctuations, as well as strong interaction with the wider solvated environment, would be expected to quench any quantum mechanical coherence effects in the EET dynamics; however, both experimental and computational studies suggest that this is not the case. A seminal set of experiments by Fleming and co-workers [3,20,40,48,52,57] has demonstrated that coherent ‘beats’ can be observed in the two-dimensional electronic spectra of FMO, corresponding to the signature of electronic coherence during the first few hundred femtoseconds of EET dynamics following photoexcitation [40]. While these measurements were initially performed at cryogenic temperatures, later experiments demonstrated that quantum coherence is also observed under physiological temperatures (albeit significantly damped relative to low temperature studies) [58,52]. The role of quantum coherence has been supported by a wide range of computational studies, performed using a variety of different methods ranging from quantum density matrix propagation [29,32,34–36,42,47,49–51,56,59,60] to atomistic simulation methodologies [30,61,62]. Here, quantum dynamics simulations have predicted the existence of quantum coherent oscillations in the electronic populations in many-chromophore systems such as FMO, while analysis in terms of quantum information theory [9,42,44,48] have enabled quantification of the essential underlying features of these quantum effects. This synergy between computational and experimental predictions remains a ‘hot topic’, and has even engendered more philosophical questions arising regarding the role of fundamental quantum-mechanical features within a chemical process which is so important to life [11,20,26].

Our own recent work [56], further expanded here, has focused on using quantum dynamics simulations to model EET in PPCs; however, rather than the nature of quantum mechanical fluctuations themselves being of interest, our research has instead focused on the role of structural organization in PPCs. In particular, we note that, despite the underlying EET and chemical processes of photosynthesis being common to a wide range of biological systems, there is great diversity in the structural motifs adopted by light-harvesting complexes of different organisms. This is in surprising contrast to the RC, which is known to exhibit strong similarities across plants, algae and bacteria [63]. For example, the EET network of FMO is a seemingly unstructured arrangement of eight Bchl molecules [18–25], while the light-harvesting complex found in purple bacteria comprises a highly symmetric ninefold ring structure [64]. Biological PPCs also span a reasonably broad range of complexity. The FMO complex of green sulfur bacteria contains just eight chromophore molecules, while spinach LHC-II is a trimer of proteins comprising 42 chlorophylls [12]. Finally, it is interesting to note that biological PPCs also incorporate a varied range of auxiliary chromophores beyond the primarily light-harvesting pigment molecules (typically chlorophyll). For example, carotenoids are particularly prevalent co-binding chromophores in PPC systems. Given the homology in photosynthetic RCs across species separated by millennia of evolution, there are clearly other factors at play in driving the selection of light-harvesting antennae complexes and EET networks in biological systems [63,65,66].

So, while much of the contemporary experimental and theoretical interest in PPCs has focused on the role of quantum-mechanical effects in the mechanism of EET, the aim of this research is to address a more fundamental question: why are PPCs built the way they are? Of course, addressing such fundamental questions is not possible using approximate quantum simulations of highly simplified models of PPCs. For example, given evolutionary time scales, it is tempting to ask whether PPCs are in any way optimized for their role in photosynthesis [34,36,45–47,49,56,67] or whether external factors, such as ease of biological synthesis [21], are instead responsible for the observed molecular organization in PPCs. Such questions are impossible to address from the standpoint of molecular simulations alone, but addressing these issues and better understanding the rationale behind natural PPC structural motifs could, in turn, lead to a better understanding of how to construct improved artificial photosynthetic systems for water-splitting or more efficient light-harvesting devices for solar energy applications [1,68]. Understanding the ‘design rules’ which underpin molecular organization in biological PPCs is one of the goals of our current research, and our efforts in this direction are highlighted in this article.

### (a) Photosynthetic pigment-protein complexes as networks

The focus of recent work from our group, as well as this article, has been on how one can use simple ideas emerging from the theory of networks (or graphs), along with quantum dynamics simulations of EET, in order to interrogate and, ultimately, rationalize the structural motifs adopted by biological photosynthetic PPCs. The connection between PPCs and networks arises due to the fact that the chromophore molecules in PPCs are held in a reasonably rigid relative spatial orientation, with the protein scaffold serving to enforce this arrangement while also, as noted below, modifying the local environment of each chromophore. This is emphasized in figure 2, which shows crystal structures determined for four different PPCs, highlighted to emphasize the spatial arrangement of the chromophore molecules. As a result, a simple coarse-grained view is to imagine a PPC as a collection of chromophores arranged in space in fixed relative positions and orientations, with the protein simply viewed as an effective environment. Thus, the chromophores may be viewed as forming the vertices of a network. As elaborated later in terms of the Frenkel exciton Hamiltonian, the edges of the network can be viewed as electronic couplings between the different chromophores. Typically, all chromophores interact with each other to some extent (albeit negligibly weakly in some cases), such that the chromophores form a so-called ‘fully connected’ network.

By reducing a complex PPC system, possibly containing many thousands of atoms, to a few localized excited states and the electronic couplings between them, this network-based view offers a great deal of scope in addressing some of the fundamental questions and features of biological PPCs. For example, is there any inherent advantage to symmetric chromophore network organization compared with random chromophore network arrangements, as found in FMO [46]? Are there particular features of the electronic couplings, spatial arrangements or orientations of chromophores in biological PPCs which are common across different PPC structures? Which features give rise to directional EET in many-chromophore systems, and are these features common across different PPC structural motifs [33,34,36,42,44,47,50,55]? Finally, and perhaps most importantly for this work, what can be learned about chromophore arrangements in PPCs to enable the design of better artificial light-harvesting and EET systems relevant to solar energy generation, photocatalysis and sensing? All of these questions can be approached by adopting a coarse-grained, network-based view of EET in PPCs.

Of course, we must note that this network-based viewpoint has featured prominently in PPC research over the last decade or so. Most importantly, Silbey and co-workers have used the idea of PPCs as networks to explore the role of a number of important physical features in EET dynamics [34,36,47]. For example, by developing network-based models of EET in multi-chromophore systems which can be solved analytically, these models have demonstrated the subtle interplay between environmental noise and chromophore coupling strengths, demonstrating that typical parameters observed in biological PPCs lie well within the regime in which EET is predicted to be highly efficient. Furthermore, the role of symmetry in PPC trimers was also considered [46]; here, a network-based viewpoint helped to demonstrate that symmetric ring-structures of the order of those found in nature are predicted to enable optimized and minimally frustrated EET dynamics and optimal packing density. Finally, we highlight interesting results which focused on mapping a quantum-mechanical density matrix-based model of a PPC onto a classical kinetic master equation [36]. By comparing the results of quantum and effective-classical models of the same PPC, this research demonstrated that the role of quantum coherence was not particularly significant with regards to EET efficiency or trapping times. However, it should be noted that this quantum-to-classical mapping was performed for a simple quantum model which modelled the environment as a Markovian ‘white noise’ which gives rise to pure decoherence (in fact, this is the same model we employ below); while this model appears to give results which are physically consistent with more complex simulation methods, the role of spatio-temporal correlation in the PPC and environment remains unclear. As a final point, it is worth noting that the idea of PPCs-as-networks has also been extended to other non-biological systems which can similarly be viewed as a collection of coupled chromophores. For example, organic photovoltaic materials have been modelled using system Hamiltonians similar to the Frenkel exciton model studied below, accounting for the role of intermolecular electronic coupling and environmental fluctuations [71].

By building on the idea of PPCs-as-networks, the work presented here extends on our recent investigation of biological PPCs by comparing and contrasting the network properties of two different PPCs; the FMO complex considered in our previous work, and LHC-II. Here, the aim is to assess to what extent we begin to observe common network features associated with PPCs exhibiting efficient EET. As a result, our simulation approach is similar to that adopted in our previous research [56], focusing on combining quantum dynamics simulations using a simple yet physically motivated methodology with a network-based analysis of calculated EET efficiencies. By comparing results for FMO monomers and with LHC-II (monomers and trimers), and using simple insights from network analysis calculations, we find that a common picture emerges, with PPCs exhibiting robust and efficient EET predominantly as a result of the high degree of connectivity in the chromophore network.

The remainder of this article is organized as follows. In §2, we highlight recent developments in computational methods for simulating EET in PPCs, before highlighting the methodology adopted in our work; similarly, we outline the main network-based tools employed to understand the role of inter-chromophore connectivity and environmental noise on EET dynamics. In §3, we present a detailed comparison of the EET properties of FMO and LHC-II, focusing on correlating the observed EET efficiency trends with network features observed in biological PPC monomers and trimers; furthermore, we seek to determine whether common network features can be identified which can be associated with efficient EET. Finally, in §4, we conclude and highlight some avenues for further work.

## 2. Theory

In this section, we outline some of the main computational approaches taken to date in modelling EET dynamics in biological PPCs. First, reflecting the focus of our recent work, we outline model Hamiltonian approaches for modelling inter-chromophore interaction and chromophore–environment coupling. Second, we describe common methods for modelling EET dynamics, again reflecting the density-matrix-based methods adopted in our recent work. Third, we highlight several network-based analysis tools which can be used to interrogate the properties of networks, paying attention to how these methods can be modified to interrogate biological PPCs. Finally, we outline the computational approach taken in the new results presented herein.

### (a) Hamiltonian description of biological pigment-protein complexes

Biological PPCs pose many of the same challenges to molecular simulations as other biomolecules such as proteins and enzymes; they are typically large (usually several thousand atoms in the simplest cases), require explicit treatment of the condensed-phase solvent environment, and exhibit long time-scale molecular motions. In addition, by the nature of the EET process, one is also faced with the challenge of calculating excited-state molecular properties, namely excitation energies and inter-chromophore couplings, in biological PPCs. As a result, so-called ‘on the fly’ simulations of EET dynamics in biological PPCs, coupling accurate *ab initio* treatment of electronic structure with simultaneous evolution of nuclear dynamics, are an enormous computational challenge which is only just beginning to be explored [61].

Instead, the majority of simulations of the dynamics of EET adopt the approach of mapping the complex many-atom biological PPC, and the associated chromophore electronic states, onto a simpler system-bath-type quantum Hamiltonian which is more amenable to simulations [3,20,29,34–37,42,47,51,60,67]. The most common of such models is the Frenkel exciton (FE) Hamiltonian, describing the electronic sub-system comprising chromophore electronic states, the interchromophore couplings, and an environment (typically taken to be a harmonic ‘bath’); the interaction between the electronic and environmental systems is also included in this model.

The FE Hamiltonian can be written as
*ϵ*_{i} is the electronic excitation energy of site *i* and *J*_{ij} is the electronic coupling element between sites *i* and *j*. Formally, the excitation energies are determined by solving the electronic Schrödinger equation for each chromophore molecule in the protein, but in the absence of the other chromophore molecules [10,31,53,54,61]. Both the excitation energies and coupling elements can be calculated using *ab initio* electronic structure methods [10,53,61,72] although the large scale of such calculations for PPCs demands high-performance computing strategies, such as linear-scaling methods [72] or GPU-enabled approaches [61]. Instead, a more common route to determining the parameters of the electronic sub-system Hamiltonian is to use an approach such as the dipole-dipole approximation [54,35], derived by assuming that the transition densities can be treated within a point approximation. Furthermore, we note that while the excitation energies *ϵ*_{i} can be calculated by *ab initio* methods, an alternative approach is to use the known experimental fluorescence spectrum of the PPC complex [31,73]. Whichever method is used to extract the electronic parameters, there is, of course, an inevitable uncertainty attached to each parameter which must be borne in mind throughout any further analysis. Our simulations below adopt the common procedure of using parameters derived from both experimental excitation spectra and the transition-dipole-cubed method. In this manner, we can straightforwardly compare our simulation results to previous calculations employing the same parameters within different quantum simulation protocols [31,35,56,73,74].

At this point, the connection between the viewpoint of the PPC as a network and Hamiltonian dynamics should now be apparent. The electronic Hamiltonian, *N*-chromophore PPC system is an *N*×*N* matrix when evaluated in the basis of chromophore excited states, with excitation energies *ϵ*_{i} occupying the diagonal elements and coupling elements *J*_{ij} occupying off-diagonal elements. The magnitude of these off-diagonals can be viewed as representing connection *weights*. Because each chromophore is coupled to every other chromophore (albeit with potentially weak interaction strengths), the electronic sub-system can be viewed as a fully connected network comprising chromophores as network vertices, with the inter-chromophore couplings defining the network edges. This picture of the electronic sub-system of the PPC acting as a connected network will be employed later in analysis of EET dynamics in FMO and LHC-II.

The bath (environment) Hamiltonian is most commonly assumed to represent a set of harmonic oscillators with a frequency spectrum chosen to model the PPC environment at the relevant physical conditions. In other words, we have
*ω*_{k} is the vibrational frequency of normal mode *m*_{k} is usually, by definition, equal to one. The vibrational spectrum of the bath, like the electronic parameters in

The final element of equation (2.1) is the system-bath coupling, *electronic* sub-system and the *vibrational* bath; as noted later this coupling is essential in the phenomenon of environmental noise-activated quantum transport (ENAQT) which has been developed by several groups to interpret the interaction of environment and EET [37,38,42,44,50,60]. Most commonly, a linear coupling is adopted in *i*, or site-dependent; furthermore, the above equation can also be modified such that a different set of bath modes is available to each chromophore site. While the parameters of the chromophore sub-system and vibrational bath can be determined through a combination of *ab initio* calculations and classical molecular dynamics, the coupling parameters

As a final point, we note that two additional contributions to the Hamiltonian of equation (2.1) are usually incorporated into the full model before simulating EET dynamics [37]. First, in the biological PPC, the electronic excitation energy ultimately traverses the PPC structure and passes into the reaction centre (RC), where charge separation is initiated; from the point-of-view of simulating the EET dynamics, this requires that population density exits the region of space defined by the PPC Hamiltonian of equation (2.1). Technically, this feature of the EET dynamics is captured by adding a non-Hermitian trapping operator at chromophore ‘sink’ site *m*, typically of the form
*κ*_{t} is the trapping rate. The upshot of this term is that population density exits the PPC at a rate *κ*. Experimental estimates suggest that the value of *κ* is typically of the order of 1 ps^{−1} [37,39,76]. In addition to the trapping term, a further non-Hermitian term is added to the Hamiltonian of equation (2.1) which is of the form,
*Γ*, is typically chosen to be 1 ns^{−1}, again following previous experimental and theoretical estimates [37,39,76].

### (b) Simulating excitation energy transfer dynamics

Although the Hamiltonian of equation (2.1), as well as the extensions to account for recombination and trapping, represents a physically motivated model with well-defined parameters, full quantum treatment of the dynamics of such a Hamiltonian is usually beyond the system-size treatable by accurate quantum dynamics simulation methods. While system-bath methods have been developed which enable one to model harmonic environments coupled to electronic sub-systems, such as reduced density matrix methods [77] and influence functional approaches [29,78], these are generally too computationally demanding to be applied to electronic sub-systems containing many electronic states (i.e. PPCs) and converged bath sizes. However, it is worth noting that the hierarchical equations-of-motion (HEOM) method has been developed specifically to accurately model PPC-type systems with enormous success [48,52]. For a harmonic bath, HEOM simulations can be viewed as highly accurate solutions of the true EET dynamics of PPCs, providing an important (yet, again, computationally demanding) benchmark method. A further method which has been employed to study PPCs is the nonlinear network model [79] which has been shown to be a less computationally expensive method for including the microscopic effects than many of the above methods while reproducing experimental data [80]. In this way both the FMO complex and its protein environment may be treated under a network-based model.

The density matrix formulation of quantum dynamics is particularly suited to modelling EET phenomena [81]; specifically, the reduced density matrix formalism allows one to focus attention on determining system properties of interest only, rather than focusing on the dynamics of the full system. Furthermore, the connection to thermal time-dependent observables is more straightforward within the density matrix picture, while the properties of the density matrix themselves are readily interpreted as site populations (diagonal elements) or site–site coherences (off-diagonals). This density matrix picture has been adopted in a wide range of simulation approaches, all of which attempt to model the essential EET dynamics of the electronic sub-system while also accounting for the role of environmental fluctuations to differing extents. For example, the well-known Bloch-Redfield method, derived using a short-time expansion of the density matrix propagator in combination with a Markovian environment assumption, leads to a computationally tractable approach. However, the Markovian assumption implicitly neglects time correlation in the environmental response of the system, while the weak-coupling assumption (embodied in the assumed separability of the density matrices for the system and environment) can lead to underdamped dynamics in simulations of EET in PPCs, where coupling between system and bath is often of comparable magnitude as that between different chromophore sites [48]. Similar ideas have led to the development of a wide range of different EET simulation methods, including the reduced density matrix hybrid method [82,77], the partially linearized density matrix method [83] and the HEOM approach noted above [52].

Given the enormous challenge associated with accurate propagation of the density matrix, our own work, reported previously and expanded in this article, has adopted a more pragmatic approach to modelling EET dynamics which moves away from explicit treatment of the environment in favour of a computationally simple approach, for reasons explained below. Specifically, we employ the Haken-Strobl (HS) model [32,37,59] to describe EET dynamics of the electronic sub-system, while also accounting for the environment as a Markovian ‘noise’ term. By assuming that the site-energies of each chromophore experience Markovian fluctuations, it is possible to derive an equation-of-motion for the *electronic subsystem* density matrix as
*γ*_{i} is the dephasing rate at site *i*; in essence, the environmental dephasing factor acts in a similar way to the trapping and recombination factors described above, although the dephasing effect is due to interaction with the environment rather than escape *via* the ‘sink’ site, for example. A typical value for *γ*_{i} is around 100 cm^{−1} although, as shown in previous work, the phenomenon of ENAQT means that robust EET is observed for a wide range of dephasing rates, from around 40 cm^{−1} to 3000 cm^{−1} in FMO. Throughout this article, we assume that each chromophore site experiences the same global dephasing rate; the influence of local variations in *γ*_{i} have been discussed in our previous work [56], which concluded that the balance between exciton transport and local exciton trapping is key to efficient EET, as similarly noted in other previous (analytical and numerical) studies [34,36,47].

Given the relative simplicity of the HS model in investigating EET dynamics, particularly when compared with more sophisticated treatments of the environment, it is worth commenting on its application here. As we have shown previously, the HS model qualitatively reproduces the key features of the EET dynamics in FMO when compared with methods such as HEOM [52]. Generalised Bloch-Redfield ((GBR) [35]) and Linearized Semiclassical Initial Value Representation ((LSC-IVR) [51,84]). Furthermore, as in our previous work, we emphasize that our main interest is in the network features of the electronic sub-system; we expect that the impact of changes to the electronic Hamiltonian ^{3} independent simulations for each of the systems considered here); the computational efficiency and simplicity of the HS model is essential in performing this large number of simulations.

Of course, it is always worth bearing in mind the approximations being made in any given simulation approach, and this work is no different. The two key features of the HS treatment of the environment are (i) the absence of explicit vibronic coupling as a result of the uncorrelated ‘white noise’ nature of interaction between electronic subsystem and environment, (ii) the use of a single dephasing rate parameter to describe environmentally induced decoherence. By comparing HS simulations to previous simulation approaches using explicit vibronic coupling (i.e. between a parametrized harmonic bath and the electronic subsystem, as in the HEOM, GBR and LSC-IVR methods), our results indicate that the HS model is capable of calculating the *qualitatively* correct dynamics in FMO [56]. This gives us confidence in our overall simulation approach, but means that we cannot address questions of interest regarding explicit environmental coupling. The introduction of a single white-noise dephasing rate could also, in principle, be addressed, and the influence of site-dependent dephasing rates within the HS has already been addressed in our previous work [56]. Furthermore, work by Coker and co-workers, already noted above, has shown that site-specific environments can appreciably change the dynamic characteristics, suggesting that simulations performed using the HS model are unlikely *quantitatively* correct. That said, and as emphasized previously, the HS model provides a qualitative view of EET dynamics in generic quantum models such as that used to represent FMO here. Our future work, expanding the study presented below, will focus on using molecular dynamics simulations as a route to introducing explicit vibronic coupling. The influence of network robustness remains to be seen, but is certainly an interesting topic.

### (c) Network analysis tools

The viewpoint adopted herein, where the quantum Hamiltonian is viewed as a network, is advantageous for several reasons. First, it enables an enormous reduction in system complexity, where we move from an atomistic model of the PPC, as provided by the crystal structure, to a much simpler system defined by a small number of vertices and well-defined connections. Second, the effect on EET dynamics of modifications to the quantum network can be readily visualized and directly understood in terms of modifications to inter-chromophore coupling elements or site energies. Finally, by mapping a quantum Hamiltonian of a complex molecular system onto a much simpler quantum network, it might be hoped that the results and trends obtained are more general than might be obtained by considering the explicit EET dynamics of, for example, an atomistic model of FMO. For example, the appearance of multiple EET pathways and their impact on robustness of the quantum transport network is easily observed using a simplified network picture, and allows one to draw conclusions which might be readily expanded and tested in other related quantum networks.

One approach to analysing the properties of networks, and the main approach we employ here, is to ‘take apart’ the underlying network structure of the system and measure the response of observables to any changes [56,85,86]. In this regard, the main property which we are interested in for EET in PPCs is the EET efficiency; following previous work, we define the EET efficiency as the time-integral of the population density at the selected ‘sink’ site relevant for the PPC structure of interest. Thus, we have
*m*, the site population *p*_{m}(*t*)=*p*_{mm}(*t*), and the integral is calculated up to some maximum time which is sufficiently long to capture the dominant population dynamics at the chromophore sites; throughout this work, we employ *t*_{max}=10 ps, which is much longer than the natural time scales for EET dynamics in PPCs (∼1/*κ*_{t}=1 ps) [39]. The efficiency *η* gives a simple indicator of the rate of population transfer throughout the chromophore network; given that the ultimate role of a PPC is to enable EET, *η* is a good measure of the response of the network to perturbations.

Our approach to investigating the role of network connectivity is relatively straightforward, although quite time-consuming (hence, the use of the simple HS approach in modelling EET). In short, we perform a large series of ‘knock-out’ studies in which a given chromophore site or inter-chromophore coupling element is removed from the system and the EET efficiency is subsequently evaluated using the HS model for the remaining ‘disrupted’ PPC. This process is repeated for all possible choices of chromophore removals and coupling removals, calculating EET efficiency each time. Furthermore, this process of removal and evaluation is then repeated for combinations of chromophore pairs, chromophore triplets, and so on; as a result, our approach requires a combinatorial number of EET efficiency evaluations, thereby justifying our previously mentioned requirement of a fast EET dynamics method. We note that, in the case of the smaller PPCs considered below, the number of possible combinations of chromophore removals and coupling removals are often such that we can evaluate EET efficiency for all permutations; however, if the number of possible permutations is deemed to be too large, we resort to sample a (still large) number of realizations in order to achieve adequate sampling. In such cases, the sampling is entirely random, with each realization independent from each other. The outcome of this ‘knock-out’ analysis is that we can interrogate the impact of disrupting the chromophore network at the level of the population dynamics, providing direct insight into how the electronic properties of each PPC system influence efficient EET.

As a further analytical approach, we note that the idea of using the Hamiltonian

### (d) Simulations details for this work

Our results below consider EET dynamics, and the influence of network organization, for two different PPCs, although the EET dynamics simulations are performed in the same way in each case. In particular, we use the HS model to perform propagation of the electronic density matrix, with trapping and exciton recombination terms included in each case. A standard fourth-order Runge–Kutta integration algorithm was used to propagate the density matrix equations-of-motion with a time-step of 0.5 fs. Unless otherwise stated, all simulations used an exciton recombination rate of *Γ*=1 ns^{−1}, consistent with estimated values for the *in vivo* FMO system [37,39,76]; similarly, the exciton trapping rate at the chosen sink site in each system (see below for details) was *κ*_{t}=1 ps^{−1}. Finally, unless otherwise stated, the environmental dephasing rate on all sites was *γ*=100 *cm*^{−1}, again typically sitting in the response regime in which efficient EET is observed in biological PPCs such as those considered here.

## 3. Results and discussion

In general, as discussed above, PPCs consist of an array of chromophores, which are responsible for the efficient transfer of electronic energy to the RC. One (or indeed more) chromophores are first excited by the absorption of a photon; such chromophores are referred to as ‘source’ sites hereafter. Intermediate chromophores then transfer the electronic energy towards an energy sink in the array, namely one chromophore (or perhaps more than one chromophore) which resides towards the bottom of the potential gradient of the chromophore array, where energy flows out of the PPC and into the RC. These ‘exit’ chromophores are referred to as ‘sink’ sites hereafter.

In this article, we investigate EET in two exemplar PPCs, the first being the FMO complex found in green sulfur bacteria (GSB). The FMO complex consists of a trimeric system of Bchls inside a protein scaffold, whose spatial distribution is well conserved across GSB species [89,90], with recent crystal structures reporting the presence of eight Bchls in each monomer. As is well known from previous experimental and computational studies, each monomer within the trimer is identical and behaves essentially independently, thus only the monomer is considered throughout this article [48,37]. The *eighth* Bchl was only recently reported due to the propensity of its loss in the crystallization processes [91], and is considered to be the linker between the FMO complex and the chlorosome baseplate [70,92]. Thus, both the 8-Bchl model and the previously known 7-Bchl model have been well studied in the literature. Here, we restrict our discussion to the 7-Bchl model to draw comparison with previous work on the 8-Bchl model [56]. The spatial distribution and the labelling scheme for the Bchl in FMO are given in figure 3*a*. Given the relative positioning of Bchl 1 and 6, adjacent to the light-harvesting antenna complexes of the chlorosome baseplate, these are considered the source sites. There is great interest in whether these behave as single excitation sites or if energy is delocalized across both sites [92], however in this work they are considered isolated single excitation sites. Multiple pathways exist which funnel this energy towards the sink, taken to be Bchl 3 given its lower excited state energy and its position close to the reaction centre [92–95].

The second PPC considered here is LHC-II, found in higher plants. The specific LHC-II we consider is from the spinach plant (*Spinacia oleracea*), and consists of a trimer of 14 Chls distributed in ring-like structures embedded inside a protein scaffold (figure 3*b*, [96]). Excitation most likely takes place at Chl 1 or Chl 8, based on their positioning near the top of the protein backbone of the complex. Electronic energy then flows towards Chl 4, the sink site, based on the one-photon excitation energies of each Chl, and their positioning relative to the reaction centre [96].

The electronic Hamiltonian for each of the PPCs considered here were the same as in previous studies on the FMO complex [97], and the LHC-II monomer and trimer [96], and are explicitly given in the electronic supporting material (SI). All simulations are carried out using the methodology discussed in §2, unless stated otherwise.

### (a) Population dynamics

The most fundamental simulation one can perform on these PPCs is to simply excite (populate) the source Bchl, and calculate the total population which reaches the sink Bchl at each time step of the simulation. The results are the population dynamics for each Bchl, with representative examples given in figure 4. For FMO, initial excitation at Bchl 1 (figure 4*a*) leads to a rapid decay over the first few 100 fs, after which population oscillates from Bchl 1 and the strongly coupled Bchl 2 neighbour. Within the first 3 ps, the majority of the population has moved away from these two Bchls and has distributed across the entire FMO complex, steadily reaching the sink Bchl 3 site where the population is trapped. A very similar observation is made for initial excitation at Bchl 6 (figure 4*b*). A rapid decay of population occurs within the first few 100 fs, however noticeably, the population on Bchl 3 grows quicker, which can be explained by the multiple strongly coupled pathway to the sink site (figure 3). Importantly, these results agree qualitatively with other methods which treat the protein environment more precisely [95], with oscillatory Bchl populations and an overall population decay towards the sink. As noted previously [56], these simulations serve as a ‘reality check’ to justifying the use of the simple HS model in rest of the simulations discussed in this article.

Moving to the larger, more complicated PPC, the population dynamics of the LHC-II monomer for initial excitation at Chl 1 (figure 4*c*) and Chl 8 (figure 4*d*) display a much slower population decay from the initial excited Chl; it takes several picoseconds for the majority of the population to move away from the initial Chl. Compared with the FMO monomer, this may be understood as arising due to the much smaller average coupling strength between Chls (|*J*_{ij}|∼8 *cm*^{−1} in LHC-II compared with |*J*_{ij}|∼24 *cm*^{−1} in FMO), and the larger network diameter of LHC-II, as discussed in §b. Adding in the other two monomers to form the native LHC-II trimer shows a noticeable change to the network. For excitation at Chl 1 (figure 4*e*), the population decays much faster than the isolated monomer, a feature which can be understood by noting the introduction of an additional strong couplings to Chl 9 of neighbouring monomers (figure 4*e* (purple dashed line) and SI). By contrast, for initial excitation at Chl 8 (figure 4*f*), the additional inter-monomer couplings are generally smaller than the intra-monomer couplings and so the dynamics appear essentially unchanged relative to the LHC-II monomer.

### (b) Robustness to network change

The population dynamics simulations are useful to understand the overall picture of EET in these PPCs, for example enabling visualization of the dominant pathways which funnel energy towards the sink, and evaluation of the time scale of such EET. However, these simulations are not necessarily informative in answering questions about a PPC's robustness to modifications; instead, breaking the network up and measuring the effect on EET relative to the original network gives us information on what each piece contributes to the EET of the overall network.

The question of what happens to EET when a PPC is modified or perturbed in some way remains important because such modifications might result from genetic mutations of the PPC, chemical damage to the PPC such as the loss of a Bchl/Chl, or the changing environment of the PPC, all of which might occur in the photosynthesizing organisms. To model this within a network-based approach is straightforward; Bchls/Chls (i.e. the network vertices) and the couplings between the Bchls/Chls (i.e. the network edges) are removed in a systematic manner, and the EET efficiency as given in equation (2.10) is calculated and compared to the unperturbed network efficiency. This is the rationale for the following investigations. While simulations allow direct access to such structural modifications, experiments are much more challenging. For example, while genetic modification has allowed production of FMO with isotopically substituted atoms (thereby altering vibronic structure) or modified chromophores [98], deletion of pigment molecules is, to the best of our knowledge, difficult to achieve experimentally.

We first investigate the influence of removing Bchl and Chl sites completely from the network (including all couplings to that site). Bchls and Chls are deleted systematically, exploring the full combinatorial space (i.e. all pairs, trimers, quartics, etc. of Bchls) and the EET efficiency is then recalculated using quantum propagation within the HS model. We note that the identified source and sink Bchls in each PPC are protected, and never removed from the electronic Hamiltonian. The results of this ‘chromophore knock-out’ study for FMO and LHC-II monomers are shown in figure 5.

We consider the seven-site FMO monomer first (figure 5*a*). FMO, for initial excitation at Bchl 1 or Bchl 6, displays a surprisingly high degree of robustness to the removal of Bchls, as noted in our previous work on the eight-site FMO model [56]. For instance, when three Bchls are removed, representing 60% of the available Bchls to be removed (excluding source and sink), EET efficiency only drops by ∼25%. This result is in close agreement with similar ‘knock-out’ studies the 8 Bchl FMO model [56], which suggests that the eighth Bchl site does not contribute significantly to additional robustness; we note that this is also consistent with the relatively weak electronic coupling between the eighth Bchl site and the rest of the Bchl sites in the FMO complex. Furthermore, we note that both the simulations of the population dynamics and the chromophore removal simulations display the same trends as discussed here when repeated using a different proposed parametrization of the FMO Hamiltonian, suggesting the results are not specific to our chosen Hamiltonian parameters (see the electronic supplementary material).

Similar observations are made for the LHC-II monomer, shown in figure 5*b*. For initial excitation at Chl 1, half of the available Chls can be removed with little change to the EET efficiency, reflecting the trend observed for FMO. However, for initial excitation at Chl 8, EET efficiency drop almost linearly with the removal of Chls, demonstrating a much lower degree of network robustness for this initial excitation. This difference in behaviour is most likely related to the fact Chl 8 is positioned close to the sink Chl 4, combined with the generally small inter-chromophore couplings to neighbouring Chls, meaning that EET is dominated by only a few efficient pathways following excitation at this Chl; as a result, the EET efficiency is more sensitive to Chl knock-outs. By contrast, for initial excitation at Chl 1, there exists many possible EET pathways to the sink site (see below), each of which is broadly similar in EET efficiency; an increased number of possible transport pathways implicitly increases the robustness to chromophore knock-out because EET can simply be ‘re-routed’ towards the sink.

These knock-out studies can also be performed for the edges of these networks; in other words, we proceed to remove inter-chromophore coupling elements between the Chls. Here, a number of couplings, *N*≤20 for FMO or *N*≤90 for LHC-II, are removed randomly and the EET efficiency is recalculated. To circumvent the combinatorial increase in the number of simulations required to sample *all* edge removals, an ensemble of 1000 repeats of this procedure is averaged to sample the total combinatorial space. For the FMO monomer (figure 6*a*), we see that EET is much less robust to the removal of the couplings (or edges) than to the removal of individual Chl sites (vertices). The removal of 15 couplings, comparable to the loss of three Chl sites, causes EET efficiency to decrease by around 75%; this is in strong contrast to the results of figure 5, which demonstrated that removing three Chl sites decreased the efficiency by around 20%. The same trend is also observed for excitation at Chl 8. Furthermore, this observation is also apparent for the LHC-II monomer (figure 6*b*). Specifically, for excitation at Chl 1, the removal of 63 couplings, equivalent to removing six chromophores, results in a drop in EET efficiency by approximately 60%. The same simulation for initial excitation at Chl 8 suggests a slightly large drop in EET efficiency which can be attributed to the observed decreased robustness as discussed for chromophore knock-outs.

These results for the seven-site FMO PPC and LHC-II are in clear agreement with our previous investigations of the eight-site FMO complex. In particular, we find that EET efficiency in these systems is much more robust to the removal of network vertices (Bchl/Chl sites) than to the *random* removal of network edges (inter-chromophore couplings). As we have noted previously, these results can be viewed as implying that these PPCs exhibit a collection of EET pathways from source to sink sites which have similar coupling strengths along their length; random removal of edges would be expected to have a much stronger impact on such a network because it is possible that many of these strongly coupled EET paths are simultaneously ‘switched off’ by random vertex removal when compared with random node removal. Furthermore, this robustness to node removal is also reminiscent of the features of ‘small world’ networks, such as those found in social media networks; this point is further discussed below.

We can extend these studies of the LHC-II monomer to consider the effect of the additional EET pathways available when other LHC-II monomers are available, as in the LHC-II trimer. Here, we assume that initial excitation still takes place on either Chl 1 or Chl 8 of just one of the monomers; however, *each* LHC-II monomer has a sink site at Chl 4 with identical trapping rates. Thus population is allowed to flow from the initially excited monomer, to the other two monomers, which can trap population at their respective Chl 4 sinks (figure 3). The reported EET efficiency is calculated as the sum of the individual EET efficiencies of each monomer. Here, an increasing number of edges (*N*) are removed from the full LHC-II trimer electronic network, and the EET efficiency is recalculated, as in our previous investigations of FMO and the LHC-II monomer. However, instead of averaging 1000 simulations of randomly removed edges to sample the total combinatorial space, we apply an additional restriction; in particular, the largest coupling is removed first, followed by the two largest and then the three largest, and so on, continuing up to the *N* largest couplings. The results of these simulations are shown in figure 7.

Specific increases and decreases in the EET efficiency are observed for the particular removal of couplings between the Chl of one monomer and the Chl of another, for example when *N*=1,6,25 or 41 highlighted in figure 7. In these simulations, the effect of removing specific strong coupling elements on the EET efficiency can be directly mapped onto the LHC-II trimer structure. In the case of the removal of the largest coupling (i.e. *N* = 1) from the electronic Hamiltonian matrix which couples two of the LHC-II monomers (*N*=41), results in a large drop in EET (blue line in figure 7). This decrease in EET efficiency is traced to the removal of the coupling between Chl 9 of the first monomer and Chl 6 of the second monomer. This pathway provides a route to the sink site of the third monomer, so its removal reduces EET efficiency. Similar observations are seen for other values of *N*, as illustrated in figure 7.

A further network-based analysis of EET dynamics is the evaluation of the network path-length and diameter, which was briefly discussed in terms of the population dynamics in §a. As noted above, inter-node ‘distances’ are defined as the reciprocal of the inter-Bchl coupling matrix elements of the electronic Hamiltonian (equation (2.11)). Following this, one can remove Bchls in a systematic way, just as in the start of this section, and recalculate both the path-length and network diameter to understand the extent to which network connectivity has changed.

Figure 8 shows the calculated network diameters and average path-lengths for FMO and the LHC-II monomer. Considering the average path-length first, figure 8*a*,*b* for FMO and LHC-II monomer, respectively, we see that these measures are reasonably insensitive to the removal of the Bchls/Chls in each network for both excitation sites considered; many sites must be removed before any significant change in the average network path-length is observed. This observation is consistent with our previous work [56] and our results above, indicating that there are several strongly coupled EET pathways in these PPCs; as a result, removing a few network nodes will ‘turn off’ some EET pathways, but other strongly coupled paths remain, leading to inherent robustness, as observed in figure 5. We note, however, that the mean path-length for the LHC-II monomer is always greater than FMO, which implies EET is generally less efficient. A similar observation can be made for the network diameters, although excitation at Bchl 1 in FMO shows more sensitivity to chromophore knock-out than excitation at Bchl 6. The LHC-II monomer, on the other hand, shows marked sensitivity beyond the removal of nine chromophores, an observation which may be understood based on the weak average inter-Chl couplings in LHC-II displays compared with those in FMO. The network diameter, like the mean path-length, highlights the presence of less strongly coupled Chls in LHC-II compared with FMO; in particular, LHC-II always has a larger network diameter than FMO. However, the conclusion from all of these sets of simulations is quite transferable; the average path-length and network diameter of both PPCs considered here are found to be reasonably insensitive to the removal of Bchl/Chl sites, up to a certain point. This relatively constant nature of EET paths through the network emphasizes the existence of multiple EET paths with comparable couplings; in short, our results are consistent with the idea of inherent network redundancy in both PPCs considered here, in line with our previous findings for the eight-site FMO model [56].

### (c) Robustness to environmental fluctuations

Until now, these simulations have considered changes to the network itself, through the loss of Bchl/Chl sites or inter-Bchl/Chl coupling elements. However, another source of perturbation is fluctuations due to the environment itself, for example arising due to thermal motion. This can be investigated here by altering the global dephasing rates (equation (2.9)) experienced at each site, and subsequently recalculating the EET efficiency; such simulations can give insight into how susceptible a given electronic Hamiltonian is to perturbations in the environment.

The results of these simulations are shown in figure 9. Notably, the width of the EET response curves for FMO (figure 9*a*) are much wider than for the LHC-II monomer (figure 9*b*). This can be explained by the average coupling strength of FMO being much stronger than in the LHC-II monomer (approx. 24 cm^{−1} versus approx. 8 cm^{−1}, as noted above). Both systems, in general, display broad response curves to a changing environmental dephasing rate; this response, which arises due to the interaction between exciton transport and trapping, as noted previously, demonstrates that these PPCs exhibit efficient EET across a broad range of dephasing rates, again demonstrating an inherent robustness to environmental fluctuations. Notably, the typical range of dephasing rates over which efficient EET is observed is consistent with the experimentally observed dephasing rates of typically 100 cm^{−1}. As an aside, we note that the broader environmental response curve of FMO is consistent with the fact that FMO complexes in GSB operate in a variety of thermal environments [18,21,99]; however, our simplified and artificial models of EET dynamics are not an appropriate tool to probe such questions further.

### (d) Implications for pigment-protein complexes-as-networks

The results presented above highlight an interesting connection to the theory of networks. In particular, the robustness properties demonstrated for FMO and LHC-II, whereby efficient EET is maintained despite significant removal of chromophore sites from the electronic sub-system, is comparable with the robust behaviour observed in ‘scale-free’ networks [85]. In such networks, the degree of connectivity *k* of each node follows a power-law distribution, *P*(*k*)≃*k*^{−α}; an important feature of such networks is the appearance of nodes which are connected to a large number of other nodes, albeit with most nodes connected to a lesser degree. Indeed, the results of figure 8 are strongly reminiscent of similar ‘knock-out’ analyses of model scale-free networks investigated previously [85].

While this connection to scale-free topology is certainly interesting, a word of caution is in order. In particular, we expect that the PPC networks considered here will not follow the power-law distribution of connectivity expected for scale-free networks; this is simply because scale-free networks exhibit a large-number of vertices with just a few (one or two) connections, whereas typical PPCs demonstrate structures in which very few chromophores are ‘singly connected’. In other words, most chromophores in PPCs, at least for those considered here, exhibit connections to more than one other chromophore site, such that the distribution of network connectivities is peaked around some value which is higher than one. This is confirmed in figure 10, which was calculated assuming that two chromophores are ‘connected’ if they have a coupling strength greater than 5 cm^{−1} in magnitude. Clearly, for both FMO and LHC-II monomer, we find the connectivity distribution is peaked at around four or five, with few singly connected nodes observed. So, the interesting conclusion is that, although PPCs exhibit scale-free-type robustness, their connectivity is somewhat different from the expected power law. Expanding our studies to other, larger PPCs will go some way to improving our understanding of this interesting feature of PPC systems.

## 4. Conclusion

In this article, we have expanded on our previous investigation of EET dynamics in eight-site FMO by considering similar network-based analyses of EET efficiency and robustness in a seven-site FMO model and in monomeric and trimeric LHC-II. Our new results are in line with our previous findings. In particular, we find that the PPC studied here exhibit a high degree of robustness to disruption of the electronic sub-system (chromophores). We find that removal of chromophores sites usually has only a small impact on calculated EET efficiency, up to some limiting point. By comparing targeted removal of chromophore sites to random removal of inter-chromophore coupling elements, as well as by considering network features such as average network path-length and network diameter, the picture which emerges is one in which PPCs exhibit several EET pathways between source and sink sites which exhibit comparable electronic coupling strengths along these paths. In other words, the PPCs considered here appear to demonstrate an inherent redundancy in EET pathways, a feature which is clearly beneficial in maintaining a robust and efficient system for controlling energy transport to the RC.

We have also noted that these robustness features of the PPCs considered here are comparable to scale-free networks, although the PPCs here do not really demonstrate a scale- free-type connectivity. The PPCs studied here are too limited in size and scope to truly investigate the connection between scale-free behaviour and EET dynamics; this is currently a goal we are working towards. Furthermore, it is clear that, while the HS model is essential in our work due to its computational simplicity, the approximations inherent in this model of EET dynamics could potentially modify the *quantitative* features of network robustness. Improved treatment of coupling between the electronic subsystem and the vibrational environment in PPCs will certainly modify EET dynamics, though the exact impact on robustness and EET efficiency is less clear; again this is an area which we are actively pursuing. This broader overview of EET dynamics, and perhaps better treatment of environment, will play an important role in transforming our initial studies here into tangible ‘rules’ for the design of new artificial light-harvesting systems.

## Data accessibility

Data for figures 4–9 can be found at http://wrap.warwick.ac.uk/85989.

## Authors' contributions

S.H and L.A.B designed the study. L.A.B performed the simulations. S.H and L.A.B wrote the paper.

## Competing interests

The authors declare no competing interests.

## Funding

The authors are grateful to the Engineering and Physical Science Research Council (EPSRC) for providing a studentship through the Molecular Organisation and Assembly in Cells Doctoral Training Centre (Grant No. EP/F500378/1), and to the Scientific Computing Research Technology Platform at the University of Warwick for providing computational resources. This research utilized Queen Mary's MidPlus computational facilities, supported by QMUL Research-IT and funded by EPSRC grant EP/K000128/1.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3782228.

- Received February 15, 2017.
- Accepted May 4, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.