Predicting the future of excitation energy transfer in light-harvesting complex with artificial intelligence-based quantum dynamics
Exploring excitation energy transfer (EET) in light-harvesting complexes (LHCs) is essential for understanding the natural processes and design of highly-efficient photovoltaic devices. LHCs are open systems, where quantum effects may play a crucial role for almost perfect utilization of solar energy. Simulation of energy transfer with inclusion of quantum effects can be done within the framework of dissipative quantum dynamics (QD), which are computationally expensive. Thus, artificial intelligence (AI) offers itself as a tool for reducing the computational cost. Here we suggest AI-QD approach using AI to directly predict QD as a function of time and other parameters such as temperature, reorganization energy, etc., completely circumventing the need of recursive step-wise dynamics propagation in contrast to the traditional QD and alternative, recursive AI-based QD approaches. Our trajectory-learning AI-QD approach is able to predict the correct asymptotic behavior of QD at infinite time. We demonstrate AI-QD on seven-sites Fenna-Matthews-Olson (FMO) complex.
Introduction<!>Reference quantum dynamics of the FMO complex<!>Parameters-based non-recursive training framework<!>Input Target values<!>Discussion<!>Training data<!>Data preparation<!>Training architecture<!>Input normalization and redundant time-functions
<p>From the birth of life, solar energy has been the driving force of life. Via the mechanism of photosynthesis, living organisms capture sunlight with the highly sophisticated pigments in their antenna systems and transfer sunlight energy to the reaction center (RC) in the form of electron-hole pairs (excitons), where it is stored as biochemical energy. 1 The transfer of solar energy from antenna to RC, which is also known as excitation energy transfer (EET), in the form of excitons is considered to be highly efficient with close to unit efficiency. 2 Understanding this high efficiency of the natural harvesting systems is very important because of its potential applications in designing very efficient organic solar cells and storage devices. 3 Experiments showed that the long-lasting coherence in the efficient natural lightharvesting harvesting complexes (LHCs) is preserved by the surrounding protein environments (scaffold), and this coherence may be responsible for this high efficiency. 4,5 The most well-investigated LHC is Fenna-Matthews-Olsen (FMO) complex, which is found in green sulfur bacteria. 6 The small size and simplicity of the FMO complex also makes it a testbed of simulation approaches. The FMO complex is a trimer of identical subunits, where each subunit consists of bacteriochlorophyll (BChl) molecules (system) attached to their protein environments. 7 Enormous amount of research work has been done on light harvesting processes. [8][9][10][11][12][13] Taking FMO as an example, it is easy to see that the system (BChl molecules) is not isolated from the environment (the protein) and thus, the correct simulation of FMO should treat it as an open system rather than isolated one. In addition, many experiments suggest, 14,15 that quantum effects, particularly coherence, might play an important role in the light harvesting processes and may even be responsible for achieving the high-end efficiency. Temporal and spatial simulation of EET with the inclusion of quantum effects can be done within many frameworks such as classical mapping-based approaches, [16][17][18] perturbative methods, [19][20][21] and dissipative quantum dynamics (QD) [22][23][24][25][26][27] adopted here. QD simulations can be performed using the hierarchical equations of motion (HEOM) 28 and its many improvements and extensions, 8,23,[29][30][31] the quasiadiabatic propagator path integral (QuAPI) 32 and its variant iterative QuAPI (iQuAPI), 27 the trajectory-based stochastic equation of motion (SEOM) approach, 25,[33][34][35][36][37][38][39] the multi-layer multi-configuration timedependent Hartree (ML-MCTDH) 26 and the local thermalising Lindblad master equation (LTLME). 22 The development of various quantum dissipative dynamics methods stirs from the fact that each of these methods has some limitations and hence there is no single universal method that works in all cases. For instance, HEOM is numerically exact but come with a very high computational cost at low temperatures, the SEOM has no explicit dependence on the temperature but has very bad convergence at long-time propagation, in the QuAPI approach all correlation effects are included over a finite time and correlation effects beyond this time are neglected. Most importantly, all these traditional QD approaches require step-wise propagation of trajectories and the next step depends on the previous steps, thus, QD simulation is an iterative, recursive process. Both calculations at each time step and recursive nature of QD makes it rather computationally expensive.</p><p>Alleviating the computational cost of QD became a target of a series of studies applying artificial intelligence (AI), [40][41][42][43][44][45][46] inspired by advances in application of AI employing machine learning (ML) algorithms in computational chemistry and chemical physics. 47,48 AI was also applied to investigate EET in a dimer system 44 and the FMO complex. 40 Saving of computational cost by AI in above studies is impressive, however, one of the studies 40 only focused on predicting energy transfer times and transfer efficiencies rather than temporal and spatial evolution, while other related studies [44][45][46] adopted basically the same recursive nature of QD trajectory propagation.</p><p>The recursive nature of the previous AI-based QD makes it prone to error accumulation.</p><p>In recursive simulations, previously predicted values are used as an input to predict the next value. Thus, the prediction error at each time-step will accumulate, which results in deterioration of accuracy. In addition, the recursive nature of predictions does not allow us to make a prediction for any arbitrary time without predicting values before that. Finally, a short-time trajectory is needed as the seed to be generated with traditional approaches such as HEOM and then provided as an input to AI model to make prediction for the next time step and ultimately propagate the long-time dynamics. Thus, even when having AI model, we still need to spend valuable computational time to generate the short-time trajectory with the traditional approaches.</p><p>Here, we suggest an AI-QD approach to directly predict QD with AI as a function of time and other parameters such as temperature, reorganization energy, etc., completely circumventing the need of recursive step-wise dynamics propagation in contrast to the traditional QD and alternative, recursive AI-based QD approaches. Our AI-QD approach is able to predict QD at infinite time with correct asymptotic behavior and can be viewed as trajectory learning, which does not need any short-time trajectory as an input, eradicates the need of traditional approaches to generate the seed, and alleviates the problem of error accumulation. We demonstrate the applicability of AI-QD on seven-sites Fenna-Matthews-Olson (FMO) complex and show how AI-QD can be used for massive, infinite-time QD simulations and provide insights into the desired range of parameters and more efficient paths followed by the transfer of excitation energy.</p><!><p>We employ the Frenkel exciton Hamiltonian 49 to study EET dynamics in the FMO complex:</p><p>with all Hamiltonian terms given below</p><p>where H s , H env , H s−env and H reorg denote system (BChl molecules) Hamiltonian, Hamiltonian of protein-environment, system-environment interaction Hamiltonian and the reorganization term, respectively. In Eq. ( 1), n is the number of sites (BChl molecules), ϵ i is the energy of the ith site and J ij is the inter-site coupling between sites i and j. P k,i and Q k,i are respectively momentum coordinate and frequency of environment mode k associated with site i. In H s−env , each site is connected to its own environment. The c k,i is the strength of coupling between site i and mode k of its environment. The reorganization term H reorg can be seen as a counter term that emerges from the interaction of the sites with the environment. 8,49,50 It is added to stop further renormalization of the site energy ϵ i by the environment. In the reorganization term H reorg , λ i is the reorganization energy corresponding to site i, 51</p><p>where J i (ω) is spectral density of the environment corresponding to site i. As shown by Nalbach and Thorwart, 52 the effects of the discrete molecular modes on the population dynamics are largely irrelevant. As a result, it is acceptable to use continues environment spectral density such as Drude-Lorentz spectral density</p><p>where γ and λ denote the characteristic frequency (bath relaxation rate) and the reorganization energy, respectively.</p><p>In general terms, the EET dynamics in the FMO complex can be described by Liouville-</p><p>where ρ is the density matrix. Because of the many-body effects, direct propagation of Eq. ( 8) is not straightforward. Different approaches are developed to simplify and propagate Eq. ( 8) and interested readers are advised to look into the corresponding references. 25,30,32,53 We use the local thermalising Lindblad master equation (LTLME) 22 to propagate the reference QD trajectories for the reduced density matrix of the system (see the Supplementary Information), where we adopt Adolphs and Renger's Hamiltonian for seven sites per subunit 54 (see Methods). The LTLME is a coherent and complete positive trace-preserving approach, but may not be that accurate as HEOM or SEOM approaches (because of approximations used in its derivation 22,55 ), but here it is not the concern of our proof-of-concept paper.</p><!><p>In our parameters-based non-recursive AI-QD, we train ML model as a function of a parameter space D (used as the input to ML model) which depends on the system of interest and on the data from a limited number of QD trajectories. For the FMO complex, our parameter space D consists of information of sites, λ, γ and T . In addition, time also becomes a part of the input of our AI-QD model.</p><!><p>In addition, labels n = {n 1 , n 2 , n 3 , . . . n 7 } are used for corresponding rows in the density matrix and labels for sites with possible initial excitation are m = {m 1 , m 2 }. As the offdiagonal elements ρ pq, p̸ =q are complex, we separate the real and imaginary parts. Our training trajectories generated with the reference LTLME-QD approach are chosen by furthest-point sampling from the three-dimensional space of the following parameters: re-</p><p>We should also decide up to what time-length t M we should run reference LTLME-QD trajectories. Based on the prior knowledge that populations plateau in asymptotic limit, for each trajectory we choose a different time-length t M using a vanishing gradient scheme, where t M is chosen such that the gradient of population G is close to zero (see Methods).</p><p>Using the vanishing gradient scheme to find different t M for each trajectory allows us to sample more data from the training trajectories, which are hard-to-learn, while avoiding redundant sampling from trajectories, which are easy-to-learn. This also removes arbitrariness in choosing fixed t M parameter as was done in previous studies using the recursive AI-QD scheme. 44,46 Application to EET dynamics in FMO complex</p><p>As an application of our approach, we predict EET dynamics in the FMO complex with seven sites per subunit for parameters of the test set trajectories (none of which used in training).</p><p>Site-1 (BChl molecule 1) and site-6 (BChl molecule 6) are most likely to get initially excited as they are close to the photosynthetic antenna complex called chlorosome, 6 we thus present results for both cases. For predictions, we just provide the parameters of the test trajectories (characteristic frequency, reorganization energy, temperature) as an input and predict the evolution of EET. Fig. 2 shows the evolution of excitation energy in all seven sites for both cases. In Fig. 2, we show EET for both short and long time periods, demonstrating that AI-QD is able to capture the coherent EET (aka quantum beating or modulation of amplitudes) of short-time dynamics and also can predict the asymptotic limit. Fig. 3 shows the prominent off-diagonal terms (aka coherence) of the reduced density matrix for Fig. 2.</p><p>Table 1 shows mean absolute error (MAE) and root mean square error (RMSE) averaged over 600 trajectories. As AI-QD is non-recursive (non-iterative), without any trajectory propagation, we can directly predict the asymptotic behaviour. Our AI-QD performs well in all cases (from weak coherence to strong coherence, from Markovian to non-Markovian, from adiabatic to nonadiabatic situations) as can be observed for selected trajectories shown in γ and λ are in the units of cm −1 , while T is in the units of K.</p><p>Fig. S1 with corresponding errors reported in Table S1 of the Supplementary Information.</p><p>From Table S1, we observe that our AI-QD approach is comparatively more accurate in strongly coherent cases (large value of γ and small values of λ and T ) which can be seen as a consequence of the vanishing gradient scheme which may favor these challenging cases due to a larger number of training points sampled from such trajectories. AI-QD approach can even extrapolate to a good degree as its error for the test trajectories propagated with parameters outside the training parameter space is of a similar order of magnitude to the test as shown in the Supplementary Information (see additional set of test trajectories in Fig. S2 and Table S2).</p><p>It was shown, 8,56,57 that the transfer of excitation energy in the seven-sites FMO complex follows mainly two paths, i.e., site-1 → site-2 → site-3 ↔ site-4 and site-6 → site-5, site-7, site-4 → site-3, here the ↔ shows that the excitation energy equilibrates between site-3 and site-4 after site-3 is populated (see Fig. 3) Among the seven sites, the sites 1 and 6 are close to the baseplate protein, while the sites 3 and 4 are near to the target RC complex. 54,58 It has been proposed that the quantum coherence allows the FMO complex to quickly sample several routes (paths) in search for site-3. 5 In Fig. 4, we show the population of site-3 at t = 0.5 ps (500 fs) as a function of γ, λ and T . From Fig. 4(a), we observe that at room temperature T = 300, the ETT to site-3 or, in other words, to RC complex gets slow as the characteristic frequency γ increases. In contrast, the ETT to site-3 increases with the increase in reorganization energy λ as shown in Fig. 4(b). Similar trend can be observed with the increase in temperature T as can be seen in Fig. 4(c).</p><p>In order to find the optimum parameters for the fastest transfer of excitation energy, we From Eq. ( 9) (Methods), energy of the site-1 (12410 cm −1 ) is lower than the baseplate, which has been reported to be 12500 cm −1 . 59,60 This allows a quick transfer of the excitation energy to site-1 from the baseplate. However, the energy of site-2 (12530 cm −1 ) is higher than site-1 and also than site-3 (12210 cm −1 ), which on the one hand stops backward transfer from site-3, but on the other hand creates a local minimum on site-1. Despite the local minimum on site-1, the excitation energy is not trapped because of the quantum coherent wave-like motion between site-1 and site-2. Following path-2, the energy of site-6 (12630 cm −1 ) is higher than the energy of baseplate. To stop backward transfer of excitation energy from site-6 to baseplate, site-6 should quickly transfer excitation energy to other sites such as site-5, site-7 and site-4. This quick transfer from site-6 to site-5, site-7 and site-4 is only possible by the strong coupling of site-6 to site-5 and site-7, which in return are strongly coupled to site-4. T . The blue line corresponds to the case with initial excition on site-1 while the red line is for the case with initial excition on site-6. γ and λ are in the units of cm −1 while T is in the units of K.</p><!><p>In this work, we have presented a non-recursive (non-iterative) AI-QD approach for blazingly fast prediction of quantum dynamics, as predictions can be made for any time step up to asymptotic limit completely circumventing the need of recursive trajectory propagation.</p><p>This can be used, as we demonstrated here, for massive quantum dynamics simulations, for example, in search for the best conditions required for efficient energy transfer in designed photovoltaic devices. Just to put things into perspective, our AI-QD approach can predict the entire 2.5 ps trajectory within ca. 2 min on a single core of Intel(R) Core(TM) i7-10700 CPUs @ 2.90 GHz, independent of the reference method used for generating training trajectories, while the same propagation with the traditional recursive approaches such as HEOM would take hours, and the cost would exponentially increase for low temperatures. The high cost of accurate approaches such as HEOM was also a reason why we used a much faster LTLME for this proof-of-concept study to extensively test our approach (propagation of an entire trajectory takes only 3 min with LTLME on a single CPU of the above computer architecture.</p><p>It is worth emphasizing that AI-QD is embarrassingly parallel and the calculations can be further significantly sped up by using multiple CPUs or GPUs, because predictions with AI-QD for different time steps are independent of each other and different segments of trajectories can be distributed for independent calculations on many threads.</p><p>We demonstrated the feasibility of AI-QD approach on an example of the FMO complex, but this approach is general enough to be used for any other complex after retraining. It remains to be seen how well the AI-QD approach can be extended to describe several LHCs However, open question remains how successful would be such an approach to learn dynamics and in addition, how to circumvent different dimensionalities of Hamiltonians of different complexes.</p><!><p>In the seven-sites FMO complex (apo-FMO), where seven BChl molecules (seven sites) exist per subunit, the inter-subunit interaction is very small and each subunit can be considered relatively isolated. 61 Here we adopt Adolphs and Renger's Hamiltonian for seven sites per subunit 54</p><p>where energies are given in cm −1 . Each site is coupled to its own environment characterized by the Drude-Lorentz spectral density given by Eq. ( 7). Not long-ago, an eighth BChl molecule (site-8) has been discovered, 11 however as has been mentioned by Jia et al., 62 the role of the eighth BChl molecule (site-8) in the transfer of excitation energy in the FMO complex is negligible.</p><p>Trajectories for the reduced density matrix have been generated with the local thermalising Lindblad master equation (LTLME) 22 (see the Supplementary Information) implemented in quantum_HEOM package 63</p><!><p>With all the possible combinations of the parameters (D), we have 3960 total number of trajectories N traj (1980 (site-1) + 1980 (site-6), all these trajectories correspond to their respective combination of parameters in parameter space D). Using farthest-point sampling 65</p><p>in the three-dimensional space of λ, γ and T , we choose 1000 trajectories as our training space T S (500 (site-1) + 500 (site-6) , ca. 25% of space D)), 200 trajectories as the validation set VS (ca. 5% of space D)) and the rest of trajectories, we keep as the test set ST P (ca. 70% of space D). For each trajectory, we choose a different time-length t M using a vanishing gradient scheme. In this scheme, we take the gradient G of the population of each site (ρ nn , n = 1, 2, 3, . . . , 7) for 10 consecutive time-steps and if all of them remain less than the threshold value of G th = 1 × 10 −10 , we choose our t M . We find t M for all seven sites and then choose the maximum value among them, thus we keep a single value of asymptotic limit (t M ) for all seven-sites. By analyzing the gradients, we find the region of the trajectory, where the change in population of the site is very small. By knowing that, we keep the time-length of our trajectory t M up to that region, because beyond t M the change in population is very small, and ML is able to predict it. As the asymptotic limit for each trajectory is different, we have different value of t M for each trajectory. In our training, we have included t → ∞, corresponding to the asymptotic behaviour at long-time.</p><!><p>We use convolutional neural network (CNN) architecture, because the importance of convolutional layers is much explored for image analysis, where these layers extract important features such as edges, textures, objects, and scenes. When it comes to time-series data, we Non-trainable parameters: 0 are using convolutional layers in the hope to extract some important features from the data (such as the time influence). After learning those features, when we provide a test trajectory, the trained ML model will look for those features in that test trajectory. 66 Though we have used the CNN model, other ANN architectures such as long short-term memory (LSTM) is also an option. LSTM is considered to be more suitable for extracting long-time temporal dependencies in contrast to convolutional neural networks (CNNs) which are more local.</p><p>However, CNNs are easy to train and in many studies, they have outperformed LSTM for future forecasting.   1.</p><!><p>As we have multiple input elements, we need to normalize them all. In normalized input, we</p>