Spontaneous inward opening of the dopamine transporter is triggered by PIP2-regulated dynamics of the N-terminus.

We present the dynamic mechanism of concerted motions in a full-length molecular model of the human dopamine transporter (hDAT), a member of the neurotransmitter/sodium symporter (NSS) family, involved in state-to-state transitions underlying function. The findings result from an analysis of unbiased atomistic molecular dynamics simulation trajectories (totaling >14 μs) of the hDAT molecule immersed in lipid membrane environments with or without phosphatidylinositol 4,5-biphosphate (PIP2) lipids. The N-terminal region of hDAT (N-term) is shown to have an essential mechanistic role in correlated rearrangements of specific structural motifs relevant to state-to-state transitions in the hDAT. The mechanism involves PIP2-mediated electrostatic interactions between the N-term and the intracellular loops of the transporter molecule. Quantitative analyses of collective motions in the trajectories reveal that these interactions correlate with the inward-opening dynamics of hDAT and are allosterically coupled to the known functional sites of the transporter. The observed large-scale motions are enabled by specific reconfiguration of the network of ionic interactions at the intracellular end of the protein. The isomerization to the inward-facing state in hDAT is accompanied by concomitant movements in the extracellular vestibule and results in the release of an Na(+) ion from the Na2 site and destabilization of the substrate dopamine in the primary substrate binding S1 site. The dynamic mechanism emerging from the findings highlights the involvement of the PIP2-regulated interactions between the N-term and the intracellular loop 4 in the functionally relevant conformational transitions that are also similar to those found to underlie state-to-state transitions in the leucine transporter (LeuT), a prototypical bacterial homologue of the NSS.

Members of the neurotransmitter/sodium symporter (NSS) family, which include the human transporters for dopamine (DAT), serotonin (SERT), and norepinephrine (NET), are transmembrane (TM) proteins whose main physiological function is considered to be the termination of a neuronal signal by reuptake of the neurotransmitter from the synaptic cleft into the presynaptic neuronal cell.[1−7] The uptake of neurotransmitter against its concentration gradient is coupled to the transport of Na+ into the cell.[3] The reuptake function of the NSS in neuronal signaling accords them vital roles in cognitive, emotional, and behavioral functions and in various psychiatric and neurological disorders[8] making them primary targets for antidepressants as well as for abused psychostimulants.[9,10] More recently, specific de novo mutations in human DAT (hDAT) have been associated with various neurological diseases.[11−19]

The available X-ray models of DAT from Drosophila (dDAT)[20,21] show a remarkable similarity between the TM bundles of dDAT and the bacterial leucine transporter (LeuT), for which several crystal structures corresponding to various functional states (i.e., outward open, occluded, and inward open) had been determined earlier.[22−29] Their sequence homology had suggested the same 12-TM helical bundle architecture (including 10-TM core region),[30] a centrally located primary substrate binding site (the S1 site), and two Na+ binding sites, Na1 and Na2.[31−37]

The availability of structural information had prompted the formulation of a unified mechanistic model of the transport cycle in NSS proteins that is consistent with the alternating access mechanism.[38] Structure/function studies of bacterial LeuT-fold transporters conducted in vitro with double electron–electron resonance (DEER) and single molecule FRET (smFRET),[39−43] measurements of ligand binding and transport in LeuT,[44,45] as well as computational simulations of LeuT[39−42,44,46−54] and of homology models of DAT[33−37] and SERT,[55,56] have provided important mechanistic insights into ion- and ligand-dependent conformational dynamics in NSS, suggesting that the alternating access is achieved through concerted dynamic rearrangements involving specific structural motifs of these proteins. The results of such studies have further led to a proposed mechanistic model of allosteric Na+-coupled symport in which intracellular release of the S1-bound substrate is triggered by the binding of a second substrate molecule in the extracellular vestibule located ∼11 Å above the S1 site, termed the S2 site.[44]

The focused attention on LeuT has produced a large amount of structural and functional information that underscores the importance of establishing to what extent the knowledge accumulated from this bacterial homologue can explain the corresponding functional mechanisms in the mammalian NSS transporters, at the same detailed molecular level. In spite of the great structural similarity revealed for the TM domains of LeuT and the eukaryotic NSS, it is clear that the long N- and C-terminal regions found only in the eukaryotic transporters but not in LeuT, represent a major difference that is likely to be of significant functional importance. These segments have been implicated in key mechanistic elements of NSS function including regulatory phosphorylation[57−62] and the actions of psychostimulants.[57,58,63−66] Indeed, the involvement of the N-terminal domain (N-term) in amphetamine (AMPH)-induced reverse transport (efflux) of the substrate has been well documented for the neurotransmitter transporters.[57,63,64,67−70] In hDAT, specifically, the AMPH-induced efflux has been shown to be modulated by the first (distal from the TM bundle) 22 residues in the hDAT N-term.[66] Their electrostatic interactions with highly charged phosphatidylinositol 4,5-biphosphate (PIP2) lipids[71] and phosophorylation of this region at one or multiple Ser residues (see Figure S1A–B)[57] were implicated in the mechanism of the AMPH-induced efflux.

Reasoning that an understanding of the mechanisms underlying the functional role of the terminal domains can reveal important aspects of state-to-state transitions in the NSS transporters necessary for the phenotypes that were shown to involve the N-term, we have undertaken, to our knowledge for the first time, an analysis of the possible role of the N-term in these molecular mechanisms in a full-length model of the hDAT. The analysis of the dynamics and functionally relevant conformational transitions of the transporter was carried out on extensive (>14 μs in total time) unbiased atomistic molecular dynamics (MD) trajectories of the hDAT in physiologically relevant lipid membrane environments (enriched in PIP2 lipids). The results outline initial steps of the spontaneous transition of the transporter from a starting outward-facing state to an inward-facing state, examined from rearrangements in specific structural motifs previously identified in LeuT, and culminating in the release of Na+ from the Na2 site and in destabilization of the substrate in the S1 site. We determined the manner in which these steps are regulated by the membrane environment and found that the state-to-state transitions in hDAT can be triggered by PIP2-mediated electrostatic association between the N-terminal domain and the intracellular loop 4 (ICL4). Specifically, we show here that the N-term/ICL4 complex is stabilized by PIP2 and can be formed through involvement of either the distal region of the N-term (i.e., the most N-terminal end, containing residues K3 and K5) or by the region proximal to the TM bundle (containing R51). Together, the results of this computational study provide novel insights into the role of collective dynamics in functionally relevant state-to-state transitions in the transporter, at a level of atomistic detail that enables specific experimental probing of the mechanistic hypotheses regarding essential functional mechanisms of the NSS transporters.

Results and Discussion

We have recently demonstrated both experimentally and computationally that the first 64-residue long N-terminal fragment of the hDAT associates with lipid membranes through electrostatic interactions with PIP2 lipids.[71,72] To provide a structural context for these interactions in the full-length hDAT protein and to evaluate the role of PIP2-mediated intramolecular associations in functionally relevant conformational transitions in the transporter, we constructed full-length hDAT models (see Methods) and studied their conformational dynamics with extensive unbiased atomistic MD simulations (see Figure ) in an explicit PIP2-containing lipid membrane and solvent environment. The findings from these computational studies are presented in depth below and are compared to the results obtained from the analysis of the control trajectories of hDAT collected in PIP2-depleted bilayers as described in Figure .

Figure 1

Schematic representation of the atomistic MD simulations performed on hDAT3TT1 and hDATdDAT constructs embedded in lipid membranes. Two original trajectories (Run1 and Run2) for the hDAT3TT1 model in Na+Cl– solution and a single run of hDATdDAT in K+Cl– solution were accumulated in PIP2-enriched lipid membranes (simulation times are indicated in bold at the ends of the thick blue arrows). At various time-points in these trajectories (indicated with red arrows), PIP2 lipids were substituted with POPE lipids, and the control simulations (identified as ΔPIP) were restarted. The subsequent MD trajectories are correspondingly labeled (Run1ΔPIP, Run2ΔPIPa, Run2ΔPIPb, and dDATΔPIP), and their durations are shown on top of the respective red arrows.

N-Terminus of hDAT Associates with the ICL4 Loop through PIP2-Mediated Interactions

Monitoring the dynamics of the N-terminal domain during a 1.7 μs trajectory of the hDATdDAT model showed that the N-term, which was initially placed away from the TM bundle and the lipid membrane (Figure S3B), engages with PIP2 lipids and associates with the ICL4 segment (residues H442-R443-H444) over the course of the simulation. Specifically, the basic residues in the N-term gradually establish interactions with PIP2 lipids within the first 0.2 μs of the hDATdDAT trajectory, as detailed in Figure S4A. During the same time interval, the N-term moves closer to the ICL4 loop, as the K3 residue of the N-term approaches R443 (see Figure S4B). Later in the trajectory, other basic residues in the N-term, most notably R51 (but also K35, albeit more transiently), engage with the ICL4 loop as well (Figure S4B).

Since the ICL4 is itself highly polar as it contains, in addition to R443, two neighboring His residues, we reasoned that PIP2 lipids would be involved as well in the association between the N-term and the ICL4 through electrostatic forces. Indeed, quantification of the proximity of PIP2 lipid molecules to different regions of the transporter, identified ICL4 as one of the areas that was most densely populated with surrounding PIP2 lipids (Figure S5); indeed, R443 is among the specific residues contacting PIP2 most frequently (Figure C). Of note, the ICL4 region containing R443 overlaps with residue K460 in SERT, which has been suggested to be critical for PIP2 binding[73] (R443 and H444 in hDAT align in sequence with SERT A459 and R461, respectively).

Figure 2

HDAT N-term interacts with ICL4 and PIP2 lipids in a 1.7 μs long hDATdDAT trajectory (also Figure S4). (A) Number of PIP2 lipids in contact with the interacting polar residues from the N-term (K3, K5, K27, K35, and R51) and ICL4 (H442, R443, and H444) presented as the moving average over 1 ns windows. Contact is defined by Cα-Cα distance of the interacting pair (within 9 Å), and any atom of PIP2 within 3 Å of both residues. (B) Cα-Cα distances for N-term/ICL4 residue pairs averaged over the second half of the hDATdDAT trajectory. Only the pairs within 9 Å are shown. (C) Residues within 3 Å of a PIP2 lipid headgroup in >80% of the trajectory are shown in space fill format (blue for N-term; green for ICL4). The N-term is rendered in red, and the neighboring region is shown in white cartoon format. A black arrow indicates the position of the distal segment of the N-terminus (DNt).

Analysis of the time evolution of the N-term/ICL4 association and of the contacts between the interacting residues in these segments and the nearest PIP2 lipids (Figure A) showed that PIP2-mediated interactions between the N-term and ICL4 first occur at ∼0.2 μs into the hDATdDAT trajectory and rise gradually to a value of about two contacts per time frame. A similarly strong tendency for the formation of the N-term/ICL4 complex involving PIP2 was observed in the simulations carried out with the hDAT3TT1 model (Figure S6A). Interestingly, however, in the hDAT3TT1 trajectories the N-term was observed to engage with the ICL4 segment through both its proximal (harboring R51) and its distal ends (including K3/K5) (see Figures S6B and S6D).

In the control MD simulations initiated from intermediate time-points along hDAT3TT1Rn1 and hDAT3TT1Run2 (0.46 and 2.0 μs, respectively; see Figure ), in which PIP2 lipids were replaced by POPE, the N-term/ICL4 complex that was stabilized by PIP2 dissociated rapidly (Cα-Cα distances between any pair of N-term/ICL4 residues became >10 Å). Importantly, in the additional MD trajectories in which PIP2 lipids were substituted with POPE in the initial configurations of the hDAT3TT1Run2 and hDATdDAT systems (Run2ΔPIPa and dDATΔPIP trajectories in Figure ), the N-terminus remained mostly in the solvent, away from ICL4 (data not shown). Collectively, these results suggest that the association of the N-terminus with the ICL4 segment is stabilized by the electrostatic contribution of PIP2 lipids and can involve both the distal end of the N-terminal segment (K3/K5 residues) and its proximal region (R51 residue). Note that the involvement of the K3/K5 pair in interactions with PIP2 lipids in all our simulations is consistent with the experimental findings we described[71] regarding the effect of charge-neutralizing mutations in these residues (either to Ala or to Asn) which resulted in reduced binding of the N-term to PIP2-containing liposomes (see Discussion).

PIP2-Mediated N-term/ICL4 Association Destabilizes the Functionally Important IC Network of Ionic Interactions

We next reasoned that the N-term/ICL4 association described above could influence the functionally important network of the conserved residues[36,74] in the intracellular end of the transporter, near the N-term and ICL4 segments. The available X-ray structures of LeuT and dDAT, as well as functional assays conducted on hDAT,[36] suggest that the inward-closed state in hDAT is stabilized by salt bridges formed between R60 (in the N-term) and D436 (in TM8 that connects to ICL4), as well as between Y335 (in ICL3) and E428 (in TM8). Interestingly, the dDAT X-ray structure[20] suggests that the IC network in the inward closed DAT contains an additional salt bridge, not present in LeuT, formed between E428 and the conserved R445 residue adjacent to the ICL4.

In our simulations, the IC network of salt bridges becomes destabilized at the onset of the PIP2-mediated N-term/ICL4 interactions (Figure ). We observe the disruption of the R60-D436, E428-Y335, and E428-R445 contacts at ∼0.2 μs time point of the hDATdDAT trajectory (Figure A–C), when the N-term/ICL4 complex is first formed (Figure A). As the extent of the N-term/ICL4 association strengthens (between 0.2 and 0.6 μs), a new salt bridge, formed between residues D436 and R445, is stabilized (Figure D–G).

Figure 3

(A–D) Time evolution in the hDATdDAT simulation of the minimal distance between residue pairs: (A) R60 and D436; (B) Y335-E428; (C) R445-E428; and (D) R445-D436. Panels (E–G) show snapshots of the system at 0 μs (E), 0.35 μs (F), and 1.4 μs (G) time-points along the trajectory, highlighting changes in the positioning of R60 (blue), D436 (red), Y335 (green), E428 (pink), and R445 (cyan) during the simulation. The relevant IC segments of the protein are shown in cartoon rendering. TM1a, TM8, and TM9 are labeled.

Similar tendencies are observed in the hDAT3TT1Run1 and hDAT3TT1Run2 simulations. Thus, the R60-D436 and Y335-R445 interactions break in these trajectories as well (Figure S6C) when the N-term engages with the ICL4 (Figure S6A). Concomitantly, the D436-R445 ionic interaction is stabilized (D436-R445 minimum distance is <2.5 Å in these simulations). Importantly, in the two independent control MD trajectories of hDAT carried out in the PIP2-depleted membranes (dDATΔPIP and Run2ΔPIPb), in which the N-term did not engage with the ICL4 (see above), the Y335-E428 salt bridge, suggested to be the most influential for the stability of the inward closed state,[74] was essentially retained, and the IC gates remained largely closed (Figure S7A–B). Collectively, the above data suggests that PIP2-mediated N-term/ICL4 association leads to the destabilization of the IC network of ionic interactions. Since the latter has been suggested to be of critical importance for the stability of the inward closed state of the transporter,[36] we next examined the conformational dynamics of the hDAT in response to the PIP2-mediated N-term/ICL4 interactions.

hDAT Undergoes Spontaneous Outward-Facing to Inward-Facing Isomerization in PIP2-Enriched Membranes

Recent experimental studies of conformational dynamics in LeuT[39−42] and the structurally related Mhp1[43] transporter using single molecule and ensemble spectroscopy measurements, respectively, uncovered structural rearrangements that mark the functionally relevant transitions in the LeuT-like fold NSS proteins. The distance changes between specific loci in IC and EC segments in response to ion and/or ligand binding are attributable to isomerization of these transporters between different functional states (i.e., outward-facing, inward-facing, and occluded). In describing conformational transitions in our simulations, we monitored the same sets of distances between the analogous IC and EC sites in the hDAT and quantified the accumulation of water molecules at the IC and EC vestibules of the transporter during the MD trajectories.

Markers for Opening of the Intracellular Vestibule in hDAT Are Similar to Those in LeuT

Figure shows the count of water molecules inside the EC and IC vestibules during hDATdDAT simulation, which indicates a rapid (within the first 50 ns of the simulations) loss of hydration in the EC vestibule and increase in water count in the IC vestibule with the transitioning of the transporter from the initially outward-facing state to the occluded state (Figure B). This isomerization is due to a relaxation of the protein structure to occlude the substrate DA in the S1 site from the EC vestibule. After this initial equilibration, the transporter remains in the occluded state for the next ∼150 ns of the trajectory, which is followed by large-scale concerted motions of the intracellular TM1a, TM6b, and TM9 segments. Indeed, distance measurements between various IC regions of the transporter (Figure A–C) reveal that during the subsequent ∼0.2–0.6 μs trajectory interval TM1a-TM6b, TM1a-TM9, and TM6b-TM9 distances increase by ∼4.5 Å, 8 Å, and 5 Å, respectively. As a result, the TM1a and TM9 segments swing away from the TM bundle (Figure F–G) and the IC vestibule opens, allowing the large influx of water molecules (see red trace in Figure B).

Figure 4

(A) Snapshots of the hDAT TM bundle (gray cartoon) in the hDATdDAT trajectory at 0.13 and 1.3 μs time-points. Red and blue spheres represent oxygen atoms of the water molecules in the IC and EC vestibules, respectively (see Methods for description of the water count algorithm). The substrate, DA, is shown in van der Waals rendering. The TM1a segment is labeled. (B) Time evolution of the number of water molecules in the IC (red) and EC (blue) vestibules in the hDATdDAT simulation. The green arrows denote time-points at which the snapshots in panel A were taken.

Figure 5

Time evolution in the hDATdDAT simulation of Cβ-Cβ distances between residues in various TM segments: (A) I67 (in TM1a) and L447 (in TM9); (B) I67 (in TM1a) and S333 (in TM6b); (C) S333 (in TM6b) and L447 (in TM9); (D) E307 (in TM6a) and F171 (in TM3); and (E) F171 (in TM3) and K92 (in TM1b). Panels F and G depict conformations of the hDAT TM bundle (silver) in the initial and final frames of the trajectory. The IC and EC segments from panels A–E are labeled and colored in yellow. Blue arrows in panel F indicate the direction of movement of the different regions in the transition from F to G. Collective motions of TM1b and TM6a segments are highlighted by a red dotted oval.

As seen in Figure S8, during the time period of the trajectory when the IC vestibule widens (0.2–0.6 μs) several other IC segments undergo large-scale motions. Most notably, we observe the displacement of TM7a with respect to ICL1 and TM4, as well as an increase in TM1a-ICL4 and TM2-TM5 distances. Collectively, the data in Figures A–C and S8 indicate that the hDATdDAT simulation captures an inward-opening event in hDAT that follows dynamic trends, which are consonant with DEER distance measurements[39] in LeuT. These DEER measurements showed that TM1a, NT (the fragment of the N-terminus adjacent to TM1a), TM6b, and TM7a undergo the most substantial ion- or ligand-dependent movements in the IC region of the transporter.

We note that a similar isomerization event was detected in the two 4 μs long simulations initiated from the hDAT3TT1 model (Figure ). Importantly, we found that the inward opening in these simulations followed dynamic trends largely similar to those observed in hDATdDAT (see Run1 and Run2 in Figure S9).

Inward-Opening Is Related to the PIP2-Mediated N-term/ICL4 Association

As described in Methods, we have developed protocols for quantitative analysis of correlation and clustering of time-dependent variables. Using this analysis, we found, as detailed in a section below, that the isomerization event in the hDATdDAT trajectory is related to the rise in the PIP2-mediated N-term/ICL4 interactions (see Figure A), which, in turn, have a destabilizing effect on the IC network of ionic interactions discussed above (see Figure ). Similar relationships are observed in the hDAT3TT1Run1 and hDAT3TT1Run2 simulations. For example, in the hDAT3TT1Run2 (Figure S9C), the TM9 segment moves away from TM1a and TM6b within the first 0.5 μs time interval during which the N-term/ICL4 complex is formed (Figure S6A), the IC gates brake open (Figure S6C), and water starts to build up inside the IC vestibule (Figure S9B).

Importantly, an additional 2.2 μs long MD trajectory (Run2ΔPIPb in Figure ), initiated from the 2.0 μs time-point of hDAT3TT1Run2 (inward-facing conformation) after replacing PIP2 with POPE lipids, showed movement of the ICL4/TM9 segment in the opposite direction, i.e., toward the TM bundle, as the TM1a-TM9 and TM6b-TM9 separations decreased (Figure S9C). Remarkably, we find that the inward movement of TM9 in this simulation was accompanied by dehydration of the IC vestibule (Figure S9B).

Furthermore, in the two control MD trajectories of hDAT carried out in the PIP2-depleted membranes (dDATΔPIP and Run2ΔPIPa), in which the N-term did not engage with the ICL4 (see above) and the IC gates remained closed (Figure S7A–B), there occurred neither significant repositioning of the ICL4/TM9 region with respect to the TM bundle nor an overall opening of the IC vestibule (Figure S7C–D). Collectively, the conformational rearrangements and their spatial and temporal relationships suggest that by destabilizing the IC network of salt bridges, the PIP2-mediated N-term/ICL4 interactions allow a repositioning of the ICL4 and neighboring IC regions so as to trigger an opening of the IC vestibule.

Inward Opening Is Accompanied by Concomitant Movements in the EC Vestibule

Conformational changes occurring on the EC side of the hDAT during the inward opening were evaluated from changes in distances between various EC segments.[39] The time-traces in Figures D–E and S10 show that the most prominent motions on the EC side during inward-opening involve the TM1b and TM6a segments that move closer to TM3 (Figure D–F). In addition, separation between the ECL4 and ECL6 segments increases (Figure S10E), and the change is related to the repositioning of the ECL4b motif with respect to the TM1b segment whereby F391 (in ECL4b) pulls toward W84 (in TM1b), whereas P387 (in ELC4b) pulls away from W84 (see Figure S11). Notably, these changes are consonant with experimental results reported for the structurally related Na+/proline symporter PutP,[75] which suggested that a pulling motion of the residue in PutP analogous to F391 is a trigger for inward opening. Furthermore, residue W84 in hDAT is analogous to residue W33 in MhsT, the LeuT-fold Na+/hydantoin transporter, where it has been suggested to be mechanistically important for enhancing the transitioning to the inward-facing conformation.[76]

Overall in the hDATdDAT trajectory, the collective movements result in contraction of the EC vestibule during the transition, as quantified by the decrease in the number of water molecules it accommodates (blue trace in Figure B). The specific conformational rearrangements observed in the simulations are consistent with the results from DEER spectroscopy in LeuT[39] showing that the transition to the ligand-bound state involved movements of TM1b and TM6a segments toward TM3 and a repositioning of ECL4 to favor a closed conformation of the EC vestibule. Together, the data presented in Figures , S8, and S10 establish that the inward opening in the hDATdDAT trajectory proceeds through concerted motions on the IC and EC sides of the transporter and involves structural motifs that have been identified experimentally to underlie state-to-state transitions in LeuT.

Isomerization to the Inward-Facing State in hDAT Results in the Release of Na+ Ion from the Na2 Site and Destabilization of the Substrate DA in the S1 Site

Experimental studies have suggested that an inward-closed state is stabilized by a Na+ ion bound in the Na2 site (Na2 ion) in both LeuT[44] and DAT.[55] Computational explorations of dynamics of LeuT,[44,52] DAT,[33] and SERT[56] showed that the isomerization to the inward-facing state induces a destabilization of the ion in the Na2 site. Consistent with these findings, we observe that the transition to the inward-facing state in our simulations is accompanied by the spontaneous release of the Na+ ion from the Na2 site.

The detailed dynamics of Na2 ion in the hDATdDAT trajectory represented in Figure show that during the first 0.6 μs time interval, the ion in Na2 gradually loses contact with coordinating residues in this site: L418 in TM8 and V78 in TM1 (Figure A). This destabilizes the ion in the site, and Na+ leaves the Na2 site via the intracellular vestibule. The path followed by the Na+ ion released from the Na2 site is depicted in Figure C–D, highlighting several key residues along the pathway that coordinate the Na+ ion at different time-points of the simulation. Among them, we note D421, T269, and D436, all of which were shown to have functional importance for DAT.[12,36,77]

Figure 6

Spontaneous release of the Na+ ion from the Na2 site during hDATdDAT simulation. (A) Time-evolution, during the hDATdDAT run, of minimum distances between the Na+ ion originally in the Na2 site and its coordinating residues L418 (red), D421 (green), G75 (blue), and V78 (black). (B) The distance between the two Na+ ions initially in the Na1 and Na2 sites. The insets on panels A and B detail the dynamics of the respective time-dependent variables within the first 0.7 μs time-interval. (C–D) The trajectory of the Na+ ion (initially in the Na2 site) is traced by the position of the orange spheres. The Na+ ion in the Na1 site is shown in yellow and the Cl– ion in cyan. The positions of the ions, the DA ligand (in van der Waals representation), and relevant TM segments of hDAT (in white cartoon rendering) are taken from the initial conformation of the system. The rearrangement of the various residues coordinating the Na2 ion rendered in licorice is identified by their position at different stages of the simulation. The snapshots in panels C and D are related to each other by 180° rotation around the vertical axis (direction perpendicular to the membrane plane).

The destabilization of Na2 was observed as well in the MD simulations initiated from the hDAT3TT1 model. As shown in Figure S12, the release of the Na2 ion in Run2 follows largely the same pattern as that in the hDATdDAT trajectory. Interestingly, in Run1 (Figure S13), the Na+ ion left the Na2 site but remained within the IC vestibule of the transporter for >2 μs trajectory time, ∼5 Å from the Na2 site (and 12 Å away from the Na1 ion, Figure S13B), stabilized by interactions with the side chain of D421 and the catechol moiety of the substrate DA (Figure S13D). This finding is consonant with experimental and computational studies that revealed the key role of residue D421 in coordinating the Na+ ion in the Na2 site.[12]

To evaluate the effect of Na+ release from the Na2 site on the dynamics of the substrate, DA, in the S1 site, we quantified the root-mean-square fluctuations (RMSFs) in the positional coordinates of the substrate. As shown in Figure , we found that the fluctuations in DA positioning are more pronounced in the simulations with PIP2-containing membranes than in PIP2-depleted systems in which the inward-opening and Na2 release events were not observed. In particular, in the two simulations (hDATdDAT and hDAT3TT1Run2) in which Na+ left the transporter, destabilization of the DA substrate leads to its rearrangement in the S1 site so that its ring points toward the IC vestibule (along the z axis, see Figures and S14) and away from TM8 (along the y axis, Figures and S14). In hDAT3TT1Run1, the rotation of the substrate ring toward the IC side leads to weakening of interactions between DOP and the S1 site residues (Figure S14C) and allows DOP to engage with the Na+ ion released from the Na2 site, as described above (see also Figure S13). Gradual destabilization of the substrate as Na+ diffuses out of the Na2 is consistent with the free energy calculations conducted in LeuT, which suggested that the release of a sodium ion from the Na2 site lowers the energy barrier for substrate uptake.[44]

Figure 7

Root-mean square fluctuations (RMSFs) of the substrate DA in the S1 site in the different simulations. For each trajectory, the protein structure at each frame was aligned first to the initial model using the backbone atoms of the TM bundle; the RMSF of the center-of-mass of the DA ring structure was calculated in the three Cartesian coordinates (shown in the three separate panels; see also Figure S14). The data in red and black are from simulations in PIP2-enriched and in PIP2-depleted membranes, respectively.

Correlation Analysis Shows That PIP2-Mediated N-term/ICL4 Interactions Are Allosterically Coupled to the Functional Sites in hDAT Involved in Inward-Opening Dynamics

We established the relationship between the PIP2-mediated association of the N-term with ICL4 and the sequence of rearrangements leading to the inward-opening transition associated with (i) the release of the Na+ ion from the Na2 site and (ii) the destabilization of the substrate DA in the S1 site. To this end, we clustered the time-dependent variable describing PIP2-mediated N-term/ICL4 contacts (Figure A) with the dynamic quantities discussed in Figures –6, S8, and S10–S11, using mutual information clustering (MIC). To determine the temporal relationship between the structural motifs that underlie the isomerization event in the hDAT and the PIP2-controlled N-term/ICL4 dynamics, we identified the dominant clusters (Figure ) using the FMI measure (see Methods and Figure S15).

Figure 8

Cluster tree (dendrogram) constructed from the dynamic variables discussed in Figures –6, and S8 and S10–S11. The dendrogram shows a merger of two clusters (variables and limbs shown in red and blue) and a connection of the resulting large cluster to smaller branches of the tree (in black). The variables in black are weakly correlated with those belonging to the colored clusters; as assessed by the Fowlkes–Mallows Index, the variables in black do not affect the clustering significantly, (see Methods and also Figure S15).

The resulting “cluster tree” (dendrogram) shown in Figure is constructed from the dynamic variables described in Figures –6, S8, and S10–S11. It shows a large cluster of highly correlated variables (rendered in red) consisting of distance changes on the IC and EC sides of the transporter. The cluster identifies strong temporal correlations between the dynamics of the TM1a/TM6b/ICL4/TM9 segments on the IC end, and the distance changes related to the movements of TM1b/TM6a regions on the EC vestibule combined with the movement of the ECL4b region with respect to TM1b (involving residues W84, F387, and F391). Connected to this limb of the dendrogram is another large branch (shown in blue in Figure ), which consists mostly of additional dynamic variables that describe the opening of the IC vestibule during isomerization. This branch also contains the PIP2-mediated N-term/ICL4 contacts (“NT/ICL4 cont.” in Figure ), as well as the dynamics of the Na2 ion.

Clustering analysis on the same quantities extracted from hDAT3TT1Run1 and hDAT3TT1Run2 trajectories (Figure S16) again revealed a merger of two clusters (in red and blue colors) containing various IC and EC structural motifs that describe (i) inward-opening of the transporter (most prominently, dynamics in TM1a/TM6b/ICL4/TM9 segments), (ii) movement of the Na2 ion, (iii) the extent of PIP2-mediated N-term/ICL4 contacts, and (iv) the dynamics of the ECL4b. We find that in all cases of inward-opening, there is one dominant cluster composed of interhelical distances between TM1a, TM6b, and TM9 (colored in red), which is strongly coupled to a cluster containing the N-term/ICL4 contacts (colored in blue). In this manner, the clustering analysis establishes quantitatively the coupling between the PIP2-mediated N-term/ICL4 association and the structural hallmarks related to the inward-opening transition in the hDAT.

The prominent role played by the ICL4 region in the inward opening transition prompted a deeper analysis of the manner in which dynamics in the ICL4 propagates to the functional sites of the transporter. To this end, we quantified in the simulated trajectories the extent to which collective motions within the ICL4 segment were allosterically coupled to the collective motions composed of (i) the residues that line the ion binding sites Na1 and Na2 and (ii) the residues in the primary S1 and the presumed secondary S2 substrate binding sites in hDAT. We used a measure of information transmission between multivariate systems with correlated collective fluctuations known as the total intercorrelation, which was developed as part of the N-body Information Theory (NbIT) framework we have described previously (see refs (78 and 79) for details).

The time evolution of the intercorrelation coefficient (r), which quantifies the extent of coupling between the collective motions of the ICL4 region and various functional sites, was calculated from the hDATdDAT, hDAT3TT1Run1, and hDAT3TT1Run2 trajectories and is shown in Figure . Similar plots obtained from the analysis of the simulations in PIP2-depleted membranes (dDATΔPIP and Run2ΔPIPa) are shown in Figure S17. For completeness, the Supporting Information Movie illustrates the time sequence of the corresponding intercorrelation matrix, obtained for the hDATdDAT system, quantifying the intercorrelation coefficients between all of the pairs of residues considered in the NbIT calculations.

Figure 9

Total intercorrelation coefficient (see eq ) between the residues in ICL4 and in different functional sites of the hDAT: substrate binding S1 and S2 sites (green and red, respectively); Na+ ion binding sites Na1 and Na2 (blue and black, respectively). The intercorrelation coefficients are shown separately for simulations Run1 (A), Run2 (B), and hDATdDAT (C) and were obtained as averages over 500 ns time intervals by sliding the analysis windows by 50 ns (see Methods for more details).

Figure shows the strong coupling (large r values) between the dynamics within the ICL4 segment and the dynamics in the functional sites of the transporter. The values are highest in the initial stages of the hDATdDAT, hDAT3TT1Run1, and hDAT3TT1Run2 simulations, as the initial 0.5 μs interval coincides in time with the inward-opening in these simulations (see above), but they decrease after the completion of the isomerization event. The results suggest that the isomerization event is associated with highly coupled motions in the ICL4 and the functional regions. In fact, the same analysis performed on dDATΔPIP and Run2ΔPIPa trajectories, for which the initial protein models were the same as in hDATdDAT and hDAT3TT1Run2 simulations, respectively, but in which the hDAT did not transition to the inward-facing state, revealed relatively small r values (Figure S17), comparable to those after isomerization in the hDATdDAT, hDAT3TT1Run1, and hDAT3TT1Run2 simulations. This leads to the inference that the high correlations measured in the trajectories collected in PIP2-enriched membranes are indeed related to the transition to the inward-facing state observed in these systems.

Especially notable is the higher value of r for the S2 site (red traces in Figure ). The S2 site in hDAT includes residues W84, P387, and F391, which as described above (Figure S11) are involved in the conformational rearrangements accompanying the inward opening. In these rearrangements, the specific motion in ECL4b results in a pulling of F391 toward W84, while P387 moves away from W84. Thus, we find that the observed strong correlations between the ICL4 and the S2 site is primarily due to highly coupled motions in the ICL4 segment and the S2 residues from the ECL4b loop, adjacent to P387, that participate in the described pulling motion analogous to that observed in PutP.[75] This is illustrated in the Supporting Information Movie and identified by the clustering analysis.

Figure also shows that after the initial decrease in the correlations, the coupling between the ICL4 and some of the functional sites exhibits positive fluctuations in the hDATdDAT, hDAT3TT1Run1, and hDAT3TT1Run2 systems. The timing of the increase in correlations shown in Figure coincides with the event of Na+ release from the Na2 site (at ∼0.75, 1.25, and 1.0 μs time-points in the 3 simulations hDATdDAT, hDAT3TT1Run1, and hDAT3TT1Run2 systems, respectively). After the release is complete, we observe a slow relaxation of the correlations (see hDAT3TT1Run2 in Figure ), with the coupling between ICL4 and the S2 site remaining the highest among all the correlations considered. Collectively, the NbIT analysis demonstrates and quantifies the allosteric coupling between the ICL4 and the functional sites during the PIP2-mediated inward-opening of the hDAT, and identifies the S2 site as the region with the strongest coupling to the ICL4, suggesting that this distant communication is mechanistically important for the inward-opening transition.

N-Terminus as an Important Allosteric Modulator of the Functional Inward-Opening in hDAT

The results presented above addresses, to our knowledge for the first time, the mechanistic involvement of the N-terminal region of the hDAT in the functionally relevant conformational transitions of the transporter involved in the inward-opening of the hDAT. Thus, novel molecular-level mechanistic insights about the PIP2 lipid-mediated dynamics of the hDAT N-terminus in functionally relevant state-to-state transitions in the transporter suggest specific steps of the spontaneous transition of the hDAT from the initially outward-facing conformation to the inward-facing state. Interestingly, the rearrangements observed in specific structural motifs parallel those identified previously in LeuT and are found to culminate in the release of Na+ ion from the functional Na2 site and the destabilization of the substrate in the S1 site.

The results show that the conformational isomerization triggered by the strong tendency of the N-term to associate with the ICL4 segment through PIP2-mediated electrostatic interactions involves the H442/R443/H443 triad of polar residues in the ICL4 and the set of basic residues in the N-term. The mechanistic consequences of the PIP2-mediated N-term/ICL4 association that emerge from this analysis are the disruption of a conserved IC network of ionic interactions, which triggers the inward-opening by destabilizing the IC network of ionic interactions, and the associated release of the Na+ ion from the Na2 site causes destabilization of the substrate DA in the primary S1 site.

Taken together, the simulation results suggest how the N-term/ICL4 association driven by the interaction with PIP2 facilitates the isomerization of hDAT to the inward-facing state. The structural hallmarks of this isomerization in our simulations are similar to those identified experimentally to underlie functionally relevant transitions in the LeuT-like NSS proteins. Moreover, results implicating interactions between the N-term and ICL4 segments in the inward-opening isomerization of the hDAT construct are consonant with experimental studies in the GABA transporter (GAT1) in which the same two regions have been shown to modulate forward as well as reverse substrate transport.[80,81]

The atomistic details of the computational modeling show that the mechanism for the N-term interaction with the PIP2-containing membrane involves charged residues both at the distal (i.e., the very N-terminal part) end of the segment (residues K3 and K5 in particular) and the proximal end (near TM1 of the transmembrane bundle) of the segment (and especially R51). It is known that specific charge-neutralization mutations in the N-term, either K3A/K5A or R51W, do not affect regular substrate transport but do interfere with amphetamine-dependent efflux.[11,71] The simulations identify the contacts between PIP2 and R51 as leading to the association of R51 with the ICL4 residues and show that the interactions between PIP2 and K3/K5 lead to an anchoring of the distal part of the N-term to the membrane surface (as in the hDATdDAT and hDAT3TT1Run2 simulations), which has been suggested to enable the N-term phosphorylation required for AMPH-induced efflux.[71] In the hDAT3TT1Run1 trajectory, this drives N-term/ICL4 complex formation.

These detailed observations from the multiple simulations can explain why the experimentally measured substrate uptake is not necessarily impaired when efflux is reduced by mutation. In the hDAT3TT1Run1 trajectory, it appears that the charged K3/K5 distal segment can also drive N-term/ICL4 complex formation in the absence of phosphorylation so that the PIP2-mediated N-term/ICL4 association may be sustained in the hDAT construct bearing either the K3A/K5A or the R51W mutation separately and that both the K3A/K5A hDAT and the R51W hDAT constructs could transition toward the inward-facing state through this mechanism.

The consequences of these interactions for the functional mechanism of the transporter are underscored by our findings showing that inward opening is accompanied by concomitant movements in the EC vestibule and that isomerization to the inward-facing state in hDAT results in the release of the Na+ ion from the Na2 site and the destabilization of the substrate (DA) in the S1 site. Indeed, our analysis using the tools offered by the NbIT framework[78,79] identified the allosteric connection showing that collective motions triggered by the N-term/ICL4 association on the intracellular side are strongly coupled to collective motions in the extracellular vestibule and in the substrate and ion binding sites. Further substantiating the mechanistic importance of the PIP2-mediated N-term/ICL4 interactions, their effect on the intracellular side emerged clearly from the MIC results showing that they are allosterically coupled to the functional sites in hDAT involved in inward-opening dynamics. Thus, the MIC revealed a strongly coupled helical bundle (composed of TM1a, TM6b, and TM9) in the intracellular side that was highly correlated to the N-term/ICL4 association, suggesting how the N-term can modulate the stability of this bundle and thus modulate intracellular gating. In combination, these results identify the N-terminus as an important allosteric modulator of the functional inward-opening and ion/substrate release in hDAT.

Direct effects of the N-term interaction with PIP2 in the function of hDAT had been inferred from the various experimental studies cited throughout. Some of these studies have demonstrated that the PIP2-sequestering basic peptide (palpeptide) competes with the hDAT N-term to bind PIP2-enriched liposomes[71] and that cell perfusion with palpeptide impairs AMPH-induced DA efflux but not the regular DA uptake.[71] However, the interpretation of such results is rendered somewhat difficult by the possible involvement of other proteins and various plasma membrane components that are known to regulate functional mechanisms in hDAT; specifically, palpeptide treatment is very likely to affect as well the spatial organization of such membrane components. Indeed, colocalization of hDAT with effector proteins such as Flot1 and STX within plasma membrane microdomains is known to be important for AMPH-induced substrate efflux but not for the uptake.[70] In turn, STX in neuronal cells is found to reside inside PIP2-enriched plasma membrane compartments,[82] suggesting that the effects of PIP2-sequestering palpeptide may stem, at least in part, from disruption of these functional domains. These considerations bring to light the importance of future attempts to broaden the scope of structure–function studies to include rigorous molecular-level mechanistic investigations of DAT and the NSS molecular machines within the rich environments that are relevant to their function.

Methods

Molecular Modeling of hDAT

Homology Model of the hDAT Transmembrane Bundle

Several molecular models of the full-length hDAT (residues 1–620) were prepared for all-atom MD simulations in explicit lipid membranes and water environments. Briefly, we used Modeler 9v10[83] and a previously published sequence alignment of the NSS-family proteins[31] to first construct homology models for the TM part of the hDAT (contained in residues 57–590) based on either the outward-facing dDAT structure (PDB code: 4M48),[20] or on the high resolution outward-open X-ray structure of the bacterial member of the NSS family, LeuT (PDB code: 3TT1).[22] As described in detail,[12,55] the models included the substrate, dopamine (DA), docked in the S1 site, two Na+ ions, positioned equivalently to those in the LeuT crystal structure, and a Cl– ion coordinated by residues Asn82, Tyr102, Ser321, and Asn353 of hDAT, based on the chloride binding site described previously.[32,84]

PROPKA 3.0[85] was used to determine the protonation states of ionizable residues of hDAT and resulted in Glu491 being protonated. Furthermore, a disulfide bond was introduced in EL2, between Cys180 and Cys189, according to available structural and functional information. More details about the models can be found in refs (12 and 55).

Ab Initio Modeling of hDAT Termini

Structural information regarding the terminal domains of the DAT or of other NSS proteins is lacking. Indeed, the X-ray models of the dDAT does not include the N-term, whereas the C-term is only partially resolved. Therefore, we generated and tested 3D models of the structurally unknown N- and C-terminal domains of the hDAT (fragments 1–57 and 591–620, respectively, that lack sequence homology to proteins of known fold) as described earlier,[72] using Rosetta-based ab initio structure prediction algorithms. Briefly, different fragments of the termini were subjected to the Rosetta ab initio fold prediction routine, and for each construct, the predicted structures were clustered under various residue exclusion conditions. Clusters containing the majority of structures were identified. The conformations in the top clusters were evaluated with the RMSDTT iterative fitting algorithm[86] to find regions with the highest structural conservation within each cluster, and the folds with the lowest scores (from the Rosetta energy function) in each cluster were selected.

As described in our recent work,[72] the protocol described above identified two structured elements in this fragment (Figure S1), an α-helical region and a β-sheet motif, that were stable in long time-scale (∼2.2 μs) MD simulations. When the same set of computational tools was applied to the folding of the C-terminus of hDAT, α-helices were predicted in this region as well (Figure S2), consistent with the X-ray structure of the dDAT in which a large part of the C-terminus that is present is helical. More details of the ab initio modeling of the terminal domains are presented in ref (72).

The predicted structures for the N- and C-termini were docked (using Modeler 9v10[83]) onto the two models of the hDAT TM bundle described above to complete the full-length hDAT models based on dDAT and LeuT (referred to throughout as hDATdDAT and hDAT3TT1, respectively; see Figure S3). Since the X-ray model of the dDAT does not inform on the mutual orientations of the two termini, we considered for the hDAT3TT1 model two alternative docking poses, resulting in two starting conformations for subsequent MD simulations, in which the relative positioning of the N- and C-termini were different (see Figure S3A). For the hDATdDAT model, only one docking pose was considered in which the positioning of the C-terminus closely followed that in the dDAT X-ray structure, and the N-terminus was docked so as not to contact any residue in the TM bundle (Figure S3B).

hDAT-Membrane Complexes

hDATdDAT and hDAT3TT1 models were immersed into a pre-equilibrated membrane[72] containing an asymmetric lipid distribution of 451 lipids between the two leaflets so as to resemble a lipid composition of neuronal cell plasma membranes:[87] 100:40:32:27:29 mixture of POPE/POPC/PIP2/POPS/cholesterol on the intracellular leaflet and 176:29:18 mixture of POPC/DPPC/cholesterol on the extracellular leaflet.

For each transporter-embedded membrane patch, lipids overlapping with the protein were removed so that the overall lipid fractions were retained. After solvating with TIP3P water, the transporter–membrane complexes were neutralized with either K+Cl– or Na+Cl– as specified in Results and Discussion (see Figure captions), resulting in a final atom count of ∼150,000.

Molecular Dynamics Simulations

Simulations of the hDATdDAT and hDAT3TT1 constructs in the corresponding membrane environments were initiated with the previously established multistep equilibration protocol[44] carried out with NAMD software, version 2.9.[88] During this stage, the backbone of the protein was first fixed and then harmonically constrained. The solvent was initially prevented from entering the lipid–water interface. The constraints on the protein backbone were released gradually in three steps of 300 ps each, changing the force constants from 1, to 0.5, and 0.1 kcal/ (mol Å2), respectively. This step was followed by relatively short (50–100 ns) unbiased MD simulations performed with a 2 fs integration time-step and under the NPT ensemble (at T = 310 K), using the Particle-Mesh-Ewald (PME) method for electrostatics[89] and the Nose–Hoover Langevin piston[88] to control the target 1 atm pressure, with Langevin piston period and decay parameters set to 100 and 50 fs, respectively.

After this equilibration phase, long, microsecond-scale unbiased MD simulations were initiated on the Acellera GPU cluster that runs the specialized MD simulation software ACEMD.[90] ACEMD allows computations with standard CHARMM force fields, and for all the runs (including the equilibration phases with the NAMD described above), we used the all-atom CHARMM27 force field for proteins with CMAP corrections,[91] the CHARMM36 force field for lipids,[92] the TIP3P water model, and the CHARMM-compatible force-field parameter set for PIP2 lipids.[93] For the simulations with ACEMD, we implemented the PME method for electrostatic calculations and carried them out according to the protocol developed at Acellera[90] with a 4 fs integration time-step with a standard mass repartitioning procedure for hydrogen atoms implemented in ACEMD. The computations were conducted under the NVT ensemble (at T = 310 K), using the Langevin Thermostat with Langevin Damping Factor set to 0.1. More details about MD simulations with ACEMD can be found in ref (90).

The production times for the two runs with the hDAT3TT1 model and the run with the hDATdDAT construct are listed in Figure . As detailed in Figure , at different time-points along these trajectories we depleted the PIP2 lipids by replacing them with POPE lipids. After short (5000 step) minimization, these new systems were subjected to long unbiased MD simulations (see Figure for time-scales). All the trajectories were analyzed using a combination of VMD,[94] R software,[95] and various analyses tools developed in-house.

Calculating the Number of Water Molecules in the EC and IC Vestibules of hDAT

Following an algorithm used in our previous work,[46] a water molecule is considered to be in the EC vestibule if (a) its oxygen atom is within 26 Å of the Cβ atom of F326 but not within 5 Å of lipid atoms and (b) if its z-coordinate is larger than the Cβ atom of F326 by no more than 23 Å (the z-axis is perpendicular to the membrane with its positive direction toward the EC side). Similarly, the algorithm assumes that a water molecule belongs to the IC vestibule if its oxygen atom is within 15 Å of the ligand dopamine but not within 5 Å of any lipid and if its z-coordinate is larger than that of the Cβ atom of D436 by no more than 15.5 Å.

Correlation and Clustering Analysis of Time-Dependent Variables

The temporal correlation between time-dependent variables extracted from the MD simulations was quantified by calculating the Pearson correlation coefficients between the pairs of variables. To cluster the dynamic quantities based on the strength of temporal correlations between them, we then performed agglomerative mutual-hierarchical clustering[96] on the matrix of correlation coefficients, using mutual information as the distance criterion, as described in ref (96). The linear approximation of the mutual information was calculated and corrected for dimensionality using the generalized correlation coefficient.[97] Briefly, the clustering algorithm first assigns each variable to its own branch and then calculates the pairwise generalized correlation coefficient between all branches. Two branches with the highest correlation are then merged, and the generalized correlation coefficient between the new branch and all other branches are recalculated. The algorithm continues until all variables are members of the same branch. In our application here, the dendrogram (tree) is built from the small number of medium-sized clusters consisting of highly correlated variables. The moderately correlated clusters are then merged into one large cluster. Finally, the tree is completed with the remaining small clusters that are weakly correlated to both each other and the large cluster.

In order to distinguish quantitatively between mergers of different sized clusters, we used the FMI[98] as a measure of similarity between the different clustering before and after each merger. The FMI is defined as

In the above equation, n is the number of objects being clustered, and m is a matrix with rows equal to the number of clusters in the premerger clustering, columns equal to the number of clusters in the postmerger clustering, and elements m equaling the number of common members in cluster i of the premerger clustering and cluster j of the postmerger clustering. The FMI ranges from 0 to 1, with low values indicating the merger of similarly sized clusters that significantly change the clustering and high values describing mergers that do not change the overall clustering, such as the joining of small clusters (compared to other existing clusters) to each other or to a larger cluster. Therefore, the mergers of the medium-sized clusters will have relatively low FMI, whereas the mergers of small clusters to the large cluster (as is observed at the end of the clustering) will have a high FMI.

Using NbIT to Calculate the Time-Dependent Total Intercorrelation Coefficient

In order to calculate the time-dependent total intercorrelation coefficient, we used the computational framework in the recently developed NbIT (N-body Information Theory) analysis methodology.[78,79] Briefly, subsets of the full trajectories are generated using sliding windows, beginning a new 500 ns subset every 50 ns. For each subset, the total intercorrelation between specific regions of the protein is calculated as described in ref (79). First, for each set of residues of interest, the set is considered as a 3N-dimensional multivariate distribution, where N is the number of non-hydrogen atoms in the set, and the multivariate normal approximation of the differential entropy H is calculated from the non-hydrogen atomic fluctuation covariance matrix using Carma[99] and in-house software. Next, the total intercorrelation between set A and B is calculated aswhere TC(A,B) is the total correlation of the joint distribution of A and B, and TC(A|B) and TC(B|A) is the total correlation of A conditional on B and the total correlation of B conditional on A, respectively. The total correlation is defined aswhere d is the dimension of distribution A, and A is the i-th dimension of A. The conditional total correlation is defined asThe total intercorrelation is then expressed as the total intercorrelation coefficient, r, which ranges from from 0 to 1