Exploring channelrhodopsin-2 and gold nanocluster interaction: a route to control the protein photocycle

Abstract

Hybrid light–matter platforms that combine biomolecular photochemistry with engineered electromagnetic fields offer a promising, largely unexplored way to extend optogenetic control beyond the limits of native chromophores. Channelrhodopsin-2 (ChR2), a widely used light-gated ion channel, serves as an ideal model to test whether metallic nanoclusters can interface with membrane photoreceptors without disrupting their structure or function. Using large-scale molecular dynamics simulations from multiple starting configurations, we simulate ChR2 in complex with a gold nanocluster, which does not exhibit plasmonic behavior but lies within the size regime where metallic nanoclusters begin to display plasmonic resonances. We find that these functionalized nanoclusters spontaneously form stable, physiologically relevant complexes with ChR2 while preserving its global structure, hydration, retinal binding pocket, and native dynamics. Protein network and contact perturbation analyses uncover an unreported allosteric link between nanocluster binding sites and the retinal cavity, suggesting a mechanism to tune photochemical behavior. Our results establish a molecular feasibility framework for integrating nanoclusters with optogenetic receptors, opening the way for the design of future nano-optogenetic platforms that could modulate photoreceptor photochemistry without genetic or chemical modification of the photoreceptor itself.

---

1 Introduction

The ability to manipulate neuronal activity with light has reshaped neuroscience, positioning optogenetics as a central technology for circuit interrogation and therapeutic innovation.^1–5 Yet, current optogenetic actuators remain constrained by the photophysical limits of their native chromophores and the restricted ability to externally modulate their excited-state landscapes. Hybrid light–matter systems that couple biomolecular photoactivity to engineered optical near-fields offer a promising, but largely unexplored, path to overcome these limits.

Channelrhodopsin-2 (ChR2), a light-sensitive ion channel protein naturally found in green algae, is a key player in optogenetics.^6–9 Its structure consists of two monomers, each composed of seven transmembrane (TM) helices connected by extracellular and intracellular loops. The N-terminal region of the protein is located outside the cell, interacting with the external environment, while the C-terminal region is located inside the cell, interacting with the cytoplasm. The light-sensitive part of ChR2 is retinal, a chromophore bound to the protein via a covalent bond with a specific lysine residue in TM-7 via a protonated retinal Schiff base (PRSB). Upon light absorption, the retinal undergoes photoisomerization from its all-trans form to its 13-cis form.¹⁰ Isomerization represents the initial event in the photocycle, a series of molecular processes that occur in response to light. These include proton transfer reactions and conformational changes, which result in the opening of the channel allowing ions to flow across the membrane.¹¹ ChR2 has proven indispensable for millisecond-scale neural control; however, there is still no means to externally control retinal photoisomerization, which remains fundamentally dictated by the native protein environment.

Since its inception, optogenetics has frequently employed metal nanoparticles, particularly gold nanoparticles (AuNPs),¹² primarily to exploit plasmonic effects for signal enhancement.¹³ Interestingly, localized plasmonic fields have been shown to influence the bacteriorhodopsin photocycle through near-field enhancement of the local electric field, which perturbs the potential energy surface of charge-redistributed retinal intermediates,¹⁴ and resonant enhancement of light absorption at key intermediates, which accelerates photocycle kinetics resulting in a significantly amplified photocurrent.¹⁵ Importantly, plasmon-induced effects are not restricted to larger nanoparticles but can also emerge in smaller metallic nanoclusters (NCs), which can reshape photochemical energy landscapes via strong and ultrastrong light–matter coupling,^16–18 as evidenced by recent advances in the emerging field of polaritonic chemistry.^19–28 Due to their unique electronic/optical properties, metallic NCs have already demonstrated significant potential across a wide range of applications such as catalysis,^29–31 biomedical imaging,^32,33 medical therapy,^34–36 etc. Combining metallic NCs with biomolecules like proteins has been shown to create synergistic effects, uniting the distinctive optical, electronic, and catalytic properties of the metal NCs with the inherent biological functions of the proteins.^37–39 Therefore, metallic NCs seem promising for modulating the photocycle of light-sensitive proteins under the strong coupling regime. Yet, it is still unexplored whether these optical antennas can interface productively with a membrane ion channel without ablating the structural integrity required for gating. Structural alterations of the protein host and/or uncontrolled formation of clusters at unintended locations within the protein can in fact limit the application of these materials in biomedical research.⁴⁰

Here, we use large-scale molecular dynamics simulations across 31 independent trajectories to investigate the interaction of ChR2 with a functionalized Au₂₅ nanocluster.^41,42 Although Au₂₅ does not exhibit plasmonic behavior, it falls within the size regime where metallic NCs already display plasmonic resonances.^43–46 Unlike larger plasmonic objects, however, Au₂₅ can diffuse through cellular membranes. Au nanoparticles with a 2 nm radius and hydrophobic surfaces have already been successfully incorporated into the leaflets of bilipid vesicles.^47,48 This makes Au₂₅ a suitable model system to investigate nanocluster–protein interactions from a mechanistic perspective. Furthermore, the Au₂₅ nanocluster used here is functionalized with phenyl groups, promoting interaction with the hydrophobic region of the lipid bilayer and the hydrophobic surface of the membrane protein. We analyze the propensity and stability of the interaction, identify the most probable binding sites, and assess the effects of the NC on the global protein structure, channel hydration, and key retinal interactions.

Our findings establish a molecular feasibility boundary for integrating plasmonic nanostructures, similar to the NC analyzed here, with intact optogenetic receptors, marking a foundational step toward hybrid nano-optogenetic systems in which metallic nanostructures could regulate photoreceptor excited-state dynamics – and ultimately neural activation – through engineered nanoscale electromagnetic environments.

2 Methods

2.1 Molecular dynamics simulations

We perform MD simulations of ChR2 dimeric assembly (PDB ID 6EID⁵⁰) in the dark state conformation (i.e., with retinal in the all-trans state) both alone and in interaction with a hydrophobic gold thiol-protected NC for which the X-ray structure is known.⁴¹ Specifically, the NC employed consists of a Au₂₅L₁₈⁻ cluster, where L = (S(CH₂)₂Ph) and Ph = C₆H₅ (phenyl group). It comprises a Au₁₃-centered icosahedron, with six Au₂L₃ units arranged around it on the 12 faces of the core. The NC exhibits a net charge of −1.00, is diamagnetic, and displays poor solubility in water. Previous simulations of its interactions with the protein β₂-microglobulin,⁴² showed that this NC is able to bind the protein, forming persistent complexes. Previous simulations by some of us showed that glutamic acid (Glu) 90, a key residue in the ChR2 photocycle involved in gating and ion selectivity, can adopt two different hydrogen bonding patterns in the dark state (referred to as DOWN and UP, see section 3.2).⁵¹ This observation is also in agreement with previously performed MD simulations of ChR2.^52,53 Therefore, we selected a structure extracted from our previous MD simulations⁵¹ as the starting structure for all MD simulations, in which one monomer is in the UP conformation and the other in the DOWN conformation. With this starting structure, which also included the phosphatidylcholine (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine, POPC) membrane and the solvent, we performed two 500 ns-long MD simulations of ChR2 alone (without the NC) to be used as a reference. Using the same coordinates for ChR2, we also built the starting structures for 26 MD simulations of the protein in interaction with the NC as follows. We obtained 18 starting structures by creating a 3 × 3 grid on the x–y plane of the box, with 9 points on the extracellular side and 9 points on the intracellular side. The NC was then inserted at the desired position within this grid at a distance of ≈3 nm from the protein by removing an appropriate number of water molecules (see Fig. S1A and B in the SI). The remaining 8 starting structures were obtained by placing the NC on either the extracellular or the intracellular side immediately above the membrane and in close proximity to the protein (see Fig. S1C–F in the SI and section 3.1 for the rationale of this choice). Initially, the 26 simulations were run for 50 ns. Then, the 12 simulations in which the NC interacted with the protein were extended to 500 ns together with 3 simulations in which the NC entered the membrane. The analysis of these 26 simulations shows that the NC interacts with different regions of ChR2 with varying propensities (see section 3.1). Specifically, we observe that in 7 out of the 26 instances, the NC interacts with the C-terminal region of ChR2 (IvC) while it interacts with the N-terminal region (IvN) and with the membrane-embedded region (IvM) only in 2 and 3 instances, respectively. Therefore, to improve our statistics and ensure a minimum of 5 simulations for each interaction region, we added 3 and 2 additional MD simulations with the NC interacting with the N-terminal and membrane-embedded regions of ChR2, respectively. The starting structures for these additional simulations were taken from existing IvN and IvM trajectories at t = 20 ns, and new initial velocities were assigned to ensure independent sampling in the new MD runs. These 5 simulations performed to increase the statistics of the IvC and IvN interactions were also run for 500 ns. A table summarizing the initial NC position and site of interaction between the protein and the NC is provided in the SI (Table S1).

For all 33 MD simulations (2 without the NC and 31 with the NC), the POPC-embedded ChR2 dimer is solvated with TIP3P water molecules⁵⁴ and neutralized with Cl⁻/Na⁺ with a salt concentration of about 0.15 mol l⁻¹ in a simulation box of ≈(98 × 98 × 140) ų, for a total of more than 140k atoms. Amber99sb parameters^55,56 are used for the protein (including the retinal^57,58), membrane^59,60 and counterions and previously QM-derived parameters are used for the NC.⁴² Before proceeding with the productive run, a preliminary energy minimization and subsequent relaxation for 1 ns are performed for all systems after the NC insertion. All MD simulations are performed in the NPT ensemble (T = 310 K and P = 1 bar) maintaining the temperature and pressure constant with the V-Rescale⁶¹ and Parrinello–Rahman⁶² algorithms, respectively. Periodic boundary conditions and the Particle–Mesh–Ewald algorithm are used.⁶³ A 12 Å cutoff is used for van der Waals interactions. The LINCS algorithm⁶⁴ is used to constrain all covalent bonds involving hydrogen along with a 2 fs integration time step. All simulations were performed using the GROMACS 2022.3 package.⁶⁵ In summary, 33 independent trajectories were generated (2 of ChR2 alone and 31 in the presence of the NC). Of these trajectories, 22 were initially run for 50 ns and then extended to 500 ns while the remaining 11, where no interaction between ChR2 and the NC occurred, were run for 50 ns. In total, these 33 MD simulations yielded a cumulative simulation time of more than 11 μs.

2.2 Dynamical perturbation contact network and connected component analysis

In order to achieve a more profound comprehension of the local conformational changes resulting from the interaction of the NC, we employed the Dynamical Perturbation Contact Network (DPCN) method. Some of us have previously used this method successfully to study allosteric pathways, and it can also be used to make comparisons between different MD simulations.⁶⁶

Briefly, following the procedure detailed in Gheeraert et al.,⁶⁷ for each frame of molecular dynamics simulations, we determined which atoms were in contact based on a distance cutoff of 5 Å, considered a robust choice for protein structure networks.^67–69 These contacts were used to build an atomic contact matrix for each simulation frame. The atomic contact matrix was averaged across the entire simulation trajectory. This average atomic contact matrix was then transformed into a residue contact matrix, which describes interactions between amino acid residues rather than individual atoms. To focus on inter-residue interactions, we removed self-contacts. To identify changes induced by perturbations, we subtracted the average residue contact matrix of the perturbed state (simulation with the NC) from that of the reference state (the average of two simulations without the NC), creating a perturbation contact matrix. This matrix was used as an adjacency matrix for the calculation of the properties of the DPCN and its elements correspond to the weights of the edges of the network. For visualization purposes, we added a coloring scheme to the edges: blue if the weight is bigger in the reference state and red if the weight is bigger in the perturbed state.

We further analyzed the connectivity of these networks using Connected Component (CC) analysis.⁶⁷ Two residues are considered connected if there is a path between them within the network. A connected component is a subgraph in which all residues are mutually connected and are not part of a larger connected subgraph. While protein networks often form a single connected component, removing weak interactions by applying a threshold can lead to disconnected networks with multiple connected components. Connected component analysis, thus, simplifies the interpretation and visualization of complex DPCNs.

DPCN analyses were performed using the Python PMDlearn package.⁷⁰

3 Results and discussion

3.1 Assessment of the ChR2–NC interaction tendency

To understand if the NC can spontaneously interact with ChR2, we started 18 MD simulations placing the NC at ≈3 nm from the protein. In 9 initial structures, the NC was positioned on the intracellular side and in the other 9 initial structures, it was positioned on the extracellular side (see section 2 and Fig. S1A and B in the SI). The analysis of these trajectories showed that the NC could either interact with the protein or enter the membrane. Within 25 ns, in 7 out of 18 simulations (≈40%), the NC started interacting with the protein (average time before the interaction: 7.5 ± 6.9 ns) and maintained the interaction for the rest of the MD simulations. In the other simulations, the NC entered the membrane. It should be noted that these conditions could still lead to an interaction between the NC and ChR2 because the NC could approach the protein by spontaneous diffusion within the membrane. Since diffusion in the membrane is expected to be a slow process, we mimicked the above scenario by performing 8 additional MD simulations in which the starting structures were obtained by placing the NC immediately above the membrane and in close proximity to the protein, leaving only a thin layer of water molecules between them. These initial positions promoted a membrane-mediated NC–protein interaction, in which the hydrophobic NC penetrated the membrane but remained in close proximity to the protein and interacted with it. In 5 of these 8 MD simulations, the NC interacted with the protein. However, the desired interactions (insertion in the membrane in close contact with the protein, see Fig. 1) only occurred in three of them. In the other 2 MD simulations, the NC interacted with the protein similarly to what was observed upon spontaneous diffusion (vide infra). Overall, the interaction between the NC and ChR2 occurred in 12 MD simulations out of 26 (i.e., 46%). In what follows, we classify the 12 MD simulations in which the NC interacted with ChR2 according to the interaction site (see Fig. 1): interacting via the C-terminal (IvC), interacting via the N-terminal (IvN), and interacting via the membrane (IvM). The analysis of these 12 simulations revealed that the NC interacted with the C-terminal in 7 instances, while interactions with the N-terminal occurred in 2 cases. The remaining three interactions were mediated by the membrane. Furthermore, the NC demonstrated an equal propensity to interact with either monomer 1 or monomer 2. These results suggest that the protein–NC interaction was favored when the nanocluster diffused toward the C-terminal part (intracellular). This could be ascribed to the fact that the C-terminal is rather flexible and is therefore able to undergo a structural reorganization to accommodate the NC. In contrast, the N-terminal region (extracellular) is relatively rigid because of the presence of two disulfide (SS) bridges (involving cysteine, Cys-34 and Cys-36) connecting the two monomers. To compensate for the difference in sampling between IvC and IvN/IvM interaction sites, due to the varying interaction propensities of the NC with these sites, additional simulations were conducted to enhance the statistical representation of IvN and IvM. Thus, in order to have at least five simulations for each group, three simulations were added to the IvN group and two simulations were added to the IvM group. A table summarizing the initial NC position and site of interaction between the protein and the NC is provided in the SI (Table S1).

Fig. 1 Representation of NC-ChR2 interaction sites: IvC, interacting via the C-terminal; IvN, interacting via the N-terminal; IvM, interacting via the membrane; and NoI, non-interacting. The NC is represented in ball-and-stick format (orange, gold atoms; yellow, sulfur atoms; cyan, carbon atoms; white, hydrogen atoms). ChR2 is represented in cyan. For clarity, the POPC membrane and the solvent are not shown.

Overall, the MD simulations performed showed that the hydrophobic NC was able to interact spontaneously and effectively with the protein in approximately half of the MD simulations performed. In the other cases, a favorable interaction between the NC and the membrane was observed that might lead, by diffusion through the membrane, to further interactions with ChR2.

3.2 Characterization of ChR2–NC interactions

To highlight the protein regions and specific amino acids involved in the interaction with the NC, we computed a contact map. To this aim, we define an existing contact between residue i of ChR2 and the NC at time t when the distance between any atom of residue i and any atom of the NC is less than 4 Å. This distance is chosen because, given the nature of the NC, we expect that the interaction is mainly driven by hydrophobic effects.

In the IvC MD simulations (Fig. 2A), the involvement of the C-terminal residues was clearly present. Nonetheless, the involvement of other protein regions (highlighted in green in Fig. 2A) was also observed. These regions correspond to intracellular loops 1 and 2 (ICL-1 and 2) and helix TM-5. Notably, this interaction pattern was highly consistent across all simulations. Furthermore, the binding profile appeared to respect the symmetry of the dimer, as the NC tended to interact with structurally equivalent residues on both monomers. This recurrence suggests a preferential binding motif.

Fig. 2 Color map of the contacts between ChR2 and the NC. In panels A, B and C, the results from the IvC (A), IvN (B) and IvM (C) MDs, respectively, are reported. The X in the label Iv(C/N/M)_X identifies the simulation. The color scale (from white to dark green) indicates the ratio of the number of frames in which a residue is in contact with the NC (Nf_i) to the total number of frames in the simulation (Nf_tot). The two monomers are separated by a dotted red line. Below each panel a representative snapshot of ChR2 is reported in which the protein regions in contact with the NC are highlighted in green (most frequent contacts in dark green). For the sake of clarity, five simulations for each interaction site are reported. The two additional IvC MD simulations are shown in Fig. S2 in the SI.

This was also suggested by the observations that the same interaction pattern was present in 2 out of the 8 simulations whose starting structure was chosen to promote a membrane-mediated interaction (see above). Note that, for the sake of clarity, only five out of seven simulations of the IvC group are reported in Fig. 2A. The other two are shown in Fig. S2 in the SI and feature the same interaction pattern. Fig. 2B shows that the five MD simulations in which the NC interacted with the N-terminal featured a slightly different interaction pattern. In the IvN_1,3,4 MD simulations, the NC fit into the extracellular region of one monomer, with the involvement of the N-terminal and extracellular loops 2 and 3 (ECL-2 and 3). In the IvN_2,5 MDs, the NC strongly interacted with the SS bridge, as also shown by the involvement of the first residues of both monomers. In Fig. 2C, we report the contacts observed in the MD simulations in which a membrane mediated interaction occurred. In this instance, the interaction pattern was contingent upon the height and the area where the interaction occurred. Nonetheless, regardless of the monomer involved, a preference for the protein region where helices TM-5 and TM-6 are located was observed in the IvM simulations.

It was possible to extract from the contact map the amino acids that were more frequently involved in interactions with the NC. In Fig. 3, the amino acids involved in interactions with the NC were ordered according to decreasing hydrophobicity (highly hydrophobic residues in dark blue and charged residues in dark red), showing that, despite the hydrophobic nature of the NC, polar residues were significantly involved in the NC–protein interaction in the IvC and IvN simulations. Specifically, histidine (His) residues showed a relevant fraction of contacts with the NC in the IvC simulations (Fig. 3A). As shown in Fig. S8 in the SI, this was mainly due to the abundance of His residues at the C-terminal tail. Nonetheless, the relevant role of His in the interaction pattern of the IvC simulations appeared particularly interesting. As a matter of fact, histidine is a polar amino acid with an aromatic ring and was able to interact with the NC either by aromatic contacts (stacked/edge-to-face) between its aromatic ring and that of the phenyl groups that functionalize the NC or by interacting with exposed sulfur–Au atoms (Fig. 3A, inset). The latter interaction likely arose from the balance between the propensity of His nitrogen to interact with gold⁷¹ and the repulsion between the negatively charged His nitrogen and NC sulfur atoms. Another important amino acid in the interaction pattern of the IvC MDs was isoleucine (Ile). Ile is an apolar amino acid and its interaction with the NC was driven by purely hydrophobic effects: the aliphatic side chain of Ile interacted with the phenyl rings of the NC (Fig. 3A, inset).

Fig. 3 Percentage of contacts between the NC and a specific amino acid with respect to the total number of protein–NC contacts. A: IvC MDs; B: IvN MDs; and C: IvM MDs. Representative snapshots of the interaction between the NC and the two amino acids showing more contacts are reported below each histogram. The residues are arranged in accordance with the Kyte–Doolittle hydrophobicity scale.⁴⁹ The residues are arranged in a descending order of hydrophobicity, starting from the most hydrophobic, indicated by the colour blue, located on the left, to the most hydrophilic, indicated by the colour red, located on the right.

The interaction between the NC and a polar aromatic amino acid was also present in the IvN group (Fig. 3B), in which frequent interactions with tyrosine (Tyr) were observed. The Tyr phenol side chain was observed to interact with the phenyl rings of the NC via both stacked and edge-to-face aromatic interactions (Fig. 3B, inset). The residues showing the most frequent interaction with the NC in the IvN MDs were the residues of the N-terminal SS bridges (i.e., Cys-34 and Cys-36). Inspection of the trajectories revealed that this interaction mainly involved an exposed portion of the NC core and the backbone of the residues involved in the SS bridges (Fig. 3B, inset). Note that the preferential interaction of the NC with Tyr and Cys residues in the IvN simulations was mainly due to the fact that these residues were particularly well exposed: the N-termini of the two monomers are linked by two disulfide bridges, with two Tyr residues in close proximity (Fig. 3B, inset).

In the case of the IvM simulations, the NC was in contact with a part of the protein that is distinctly hydrophobic. The most frequent interactions were therefore with apolar residues such as Ile and phenylalanine (Phe) (Fig. 3C). The latter can effectively interact with the NC via aromatic stacking. Interestingly, all MD simulations showed that the surface functional groups adapted their conformation to optimize the interaction with the protein (see, e.g., the inset of Fig. 3B).

Next, we investigated the possible effects of the above-described interactions with the NC on the structural and dynamical properties of ChR2. To understand if there was any major change in the protein conformation, we monitored its secondary structure when interacting with the NC and compared it to that along the two MDs without the NC. The results, shown in Fig. S3–S5 of the SI, showed that the ChR2 secondary structure was rather stable also when interacting with the NC.

To investigate possible effects on ChR2 dynamics, we started by monitoring protein flexibility using RMSF analysis. To better highlight possible differences induced by the interaction with the NC, we compared the 17 MD simulations in the presence of the NC with the two NC-free simulations. In Fig. 4, we report, for each protein residue, the difference Δ_RMSF = (RMSF_withNC − RMSF_withoutNC). As illustrated in Fig. 4, an increased flexibility due to the interaction with the NC could be observed regardless of the interaction site. Specifically, the N-terminal, the intra- and extracellular loops ICL-1, ICL-2 and ECL-1, and helices TM-1, TM-5, and TM-6 showed an increased RMSF. Some of these regions were directly influenced by the interaction with the NC, while the increased fluctuation of other regions was attributable to allosteric effects, a typical phenomenon in ligand binding.⁷² The C-terminal region was affected in opposite ways, either an increased or a decreased fluctuation. However, stiffening of the C-terminal region that was in contact with the NC was observed in most IvC simulations. In contrast, an increased fluctuation of the N-terminal region in contact with the NC was present in the IvN simulations.

Fig. 4 Plot reporting Δ_RMSF = (RMSF_withNC − RMSF_withoutNC) for each residue in the IvC (black), IvN (red) and IvM (green) MDs. The black dashed lines represent distinct threshold values, which have been calculated as the sum between the average of Δ_RMSF plus one, two or three times the average of the standard deviation σ_RMSF. The circle marker denotes the residues of monomer 1, and the triangle marker denotes the residues of monomer 2. At the bottom, a representative protein snapshot highlights the regions corresponding to significant Δ_RMSF changes.

Since helices TM-1, TM-5, and TM-6 showed an increased fluctuation, we focused on the behavior of all seven transmembrane helices by computing the root-mean-square deviation (RMSD) on each ChR2 helix along all simulations. The analysis indicated that the protein transmembrane region was essentially stable, as evidenced by the relatively low values observed for all helices (≈1 Å). Nonetheless, helices TM-3 and TM-5 displayed higher RMSD values in some simulations (see Fig. S7 in the SI). Inspection of these simulations revealed that the high deviation of TM-3 could be attributed to the substantial impact of the NC on the ICL-2 region, which is closely associated with the TM-3 helix. The high deviation of TM-5 could instead be attributed to a kink of the helix region in contact with the NC. Therefore, the increase in RMSD could be attributed in both cases to a direct effect of the presence of the NC.

A relevant structural feature in rhodopsins is the amount of water in the protein interior in the dark state, which can affect the protonation state of crucial residues with a major impact on the photocycle evolution.⁵³ In principle, the interaction of the NC with either the C- or N-terminal may affect water accessibility within the protein interior. Specifically, the NC could obstruct water entry or induce conformational changes in the terminal regions, potentially leading to partial opening of the channel even in the dark state, where the channel is expected to remain closed. However, the analysis of the number of water molecules in the protein interior, shown in Fig. 5, indicated that ChR2 hydration did not vary significantly, despite the NC interacting with the region near the channel entrance.

Fig. 5 A: Average number of water molecules in the protein interior N_w in the IvC (black), IvN (red) and IvM (green) MD simulations. Error bars are obtained by the standard deviation of N_w along each MD. The blue dashed horizontal line corresponds to the average number of water molecules in the two simulations without the NC, while the blue dotted lines above and below indicate the standard deviation along the same simulations. B: Representative structure showing the protein interior hydration. Water molecules within 3.2 Å from the protein are represented as a red surface, and the positions of Leu110 and Leu138 are highlighted by yellow spheres.

The average number of water molecules in the protein interior was determined by calculating, for each MD frame, the number of water molecules N_w within 3.2 Å from the protein and whose z coordinate along the membrane normal was between that of Leu110 and that of Leu138 (see Fig. 5B). Importantly, N_w represents the average hydration per monomer, computed by dividing the total number of internal water molecules in the dimer by two. A monomer-specific analysis revealed in fact no statistically significant differences in hydration variability between the UP and DOWN monomers. In the simulations without the NC, N_w = 162 ± 8. The simulations with the NC yielded comparable mean values of N_w, suggesting that the interaction with the NC did not affect the channel structure in the dark-adapted state.

As mentioned in the Methods section, in a previous MD simulation we observed two possible hydrogen bond (HB) patterns for glutamate (Glu) 90 (see Fig. 6A): the UP pattern, in which Glu-90 interacts with asparagine (Asn) 258, and the DOWN pattern, in which Glu-90 interacts with aspartic acid (Asp) 253.⁵¹ We therefore monitored along the IvC, IvN and IvM simulations the HB between Glu-90 and Asn-258 side chains. In addition, we also assessed the stability of the interaction of the protonated Schiff base (PSB) and its two counterions by monitoring the presence of a HB between the NH group of PSB and the side chain of Glu-123 and Asp-253 (the three monitored HBs are highlighted in Fig. 6A). The applied HB criteria were based on a heavy atom (O⋯O and N⋯O) distance cutoff of 3.0 Å, a hydrogen-acceptor distance cutoff of 2.5 Å, and a minimum donor-hydrogen-acceptor angle of 145°, in line with recent high-resolution protein structural analyses.⁷³ It was previously observed that in ChR2, the PSB can form a direct HB with either of the two counterions, both direct and water-mediated.⁷⁴ Consistent with this observation, we observed in our simulations, a HB between the PSB and either Glu-123 or Asp-253. Our data also suggested that the presence of the NC promoted dynamic exchange between these two configurations. Fig. 6 also shows that the HB pattern of the Glu-90 was affected by the presence of the NC: an increase in the population of the UP Glu-90 conformation by ≈20% was observed in the IvC and IvN simulations. This increase was consistent in both the UP and DOWN configurations, thus indicating that interaction with the NC stabilized the UP conformation and promoted the transition from DOWN to UP. Conversely, in the IvM simulations, the Glu-90/Asn-258 interaction was weakened by ≈20%. Notably, it was shown that Glu-90 plays a major role in ChR2 photocycle branching.⁵³ The present observation of an NC-induced change in the Glu-90 HB pattern suggested, therefore, that the photocycle branching could be affected by the presence of the NC. Interestingly, a reduction in the branching probability for ChR2 in proximity to gold surfaces was previously suggested on the basis of experimental results.⁷⁵

Fig. 6 A: Representative snapshots of the UP Glu-90 conformation, in which Glu-90 is hydrogen bonded to Asn-258 (left) and the DOWN Glu-90 conformation, in which Glu-90 points toward the retinal cavity (right). In the two snapshots, residues Glu-90, Glu-123, Asp-253, Ans-258 and the PSB are shown as licorice and the following distances are highlighted: PSB/Glu-123, PSB/Asp-253, and Glu-90/Asn-258. For clarity, only the hydrogen atoms involved in these three HBs are shown. B: Fraction of frames in which the three HBs are present.

3.3 DPCN analysis

In order to conduct a more accurate investigation of the local and global perturbations induced by the interaction with the NC, a DPCN analysis was performed. This method was chosen because it allowed for a detailed examination of the variations in protein contact networks due to NC interaction. The DPCN analysis, combined with the study of its CCs, i.e. subgraphs of nodes with highly correlated contact variations, provided a detailed map of localized perturbation effects. The limited size of the CCs confirmed that the NC influence was predominantly localized, with restricted perturbations found across the whole system but without altering the overall stability of the protein. This analysis identified the key regions most responsive to the perturbations induced by the NC binding (see Fig. 7): the direct NC interaction site at the N- or C-terminal, the retinal cavity, and the TM1–TM2 helical region.

Fig. 7 (A) Representation of the high-betweenness centrality nodes (orange spheres) identified across all interaction groups (IvC, IvN, and IvM). These 22 conserved residues form a continuous allosteric pathway connecting key functional regions: the N-terminal domain (black dotted box), the transmembrane helices 1 (TM1, red box) and 2 (TM2, blue box), and the C-terminal region (grey box). Key functional residues, including the gating-related Glu-90 and His-134, and the retinal binding pocket (brown dashed oval), are highlighted. The light red surface represents the internal water network (Quicksurface representation). (B and C) Representative snapshots of the DPCN analysis for the IvC (B) and IvN (C) interaction groups. The ChR2 dimer is shown with the nanocluster (yellow/cyan spheres) interacting at the C-terminus and N-terminus, respectively. The colored boxes highlight the localized connected components (CCs) that cluster within the functional regions identified in panel A.

Then, we investigated which specific residues were responsible for mediating and propagating these signals. To this end, we computed the betweenness centrality on the DPCN graph. Since edge weights in the DPCN represent contact perturbations, we defined the edge ‘distance’ for the shortest-path calculation as the inverse of the weight d_ij = 1/|w_ij|. This transformation ensures that the computed shortest paths, and consequently the betweenness centrality, trace the shortest routes of signal propagation. Using this metric, we identified the minimal set of nodes accounting for 50% of the total betweenness in each network. This approach enabled a direct comparison of the dominant information flow mediators across different interaction groups (IvC, IvN, and IvM), highlighting both common communication nodes and condition-specific differences in signal propagation pathways (see Fig. 7). This centrality analysis revealed a core set of 22 residues that functioned as high-betweenness nodes across all three interaction groups (Fig. S9). Spatially, these nodes formed a continuous pathway that connected the C-terminal, intracellular gate (ICG), central gate (CG), and N-terminal domain, suggesting the presence of an allosteric communication network. A comparison with existing literature showed that this set of nodes includes several residues with well-documented functional importance. For instance, Glu-82, Glu-90, and Asn-258 are known to influence the local electrostatics and affect ion permeation within the ICG/CG region.⁷⁶ Moreover, residue His-134 has been demonstrated to modulate the channel's photocurrent and closing time.⁷⁷ A notable concentration of eight nodes was also found on the TM2 helix, a region that is recognized for its role in channel gating.⁷⁸

Finally, the DPCN analysis revealed that the interaction with the NC induced localized structural perturbations within the retinal binding pocket. Specifically, CCs involving the retinal chromophore were detected in both monomers, regardless of which monomer directly contacted the nanocluster or where the interaction occurred (see Fig. 7). This finding indicates that the structural perturbation propagates allosterically to the retinal cavity in both monomers. Given the well-established sensitivity of retinal photophysics to its local environment,^79–83 such perturbations could potentially alter ChR2's function by modulating the retinal's photochemical cycle.

The DPCN analysis showed that among all the MD simulations performed, IvM₁ was the only trajectory in which the protein network was markedly affected by the interaction with the NC. In this simulation, a large CC spanning almost the entire monomer interacting with the NC was observed (Fig. 8A), in clear contrast to the localized CCs described above for the other trajectories. Notably, IvM₁ was also the only trajectory in which a pronounced kink of the helix was observed. To better characterize the dynamics of the IvM₁ simulation, we performed a principal component analysis (PCA).⁸⁴ We computed the eigenvectors of the protein's α-carbons covariance matrix on a trajectory obtained by concatenating one trajectory without the NC and the IvM₁ MD simulation. Subsequently, the two trajectories were projected separately on the plane of the first two eigenvectors obtained from the concatenated trajectory. Such a projection is shown in Fig. 8B (the trajectory without the NC is colored in black while the IvM₁ trajectory is colored according to time). The figure shows that the conformational basin without the NC remained localized in a relatively limited region of the plane while the IvM₁ conformational basin was rather different. As shown by the representative structures extracted from the IvM₁ basins (panel C), this was mainly due to the kink of TM-5. However, it is interesting to observe that as the trajectory evolved in time (panel B), the difference in the projection on eigenvector 1 with respect to the MD without the NC was partially recovered. This observation is in accordance with the time evolution of the RMSD of TM-5 (Fig. S7D in the SI), which shows that the helix kink only occurred in the first part of the trajectory. This suggests that the protein is capable of recovering its initial conformation and shows that ChR2 is rather stable and can adapt to the external perturbation exerted by the NC without undergoing major structural or dynamical changes.

Fig. 8 Panel A: connected components (CCs) identified via DPCN analysis for the IvM₁ simulation of ChR2 (cyan cartoon) with the NC (spheres, top left). A large connected component (CC₁, blue), spanning the entire interacting monomer, can be observed. Other significant CCs (CC₂ red, CC₃ green, and CC₄ yellow) are also depicted. Panel B: projection of the trajectories on the plane of the first two eigenvectors of the C_α covariance matrix for the trajectories without the NC (black) and IvM₁ (colored according to time). The red crosses mark densely populated regions from which representative structures are extracted. The representative structures, Panel C, are shown as superimposed, and the protein regions with more evident conformational differences are colored according to the time frame.

4 Conclusions

Our simulations reveal that functionalized gold nanoclusters can stably engage ChR2 without compromising its structural integrity or native dynamic ensemble, despite frequent nanoscale surface contact and membrane insertion events. Across 31 molecular dynamics trajectories, spontaneous nanocluster–protein complex formation occurs in roughly half of the configurations, primarily stabilized by aromatic-rich interfacial interactions. In other cases the nanocluster partitions into the lipid bilayer, suggesting alternative access routes to transmembrane interaction via lateral diffusion. Remarkably, ChR2 retains its global fold, water content, and retinal Schiff-base configuration, demonstrating an unexpected capacity to withstand nanoscale perturbations and recover native-like conformational states even after local distortion. Beyond structural resilience, we identify a specific mechanistic effect of nanocluster engagement: a rewiring of the Glu-90 hydrogen-bond network, an essential determinant of ChR2 photocycle branching.

Network analyses based on dynamical perturbation contact networks further show that nanocluster engagement primarily induces localized structural perturbations within ChR2, consistently involving the same set of protein regions across all simulations. These regions include the terminal domains and transmembrane helices TM-1 and TM-2. Notably, regardless of whether a stable nanocluster–protein complex forms, we consistently observe a perturbation of the retinal cavity, suggesting sensitivity of this region to external nanoscale interactions and a potential impact on the chromophore's photophysical behavior. In particular, the IvM simulations in which the nanocluster interacts laterally with the protein and approaches the retinal cavity display contact geometries that would be especially conducive to plasmon-mediated modulation of the retinal chromophore for similarly sized nanoclusters with plasmonic character. These observations suggest that metallic NCs may affect the photophysics of ChR2 in two ways: (i) indirectly, by altering the environment surrounding the chromophore, and (ii) directly, by exerting an effect on the optical properties of retinal. The latter way is unlikely for the Au₂₅ cluster investigated here, because of its molecular-like electronic excitations, but could be important for other nanoclusters in the same spatial size range (1–2 nm, such as Au₁₄₄SR₆₀ ⁴⁶ or also small silver clusters such as Ag₂₀ ⁴⁴) that exhibit plasmonic-like character.

Collectively, these results establish a molecular feasibility framework for structurally non-disruptive integration of metallic NCs with membrane photoreceptors. By providing mechanistic insights into the interaction between metallic NCs and ChR2, this work opens a conceptual path toward hybrid nano-optogenetic architectures, where plasmon-supported near-fields could re-engineer chromophore photophysics and neuronal activation dynamics without genetic or chemical modification of the receptor. More broadly, our findings may contribute to the rational design of plasmonic nanostructures for interfacing with functional membrane proteins, advancing foundational capabilities in nano-neurotechnology and protein–nanomaterial hybrid systems.

Author contributions

RM: data curation, formal analysis, investigation, software, and writing – original draft. LB: investigation, supervision, and writing – review & editing. SC: conceptualization, supervision, and writing – review & editing. LC: data curation and software. IR: conceptualization, methodology, supervision, and writing – review & editing. SD: conceptualization, funding acquisition, investigation, supervision, writing – original draft, and writing – review & editing. LZP: conceptualization, investigation, supervision, writing – original draft, and writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: schematics of the 26 MD simulation starting structures; contact maps between ChR2 and the NC for the two IvC group simulations not shown in Figure 2A; time evolution of the ChR2 secondary structure; color map of Δ_RMSF; RMSD of transmembrane helices TM-3 and TM-5; spatial distribution of histidine residues within ChR2; Venn diagram illustrating the overlap of hub residues among the IvC, IvN, and IvM groups. See DOI: https://doi.org/10.1039/d6nr00330c.

The code for the DPCN analysis can be found at https://github.com/agheeraert/pmdlearn.

Acknowledgements

This study was conducted using the DARIAH HPC-AI cluster at CNR-NANOTEC in Lecce, funded by the “MUR PON Ricerca e Innovazione 2014-2020” project, code PIR01_00022 and H2IOSC Project – Humanities and Cultural Heritage Italian Open Science Cloud funded by the European Union – NextGenerationEU – NRRP M4C2 – Project code IR0000029. R. M. and S. D thank the cluster manager of the DARIAH cluster, Dr Gian Paolo Marra, for the precious computational support. L. Z. P. acknowledges funding from MUR through PRIN 2022-WS44W4-EnvELOP (CUP B53D23014000006).

References

1. M. Carter, R. Essner, N. Goldstein and M. Iyer, Guide to Research Techniques in Neuroscience (Third Edition), Academic Press, 3rd edn, 2022, pp. 191–208.
2. M. Parisa, V. Hadi and G. P. Firooz, J. Med. Signals Sens., 2017, 7, 71–79.
3. É. Pastrana, Nat. Methods, 2011, 8, 24–25.
4. J.-A. Sahel, E. Boulanger-Scemama, C. Pagot, A. Arleo, F. Galluppi, J. N. Martel, S. D. Esposti, A. Delaux, J.-B. de Saint Aubert and C. de Montleau, et al., Nat. Med., 2021, 27, 1223–1229.
5. A. Spreen, D. Alkhoury, H. Walter and S. Müller, Iscience, 2024, 27, 109776.
6. F. Zhang, J. Vierock, O. Yizhar, L. E. Fenno, S. Tsunoda, A. Kianianmomeni, M. Prigge, A. Berndt, J. Cushman and J. Polle, et al., Cell, 2011, 147, 1446–1457.
7. O. P. Ernst, D. T. Lodowski, M. Elstner, P. Hegemann, L. S. Brown and H. Kandori, Chem. Rev., 2014, 114, 126–163.
8. A. Berndt, S. Y. Lee, J. Wietek, C. Ramakrishnan, E. E. Steinberg, A. J. Rashid, H. Kim, S. Park, A. Santoro and P. W. Frankland, et al., Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 822–829.
9. K. Deisseroth and P. Hegemann, Science, 2017, 357, eaan5544.
10. S. Gozem, H. L. Luk, I. Schapiro and M. Olivucci, Chem. Rev., 2017, 117, 13502–13565.
11. H. Kandori, Biophys. Rev., 2020, 12, 355–361.
12. M. N. Sardoiwala, A. K. Srivastava, S. Karmakar and S. R. Choudhury, ACS Chem. Neurosci., 2019, 10, 3375–3385.
13. F. Pisano, L. Collard, D. Zheng, M. F. Kashif, M. Kazemzadeh, A. Balena, L. Piscopo, M. S. Andriani, M. De Vittorio and F. Pisanello, Neurophotonics, 2024, 11, S11513.
14. A. Biesso, W. Qian and M. A. El-Sayed, J. Am. Chem. Soc., 2008, 130, 3258–3259.
15. C.-W. Yen, L.-K. Chu and M. A. El-Sayed, J. Am. Chem. Soc., 2010, 132, 7250–7251.
16. T. P. Rossi, T. Shegai, P. Erhart and T. J. Antosiewicz, Nat. Commun., 2019, 10, 3336.
17. G. Giannone, F. D. Sala and S. D'Agostino, Nanomater. Nanotechnol., 2019, 9, 1847980419856534.
18. L. Cascino, S. Corni and S. D'Agostino, J. Phys. Chem. C, 2024, 128, 2917–2927.
19. R. F. Ribeiro, L. A. Martínez-Martínez, M. Du, J. Campos-Gonzalez-Angulo and J. Yuen-Zhou, Chem. Sci., 2018, 9, 6325–6339.
20. T. E. Li, B. Cui, J. E. Subotnik and A. Nitzan, Annu. Rev. Phys. Chem., 2022, 73, 43–71.
21. M. Ruggenthaler, D. Sidler and A. Rubio, Chem. Rev., 2023, 123, 11191–11229.
22. J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi and A. Rubio, Nature, 2022, 606, 41–48.
23. J. Fregoni, F. J. Garcia-Vidal and J. Feist, ACS Photonics, 2022, 9, 1096–1107.
24. A. Mandal, M. A. Taylor, B. M. Weight, E. R. Koessler, X. Li and P. Huo, Chem. Rev., 2023, 123, 9786–9879.
25. M. Sánchez-Barquilla, A. I. Fernández-Domínguez, J. Feist and F. J. García-Vidal, ACS Photonics, 2022, 9, 1830–1841.
26. D. S. Wang and S. F. Yelin, ACS Photonics, 2022, 8, 2818–2826.
27. F. Herrera and J. Owrutsky, J. Chem. Phys., 2020, 152, 100902.
28. F. Herrera and F. C. Spano, ACS Photonics, 2017, 5, 65–79.
29. G. Li and R. Jin, Acc. Chem. Res., 2013, 46, 1749–1758.
30. Y. Zhu, H. Qian, M. Zhu and R. Jin, Adv. Mater., 2010, 22, 1915–1920.
31. X. Nie, H. Qian, Q. Ge, H. Xu and R. Jin, ACS Nano, 2012, 6, 6014–6022.
32. A. A. Sousa, S. A. Hassan, L. L. Knittel, A. Balbo, M. A. Aronova, P. H. Brown, P. Schuck and R. D. Leapman, Nanoscale, 2016, 8, 6577–6588.
33. Y. Liu, X. Lv, H. Liu, Z. Zhou, J. Huang, S. Lei, S. Cai, Z. Chen, Y. Guo, Z. Chen, X. Zhou and L. Nie, Nanoscale, 2018, 10, 3631–3638.
34. F. Ghahremani, A. Kefayat, D. Shahbazi-Gahrouei, H. Motaghi, M. A. Mehrgardi and S. Haghjooy-Javanmard, Nanomedicine, 2018, 13, 2563–2578.
35. X.-D. Zhang, J. Chen, Z. Luo, D. Wu, X. Shen, S.-S. Song, Y.-M. Sun, P.-X. Liu, J. Zhao, S. Huo, S. Fan, F. Fan, X.-J. Liang and J. Xie, Adv. Healthcare Mater., 2014, 3, 133–141.
36. V. P. Zharov, E. N. Galitovskaya, C. Johnson and T. Kelly, Lasers Surg. Med., 2005, 37, 219–226.
37. R. Antoine, D. Maysinger, L. Sancey and V. Bonačić-Koutecký, Commun. Chem., 2022, 5, 47.
38. N. Goswami, K. Zheng and J. Xie, Nanoscale, 2014, 6, 13328–13347.
39. J. Xie, Y. Zheng and J. Y. Ying, J. Am. Chem. Soc., 2009, 131, 888–889.
40. A. Soleilhac, F. Bertorelle and R. Antoine, Spectrochim. Acta, Part A, 2018, 193, 283–288.
41. M. W. Heaven, A. Dass, P. S. White, K. M. Holt and R. W. Murray, J. Am. Chem. Soc., 2008, 130, 3754–3755.
42. G. Brancolini, D. Toroz and S. Corni, Nanoscale, 2014, 6, 7903.
43. S. D'Agostino, R. Rinaldi, G. Cuniberti and D. S. Fabio, J. Phys. Chem. C, 2018, 122, 19756–19766.
44. R. Schira and F. Rabilloud, J. Phys. Chem. C, 2019, 123, 6205–6212.
45. S. Malola, L. Lehtovaara, J. Enkovaara and H. Häkkinen, ACS Nano, 2013, 7, 10263–10270.
46. C. Yi, M. A. Tofanelli, C. J. Ackerson and K. L. Knappenberger Jr, J. Am. Chem. Soc., 2013, 135, 18222–18228.
47. M. R. Rasch, E. Rossinyol, J. L. Hueso, B. W. Goodfellow, J. Arbiol and B. A. Korgel, Nano Lett., 2010, 10, 3733–3739.
48. J. D. Marshall and M. J. Schnitzer, ACS Nano, 2013, 7, 4601–4609.
49. J. Kyte and R. F. Doolittle, J. Mol. Biol., 1982, 157, 105–132.
50. O. Volkov, K. Kovalev, V. Polovinkin, V. Borshchevskiy, C. Bamann, R. Astashkin, E. Marin, A. Popov, T. Balandin, D. Willbold, G. Büldt, E. Bamberg and V. Gordeliy, Science, 2017, 358, eaan8862.
51. L. Bellucci, M. Capone, I. Daidone and L. Zanetti-Polzi, Int. J. Biol. Macromol., 2025, 305, 140977.
52. W. Zhang, T. Yang, S. Zhou, J. Cheng, S. Yuan, G. V. Lo and Y. Dou, Biomolecules, 2019, 9, 852.
53. J. Kuhne, J. Vierock, S. A. Tennigkeit, M.-A. Dreier, J. Wietek, D. Petersen, K. Gavriljuk, S. F. El-Mashtoly, P. Hegemann and K. Gerwert, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 9380–9389.
54. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
55. R. B. Best and G. Hummer, J. Phys. Chem. B, 2009, 113, 9004–9015.
56. K. Lindorff-Larsen, S. Piana, K. Palmo, P. Maragakis, J. L. Klepeis, R. O. Dror and D. E. Shaw, Proteins: Struct., Funct., Bioinf., 2010, 78, 1950–1958.
57. S. Hayashi, E. Tajkhorshid and K. Schulten, Biophys. J., 2002, 83, 1281–1297.
58. E. Malmerberg, Z. Omran, J. S. Hub, X. Li, G. Katona, S. Westenhoff, L. C. Johansson, M. Andersson, M. Cammarata and M. Wulff, et al., Biophys. J., 2011, 101, 1345–1353.
59. J. P. Jambeck and A. P. Lyubartsev, J. Phys. Chem. B, 2012, 116, 3164–3179.
60. J. P. Jambeck and A. P. Lyubartsev, J. Chem. Theory Comput., 2012, 8, 2938–2948.
61. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
62. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
63. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
64. B. Hess, H. Bekker, H. J. Berendsen and J. G. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
65. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
66. A. Gheeraert, L. Pacini, V. S. Batista, L. Vuillon, C. Lesieur and I. Rivalta, J. Phys. Chem. B, 2019, 123, 3452–3461.
67. A. Gheeraert, C. Lesieur, V. S. Batista, L. Vuillon and I. Rivalta, J. Phys. Chem. B, 2023, 127, 7571–7580.
68. J. Salamanca Viloria, M. F. Allega, M. Lambrughi and E. Papaleo, Sci. Rep., 2017, 7, 2838.
69. M. Sobieraj and P. Setny, Proteins:Struct., Funct., Bioinf., 2021, 89, 1333–1339.
70. G. Aria, pmdlearn: Protein Molecular Dynamics LEARNing, https://github.com/agheeraert/pmdlearn.git, 2023, Accessed: 2025-03-08.
71. F. Iori, S. Corni and R. Di Felice, J. Phys. Chem. C, 2008, 112, 13540–13545.
72. S. J. Wodak, E. Paci, N. V. Dokholyan, I. N. Berezovsky, A. Horovitz, J. Li, V. J. Hilser, I. Bahar, J. Karanicolas and G. Stock, et al., Structure, 2019, 27, 566–578.
73. M. Tsujimura, H. Ishikita and K. Saito, Phys. Chem. Chem. Phys., 2025, 27, 9794–9805.
74. Y. Guo, F. E. Wolff, I. Schapiro, M. Elstner and M. Marazzi, Phys. Chem. Chem. Phys., 2018, 20, 27501–27509.
75. R. Polito, M. E. Temperini, E. Ritter, L. Puskar, U. Schade, M. Broser, P. Hegemann, L. Baldassarre, M. Ortolani and V. Giliberti, Phys. Rev. Appl., 2021, 16, 014048.
76. T. Yang, W. Zhang, J. Cheng, Y. Nie, Q. Xin, S. Yuan and Y. Dou, Int. J. Mol. Sci., 2019, 20, 3780.
77. K. Feldbauer, D. Zimmermann, V. Pintschovius, J. Spitz, C. Bamann and E. Bamberg, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 12317–12322.
78. M. Müller, C. Bamann, E. Bamberg and W. Kühlbrandt, J. Mol. Biol., 2015, 427, 341–349.
79. M. Braiman and R. Mathies, Proc. Natl. Acad. Sci. U. S. A., 1982, 79, 403–407.
80. G. Bassolino, T. Sovdat, A. Soares Duarte, J. M. Lim, C. Schnedermann, M. Liebel, B. Odell, T. D. Claridge, S. P. Fletcher and P. Kukura, J. Am. Chem. Soc., 2015, 137, 12434–12437.
81. P. P. Roy, Y. Kato, R. Abe-Yoshizumi, E. Pieri, N. Ferré, H. Kandori and T. Buckup, Phys. Chem. Chem. Phys., 2018, 20, 30159–30173.
82. P. P. Roy, R. Abe-Yoshizumi, H. Kandori and T. Buckup, J. Phys. Chem. Lett., 2019, 10, 1012–1017.
83. P. Malakar, S. Gholami, M. Aarabi, I. Rivalta, M. Sheves, M. Garavelli and S. Ruhman, Nat. Commun., 2024, 15, 2136.
84. I. Daidone and A. Amadei, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2012, 2, 762–770.

---