Computational development of a phase-sensitive membrane raft probe

Derivatives of the widely used 1,6-diphenyl-1,3,5-hexatriene have been considered using a multiscale approach involving spin-flip time-dependent density functional theory, classical molecular dynamics and hybrid quantum mechanics / molecular mechanics. We identify a potential probe of membrane phase (i.e. to preferentially detect liquid-ordered regions of lipid bilayers), which exhibits restricted access to a conical intersection in the liquid-ordered phase but is freely accessible in less ordered molecular environments. The characteristics of this probe also mark it as a candidate for an aggregation induced emission fluorophore. Page 1 of 22 Physical Chemistry Chemical Physics P hy si ca lC he m is tr y C he m ic al P hy si cs A cc ep te d M an us cr ip t O pe n A cc es s A rt ic le . P ub lis he d on 1 8 M ar ch 2 02 2. D ow nl oa de d on 3 /2 1/ 20 22 9 :3 7: 33 A M . T hi s ar tic le is li ce ns ed u nd er a C re at iv e C om m on s A ttr ib ut io n 3. 0 U np or te d L ic en ce . View Article Online DOI: 10.1039/D2CP00431C


Introduction
Lipid membranes are heterogeneous structures that not only define the boundaries of cells, but are also responsible for regulating many key processes in living cells. In regions of membranes with high concentrations of sphingolipids and cholesterol, liquid-ordered (L o ) phases can occur, often known as lipid rafts. These can be characterised by thickening of the bilayer and relatively small tilt angles for cholesterol, 1,2 and having sizes between 10 and 200 nm. 3 These ordered regions are critical to many of the processes regulated by the cell membrane, including signalling, exo-and endocytosis and protein trafficking. 4,5 Small molecules can be directly trafficked through a membrane without a specific protein. In a recent study by Ghysels et al., 6 trafficking of oxygen and water molecules through membranes is made possible by transiting along the boundary regions of the different membrane phases. Direct observation of rafts in vivo is still elusive, leading to results that can seem at odds with each other. 7 Initially, extraction of raft regions was performed with detergent assays, although differing temperatures and detergent concentrations have led to very different molecular compositions of the raft domains, even from similar initial conditions. [8][9][10] While fluorescent probes which are sensitive to membrane environments have generally been very successful on reporting within general lipid bilayers, there has been less success in fluorescent labelling of specific membrane phases. 11 A simple probe of L o phase existence would therefore prove invaluable to understanding the existence and effects of raft nanodomains in a cell.
During the last decade, a new class of membrane probes that report on lateral forces present in lipid membranes have been developed by the Matile group. 12-22 These probes, based on a dithienothiophene molecular structure, exhibit large shifts in their excitation spectra upon partitioning into regions of highly-ordered lipid domains. A computational approach identified the planarization of the structure as the key element in the large spectral shifts observed between the probe in L o and liquid-disordered (L d ) regions of lipid membranes. 14 A single dihedral angle linking the two thienothiophene moieties exhibits mechano-sensitivity and rotates in response to forces exerted by the molecular environment, leading to a red-shift of the absorption spectrum and increased intensity of the emission spectrum upon ordering of the surrounding membrane environment. 14 We have previously studied contemporary and novel fluorescent probes of membrane structure and function, including di-8-ANEPPS and BODIPY derivatives commonly encountered in membrane characterisation. [23][24][25][26][27] In the current work, we explore a series of molecules inspired by the design principles of molecules which exhibit aggregation-induced emission (AIE) effect in their fluorescent behaviour. In AIE, aggregation of specially designed fluorophores leads to enhanced emission, often due to restriction of rotation at key points within molecules, leading to restriction of access to non-radiative decay pathways (e.g. conical intersections). We investigate, from first principles calculations, simple derivatives of the highly used 1,6-diphenyl-1,3,5-hexatriene (DPH) probe. 28,29 DPH itself doesn't show much sensitivity with respect to fluorescence wavelength and intensity between different phases of a lipid bilayer, although it does demonstrate fluorescence anisotropy differences between different phases. 28,29 Inspired by the simple substitutions made to 1,4distyrylbenzene which significantly enhance the AIE effect, 30 we demonstrate that inclusion of methyl groups along the polyene chain of DPH disrupts the planar nature of the molecule in gas-phase and solvated environments, while the forces involved within a raft-like membrane planarise the DPH derivatives resulting in emission being "switched on" within the raft domains due to a restricted availability to access the conical intersection. 31

Computational Details
Our overall approach to identifying and confirming a fluorescent probe sensitive to phase changes in a membrane are summarised in Scheme 1. In essence, we start with gas phase calculations to identify probes that have a conical intersection in which the minimum energy conical intersection (MECI) geometry has a twist in a polyene backbone. We then look to see if the conical intersection is thermally accessible in the S 1 electronic state. Once we have identified any such probe(s), we then perform classical molecular dynamics (MD) simulations to equilibrate the lipid bilayer and probe molecule, followed by hybrid quantum mechanics / molecular mechanics (QM/MM) simulations to ensure that the conical intersection becomes less accessible as the ordering of the lipid bilayer increases. Scheme 1. Multi-layer approach to identifying potential candidates for phase-sensitive lipid membrane probes under standard temperature and pressure (STP) conditions.

Gas phase
We investigated the DPH derivatives shown in Figure 1, along with unsubstituted DPH. The structures for the S 0 electronic states of each molecule were optimised using density functional theory (DFT) with the BHHLYP functional 32,33 and 6-31G(d,p) basis set. The S 1 geometries were optimised using time-dependent DFT (TD-DFT) using the same functional and basis set. MECI geometries were optimised using spin-flip TDDFT (SF-TDDFT) and analytical non-adiabatic coupling derivatives, 34,35 as we have previously found this approach (along with BHHLYP) to be in good agreement with XMS-CASPT2 results. 36 The starting guess for the MECI search was generated from the S 1 optimised geometry, where the planarity of the polyene backbone was disrupted by moving the carbon atom at which methyl substitution(s) occurs 0.1 Å out of the plane of the conjugated polyene.
Ab initio molecular dynamics (AIMD) simulations were performed for each of the molecules, using the same method and basis set combinations above, and were propagated on the S 1 potential energy surface for 5,000 steps with a timestep of 0.484 fs. Initial velocities were taken from a random sampling of a Maxwell-Boltzmann distribution. A Fock matrix extrapolation scheme was employed, in which the previous 12 Fock matrices are saved and extrapolated using a sixth-order polynomial to give a guess for the Fock matrix at the current iteration. 37

Classical molecular dynamics
Force-field parameters for the DPH-derivatives in their S 0 electronic state were generated using the CHARMM General Force-Field (CGenFF). 39,40 In this approach, a penalty is given to each parameter as a measure of the confidence of the quality of the parameter. Scores below 10 indicate good analogy of the CGenFF assigned parameter, while scores above 10 require further parameterisation. 40 The majority of the new parameters with penalty scores above 10 were dihedral terms. We re-optimised each of the dihedral terms with penalty scores greater than 10 using a force-matching algorithm, 41 which we briefly outline here.
The force-field is optimised by considering the set of parameters, [p], with CGenFF scores greater than 10 and minimising an objective function, : where N is the number of atoms in the molecule, M is the number of reference structures (in this case, 15 reference structures sampled at equal timesteps from the AIMD simulations were used), and f ij are the forces of atom i from reference structure j from DFT or MD. A simple gradient descent approach is then used to calculate the descent direction of the parameter set. The gradient of the objective function, G i , for each parameter, p i , is obtained numerically: where α is a scaling parameter to determine the step-size of the descent. A total of 500 reference structures taken from the AIMD simulations were used in the parameterisation scheme given above.
The optimised force-fields were used to equilibrate the molecular probe and surrounding lipid membrane system using molecular dynamics as follows: the CHARMM-GUI web service 42 was used, along with the optimised CGenFF force-field, to create a membrane system comprising a single probe molecule and 200 lipids per leaflet (total of 400 sphingomyelin (d18:1/18:0) lipids for the "non-raft" system; 30 molar % of cholesterol was used for the "raft" system, as liquid-ordered phases are known to form at this concentration 43 , giving 280 sphingomyelin (d18:1/18:0) and 120 cholesterol molecules, distributed equally between the leaflets. We abbreviate sphingomyelin as SSM from herein). Water molecules were added to a minimum thickness of 22.5 Å above each leaflet of the bilayer, giving a total of 12,500-12,900 water molecules for the different systems considered. Potassium and chloride ions were added at a concentration of 0.15 M. The system was minimised for 10,000 steps and equilibrated for 1,875 ps over six steps using the CHARMM-GUI Membrane Builder, [44][45][46][47] which slowly removes restraints (see the Supporting Information for full details of the applied restraints at each step of equilibration). Production dynamics proceeded for 100 ns with a 2 fs timestep to give a fully equilibrated system. Simulations were performed at three different temperatures: 310 K, 320 K and 330 K, all using the NPT ensemble. Constant pressure was maintained using the Nosé-Hoover Langevin-piston algorithm 48,49 and constant temperature was achieved through Langevin dynamics. Long-range electrostatics were described using the Particle Mesh Ewald method, 50 using 6 th -order interpolation. Lennard-Jones interactions were used to describe van der Waals' interactions, with a force-switching function acting in the range of 10-12 Å. 51 The CHARMM36 lipid force-field 52-54 and TIP3P water model 55,56 were used.
Electron density profiles were constructed by calculating electron densities with a slab thickness of 0.8 Å after recentring the bilayer to place the leaflet interface at Z = 0 Å.
Deuterium order parameters were calculated using the equation where θ is the angle between a given C-H vector and the membrane normal.

Gas phase
The optimised geometries for the S 0 and S 1 states, along with the S 1 /S 0 MECIs, are given in the Supporting Information, while the MECI geometries are also shown in Figure 2. The MECI geometries all share a common motif: a significant distortion from planarity close to one of the methyl groups. This is to be expected, given that the initial guess geometries were distorted at this position. In all cases, the conical intersection exhibits a peaked topology. For an isolated molecule to exhibit phase-dependent quenching behaviour, and thus potentially membrane phase-dependent emission, fast non-radiative decay pathways must be easily accessible to prevent fluorescence in the gas-phase and in solution. Critical dihedral angles for the MECI geometries from the AIMD simulations are given in Tables S1-S10 in the Supporting Information. With the exception of 2Me, none of the molecules has an average dihedral angle close to those required to access the conical intersection from the AIMD simulations, and we posit that these molecules are unlikely to exhibit environmentsensitive behaviour based on this. However, 2Me has average angles (and large standard deviations) that indicate rotation around the bonds is not energetically costly and thus the molecule is able to freely rotate towards the MECI in the S 1 state (Table 1). In contrast, DPH would be unlikely to access the conical intersection in the gas phase, and would be expected to remain fluorescent, consistent with experimental findings. 28,29 Given this, 2Me shows promise as a candidate for environment-sensitive emission behaviour. Potential energy scans of 2Me in the S 0 and S 1 electronic states are given in Figures S2 and   S3, respectively. In the ground state, the energy minimum is found with the C1-C6-C7-C8 dihedral angle at ~35°, while in the S 1 excited state, the minimum is close to zero, indicating that emission occurs preferentially from a planar geometry (see Figure S1 for atom numbering for dihedral angles). In Figure S4 is the emission energy as a function of the C1-C6-C7-C8 dihedral angle taken from the S 1 potential energy scan. The emission energy is at a maximum at ~90° and at a minimum (~3.4 eV) at ~0°, demonstrating a 0.5 eV shift between the two extremes. In Figure S5, we consider the oscillator strengths for the S 1  S 0 transition along the S 1 potential energy scan (i.e. as an approximation to emission intensity).
Where emission is the only process occurring for de-excitation, this would be a good indication of the emission intensity. In our case, we know that at a dihedral angle of ~47°, the MECI is a highly competitive process to emission, which would reduce emission intensity further still. Given this, we would expect emission to shift to higher energies and lower intensities as the C1-C6-C7-C8 dihedral angle increases. From our AIMD results in the gas phase, we would expect non-radiative decay via the MECI and hence little to no emission. In the event of an AIE-like process in a lipid membrane, the properties shown here suggest that, if 2Me is responsive to the phase change from non-raft to raft regions, then we would expect either emission energy to significantly shift (and intensity change) upon detecting a phase change (see Figure S6), or emission to switch on only when the raft-like environment causes the C1-C6-C7-C8 dihedral angle to remain near planar.
A scan of the C1-C6-C7-C8 potential energy surface for both the lowest two singlet states is shown in Figure 3. While only an approximation to the exact g and h vectors of the conical intersection seam, this dihedral is very important in accessing the MECI and has a very wide range of angles at which the S 0 and S 1 states are near-degenerate. This suggests that, at angles greater than ~20°, the MECI becomes accessible for 2Me.

Classical molecular dynamics
We 2Me, in each of the simulations. The bilayer thickness in each simulation is given in Table 2.
As expected, the bilayers containing cholesterol are thicker than the SSM-only bilayers.   (Figures S7 and S8).

QM/MM simulations
Given in Table 3  is below that at which the conical intersection is accessible. The average value in the nonraft environment is 25.6° at the same temperature, in the region in which the conical intersection becomes accessible (see Figure 3). If the conical intersection is accessed, this would then result in non-radiative decay and hence the fluorescence switches off. In the event that fluorescence still occurs, the emission shift (and change in intensity) with respect to changes in the C1-C6-C7-C8 dihedral are large (see Figure S6) resulting in substantial reduction in intensity with increasing dihedral angle, thus, the probe will consistently report on the membrane environment. Raft (S 0 ) 5.4 (± 4.1) 6.3 (± 4.7) 6.9 (± 5.0) Raft (S 1 ) 6.2 (± 4.7) 7.3 (± 5.5) 7.7 (± 5.7) Table 5. C8-C9-C10-C11 average dihedral angles for DPH taken from the QM/MM simulations. Standard deviations are given in parentheses.

Conclusions
We have considered several modifications to the widely used membrane probe, DPH, using a computational protocol. From the gas phase calculations, we determined that the 2-methyl derivative (2Me) possessed the correct characteristics to act in an AIE-like fashion, in which access to the conical intersection between the S 0 and S 1 states is easily accessible as the molecule is able to freely rotate, but which shows restricted access to the conical intersection in the L o phase of the SSM and cholesterol lipid bilayer considered. This is due to the forces imposed by the ordered molecular environment of the raft-like membrane, similar to the operation of the dithienothiophene probes. [12][13][14][15][16][17][18][19][20][21][22] Access to the conical intersection is possible in the L d lipid bilayer (containing only SSM), and as such we conclude that 2Me is capable of detecting a phase change in a lipid membrane. We predict that emission will be strong in a raft environment, and either cease within a non-raft region, or substantially shift to longer wavelengths (and have a reduction in intensity), while exhibiting similar fluorescence anisotropy properties to DPH. These properties demonstrate the potential benefits of a simple modification to the widely used membrane probe, DPH, in detecting raft regions of lipid membranes. We highly anticipate experimental confirmation of our findings.

Conflicts of Interest
The are no conflicts of interest to declare.