Molecular dynamics simulations in photosynthesis

Abstract

Photosynthesis is regulated by a dynamic interplay between proteins, enzymes, pigments, lipids, and cofactors that takes place on a large spatio-temporal scale. Molecular dynamics (MD) simulations provide a powerful toolkit to investigate dynamical processes in (bio)molecular ensembles from the (sub)picosecond to the (sub)millisecond regime and from the Å to hundreds of nm length scale. Therefore, MD is well suited to address a variety of questions arising in the field of photosynthesis research. In this review, we provide an introduction to the basic concepts of MD simulations, at atomistic and coarse-grained level of resolution. Furthermore, we discuss applications of MD simulations to model photosynthetic systems of different sizes and complexity and their connection to experimental observables. Finally, we provide a brief glance on which methods provide opportunities to capture phenomena beyond the applicability of classical MD.

Keywords: Coarse-grained; Conformational switch; Light harvesting; Molecular dynamics; Photosynthesis; Thylakoid membrane.

Fig. 1

Schematic of the thylakoid membrane and of the light reactions. A simplified scheme of the linear and cyclic electron transfer pathways (in black and cyan solid arrows, respectively) is represented together with the involved photosynthetic subunits. Each complex is labeled in the figure, together with the charge transport pathways and the cofactors and proteins involved. The primary donors of PSII and PSI are represented in yellow. The structures for PSII, PSI, and Cyt b₆f are taken from plants (Qin et al. ; Wei et al. ; Malone et al. 2019), the one for ATP from a bacterium (Morales-Rios et al. 2015). The thylakoid membrane is shown in blue

Fig. 2

Simulation techniques with different spatio-temporal resolution illustrated using examples from photosynthesis research. With increasing system size, longer simulation times are required which entails a decrease in achievable accuracy. Typical ranges of system size in terms of number of non-hydrogen atoms and simulation time are indicated on the right. The main focus of our review is on atomistic and coarse-grained MD simulations; a glimpse on the other simulation techniques depicted here is provided in Sect. 4 at the end

Fig. 3

Bonded potential terms and force field resolutions. A Representation of the principal variables inside the bonded potential terms in the example of a betacarotene molecule: the variables are distance ($d_{\mathit{ij}}$), angle (${\vartheta }_{\mathit{ijk}}$), dihedral angle (${\varphi }_{\mathit{ijkl}}$), and improper dihedral angle (${\phi }_{\mathit{ijkl}}$), and the different indices i, j, k, l refer to different interaction sites. B A betacarotene molecular structure represented at different levels of resolution: AA (all-atom resolution), UA [united-atom resolution, in particular in the case of GROMOS (de Jong et al. 2015)] and CG [coarse-grained resolution, in particular in the case of Martini (de Jong et al. 2015)]. The arrows indicate how the atoms of a UA-type of structure are grouped (“mapped”) into a CG-one

Fig. 4

Workflow of an MD simulation, as described in Sect. 2.3. The sequential steps are reported together with the associated main points that need specific attention. In the inset, an example of a simulation box for LHCII embedded in a model membrane is reported with water in cyan, lipid membrane in gray, protein in magenta, Chls in green, and Cars in orange

Fig. 5

Schematic representation of selected enhanced sampling techniques. A Schematic visualization of steered MD. A force (black arrows) is applied to one carotenoid to steer its unbinding from the LHCII binding pocket into the membrane. B During umbrella sampling, a biasing potential represented by the dashed harmonic potentials fixes the system along the reaction coordinate. The solid line shows the original potential. C In metadynamics, the biasing potential is added in the regions which were already visited along the reaction coordinate. The dashed lines indicate the changing potential energy surface during time (from light to dark red). D In accelerated MD, the biasing potential is time-independent resulting in the potential depicted by the black dashed line. E During temperature replica exchange, the system is simulated at different temperatures. Configurations are exchanged if the selection criteria are met

Fig. 6

Examples of structural insights on the function of photosynthetic protein complexes via MD simulations. Selected results from the MD work of our groups are presented here: A Snapshots of diffusive entrance of PQ (left) and exit of PQH2 (right) of PSII reaction center. The red and green dashed circles indicate the Q_B and Q_C binding sites, respectively. Figure reproduced with permission from Ref. (Van Eerden et al. 2017b). B Relative MGDG density around an LHCII trimer in the stromal and lumenal leaflet of the thylakoid membrane (top). The average volume density (bottom, magenta surface) of the LHCII trimer in thylakoid membrane obtained from 60 µs of CG simulations show an hour-glass shape. Selected MGDG lipids are depicted whose cone shape fits well the protein shape (Thallmair et al. 2019). In C, an LHCII complex from plants is reported (color scheme as in Fig. 1). Solvent, membrane, and Chls’ phytol tails not shown for clarity. In Ref. (Liguori et al. 2015) a single Chl couple (Chl a611-a612) has been found to be a disordered domain (Sect. 3.2). This disorder is measured by the changes in the excitonic coupling between the two Chls, as reported in C for one of the simulations in Ref. (Liguori et al. 2015). Such changes depend on the different conformations and organizations experienced by Chl a611-a612 along the simulated trajectory, as shown with different colors in C. In D, the pH-dependent conformational change at Helix 3 of PsbS observed by CpHMD in Ref. (Liguori et al. 2019) is reported, as described in Sect. 3.2