The binding modes of CHL in the pocket identified previously (Mathieu et al., 2002; Clark, 2020) were explored with extensive atomistic MD simulations comprised of 12 statistically independent replicates of 33.6 μs each for a total of 403 µs, starting from the crystal structure of StarD4 (PDBid: 1JSS) with a CHL molecule docked in the binding site using the Schrodinger Induced Fit Docking protocol (Sherman et al., 2006a; Sherman et al., 2006b) (see Methods for details). The simulations produced two different modes of CHL binding in the hydrophobic pocket (Figure 1B), one with the hydroxyl group located near Ser136/Ser147 (termed “Ser-binding mode”), the other with the CHL OH near Trp171 (“Trp-binding mode”). In the “Ser-binding mode”, the cholesterol hydroxyl group forms hydrogen bonds with the Ser136 and Ser147 sidechains either directly, or through a water bridge (Figures 1C, D). In the “Trp-binding mode”, the cholesterol hydroxyl group interacts either directly with Trp171, or resides between Trp171, Arg92 and Tyr117 with which it interacts through H-bonds mediated by a water molecule (Figures 1E, F). Spontaneous transitions of cholesterol between binding modes are observed in more than half of the trajectories (Supplementary Figure S1).

To reveal the sequence of conformational events that accompany the transitions between the Ser and Trp binding modes, we applied the recently developed Rare Event Detection (RED) protocol (Plante and Weinstein, 2021) as described in Methods. The mechanistically important state-to-state transitions of StarD4 in response to the translocation of ligand cholesterol (see ref. 31 for a discussion of the relation between rare events and function) were extracted from the dynamic information in an ensemble of 6 trajectory stretches (4–8 μs each, 34.4 μs total) from different trajectory replicas in which a stable cholesterol transition from one binding mode to another was observed. This ensemble of 6 trajectories served as the input training data for the RED protocol in identifying the rearrangement events shared among the different independent replicas.

The RED algorithm decomposed the trajectories into 5 components based on the evolution of the residue-residue contact map over time (see Methods). As described in ref (Plante and Weinstein, 2021), each component represents a state of StarD4 in which a particular conformation becomes dominant, at a particular time in the trajectory. This time-ordered series of events in which different components dominate the structural changes encoded in the trajectory, is shown in Figures 2A,B. Several components identified by the RED algorithm are observed to be dominant concurrently with the repositioning of cholesterol in the binding site (Figures 2A–D). This is remarkable, because the information provided to the RED algorithm consists only of the contact matrices for each frame in the trajectory, without direct information about CHL and its position. The correspondence between the binding modes and particular RED components is illustrated in the time frames by the correspondence of the times when component A (“comp A” in Figures 2A,B) is dominant, and CHL is in the Ser-binding mode of CHL (Figures 2C,D). Moreover, as the dominance of the A component wanes and is replaced by that of comp-s B and C, CHL is in the Trp-binding mode. Thus, in every independent trajectory the transition events from comp A to either comp B or comp C coincide with the translocation of cholesterol (Traj0&3 shown in Figures 2A–D, and Traj1,4,6,7 in Supplementary Figure S2A). This suggests a mechanism of coupled dynamics in the cholesterol-StarD4 complex, connecting the structural rearrangements of the protein frame with the transition of CHL in the binding pocket. The states associated with the transition of CHL to the Trp-binding mode are of particular interest because in this position the CHL is near a putative “exit gate” from the binding pocket.

We propose that the specific structural rearrangements of the protein frame constitute a preparatory step for a subsequent functional docking of the loaded StarD4 to the membrane for the release of its cargo. A key indication for the likelihood of this mechanism is the consistency in the different simulation replicas of the conformational changes coupled to the transitions between CHL binding modes.

RED analysis reveals the specific conformational changes that occur simultaneously: The conformational changes that are coupled to the CHL translocation toward the vicinity of the “gate” are identified from comparisons of the spatial arrays in the RED analysis. For this comparison we identify structure-differentiating contact pairs (SDCPs) that underlie the conformational changes that occur simultaneously during the rare event of transition from one component to the next. While numerous residue contacts included in a component may change from “on” (i.e., in contact) to “off” (i.e., no longer in contact; see protocol in Methods) an SDCP must satisfy two simultaneous criteria: the pair makes a substantial contribution to the contact feature in one component and a minimal contribution in the subsequent component in the transition. Both the substantial contribution, and the minimal contributions, are determinant elements of a particular component’s conformational features, because the change in structure depends on both the formation of new contacts and the release of contacts from a previous conformation. Thus, they collectively determine the conformational change characteristic to the rare dynamic event which makes them essential for interpreting the RED results.

Comparison of the time axis in panels A and B of Figure 2 with that in panels C and D shows the transition between components A and B, and A and C occurs at the same time as the CHL transitions from the Ser-binding mode to the Trp-binding mode. The structural motifs involved in this transition are identified by the quantification of contributions to the components shown in Figure 2E-G, and defined structurally in Figure 3.

The SDCPs in Figures 2E–G are grouped (from 1 to 8) by their location in a particular secondary structure motif. Specifically, group 1 includes contacts between β1 β2 & the C-terminal Helix (H4); group 2 includes contacts between Ω1 & H4 (see Figures 3B, E); group 3 includes contacts between Ω1 & β9 (Figures 3B, E); group 4 has the contacts between β9 & β8 (Figures 3C, F), and group 5 includes the interactions in the first turn of H4 (Supplementary Figure S4A).

The complete list of SDCPs is given in Supplementary Table S1, which summarizes the structural characteristics of the system at different stages in the trajectory. The structural analysis shows that SDCPs in groups 1 to 5 are in contact in comp A in which the CHL is in the Ser-bound mode (orange, Figure 2E), but most of them break contact during the transition from comp A to comp B or C, and thus have low feature values in comp B and comp C where CHL is in the Trp-binding mode (orange, Figure 2F, G).

Another set of SDCPs, 6 to 8, have low values in comp A (blue and purple in Figure 2E) but high values in comp B and C (blue in Figure 2F and purple in Figure 2G). They are found to represent contacts formed only after the repositioning of the CHL. More generally, comp B and comp C share a highly similar (but not identical) set of structural specifics represented in the salient SDCPs in groups 6 to 8. The highly similar pattern of motif interactions represented by these shared SDCPs include the following: in group 6: interactions between β1 β2 & H4 (Figures 3A, D); in group 7: interactions between β2 & β3 & β9; in group 8: interactions of β8 or β9 with Ω1 or H4. Thus, the summary in Supplementary Table S1 shows that groups 6-to-8 in comp B (B6-B8) and groups C6-C8 in comp C share most of the key secondary structural elements and their conformations.

The major difference between comp B and comp C is the partial unfolding of the N-terminal head of H4 (group 5; residues 199–201) which loses some contacts in comp B but remains structured in comp C. Together, the structural features of components B or C compared to those of comp A reveal the specific conformational changes that are captured as “rare events”.

The main rare events identified by the transition from a conformation in which one component is dominant to a new conformation with a different dominant component, are the A→B and A→C transitions which are concurrent with the CHL translocation (Figures 2A–D). The comparison of the structural features in comp A and in comp-s B&C shows that the conformational changes concurrent with the CHL translocation are the opening the access to the hydrophobic pocket. These conformational changes include (1) the opening of a front gate by the distancing of the H4 N-terminus from the Ω1-loop and strengthening interactions with β1 (see Figures 3A, B, D, E, and the detailed residue interaction data in Figures 8A,B); the opening of a back corridor between the β9 tail and the β7β8-loop (Figures 3C, F).

To facilitate analysis of the metastable states and transition pathways encoded in the MD trajectory data, we built the conformational landscape of the cholesterol-bound StarD4 by applying the dimensionality reduction approach tICA: time-structure based independent component analysis. To describe the dynamics of cholesterol binding modes and of the protein conformational changes (see Methods and Supplementary Table S2) we chose a set of 7 collective variables (CVs) to define the tICA space. The analysis showed that the first 3 tICA vectors describe 81% of the total dynamics of the system in the 403 μs production run trajectories (Figures 4A,B). The projections on the space spanned by these three vectors are shown in Figure 4C.

To create a Markov State Model (MSM) and analyze the transitions among the metastable states visited by the StarD4 protein, the 3D tICA space was first discretized into 200 microstates using k-means clustering (see Methods and specifics in Supplementary Figure S5). The 200 microstates were then segregated according to similarity in kinetic properties, resulting in the identification of 9 macrostates. The 5 macrostates with high population are shown in Figure 5A, and representative conformations obtained from structural analysis are shown in Figures 5B–D and Supplementary Figures S6, S8.

The structural interpretation of macrostates 0–4 in Figures 5B–D and Supplementary Figures S6, S8 indicated that the differences reflect the dynamics of gate, corridor, and cholesterol we detected in the stable transition events, as indicated by the data in Figure 5E showing the relevant CVs, and underscoring that these differences prevail throughout the conformational space.

Specifically, macrostate 1 in Figure 5B represents conformational states where cholesterol is bound in the Ser-site mode (with low “CHL-Ser-dist” and high “CHL-Tyr-dist”, seen in the first and second panel of Figure 5E), StarD4 is in the crystal-like conformations with a closed H4-Ω1 gate (high “β1-H4-dist”, low “H4-Ω1-dist” in the third and fourth panel of Figure 5E) in which the β8-β9 corridor is closed (shown by low “β6-β9-dist” and high “β9-β2-β3-dist”&“β23loop-H1β8-dist” in the fifth to seventh panel in Figure 5E). Macrostate 0 also represents conformational states with Ser-binding CHL in crystal-like StarD4 with gates closed, but with a conformation of the β23 loop that is extended away (Supplementary Figure S6A) (high “β23loop-H1β8-dist”, seventh panel, Figure 5E) which may hinder the opening of β8-β9 corridor. These Ser-binding/gates-closed conformations in Macrostates 1 and 0 contribute 30% of the populations in the conformation space.

The most populated macrostate is Macrostate 3, with 47% of the population. It contains conformational states in which CHL is in the Trp-binding mode with both the H4-Ω1 gate and β8-β9 corridor open, and with β9 leaning towards β2&β3 which rearrange into narrower conformations (Figures 5D, E). Notably, all these observed rearrangements are consistent with the conformational changes detected by the RED analysis of the corresponding RED event. In Macrostate 4 the cholesterol is in a Trp-binding binding mode or even further down towards the H4-Ω1 gate that opens wider (Figure 5E, Supplementary Figure S6B). As indicated by the analysis, Macrostates 0,1,3,4 constitute 88% of the conformational space containing the concurrent dynamic rearrangements of the cholesterol binding mode, the opening H4-Ω1 gate, and the β8-β9 corridor.

The other structural motif movements that were most important in the Comp A to Comp B/C transition event detected with RED are also evidenced by the comparisons of CVs related to the transition between the Ser- and Trp-binding modes of CHL. These are quantified in Supplementary Figure S7 and include the unfolding in H4 as well as the straightening of H4, the narrowing of β2 and β3, and residue interactions between β1-H4 and disassociations between β8-β9. Thus, the conformational changes of the StarD4 protein that we found to occur simultaneously with the transitions between CHL binding modes, yield conformational states that are specific to the position of CHL in the binding pocket.

The time sequence of cholesterol translocation in the binding site is embedded in the MSM transition pathways starting from Macrostate 1 (crystal-like StarD4 with Ser-binding cholesterol) and ending in Macrostate 3 (gates-opened StarD4 with Trp-binding cholesterol). The fluxes along the most probable pathways (contributing ∼95% of the transition flux) are summarized in Table 1. Each of the top two most probable paths contributes ∼30% of the total flux. One is the direct transition from Macrostate 1 to Macrostate 3 which encompasses the simultaneous translocation of CHL and conformational change of StarD4, and the other involves an intermediate, the less populated Macrostate 2 in which CHL is in Ser-binding modes but with both H4-Ω1 gate and β8-β9 corridor open (Figures 5C, E). These results are consistent with the structural analysis in a report (Tan et al., 2019) showing that conformational changes at the H4-Ω1 gate are required for cholesterol to bind lower in the hydrophobic pocket, due to steric hindrance.

The less populated intermediate macrostates are summarized in Supplementary Figure S8. They indicate the flexibility of the StarD4-cholesterol complex and support the relation between CHL position and the gate and corridor dynamics. Thus, with CHL in the Ser-binding mode, the gates-closed crystal-like conformations represent 30% of the population, much more than the ones with either the H4-Ω1 gate or the β8-β9 corridor open (10%). Most likely this is due to the energy cost of exposing the StarD4 hydrophobic pocket and cholesterol’s hydrophobic tail to water. In contrast, when CHL is in the Trp-binding mode the H4-Ω1 gate is always open, and the conformations with both gate and corridor open amount to 58% compared to those with β8-β9 corridor closed, which are very unlikely (1.7%).

Having identified a set of concurrent conformational changes in several structural motifs that include the cholesterol binding site, we applied NbIT analysis (see (LeVine et al., 2014; LeVine and Weinstein, 2014) and Methods) to reveal the allosteric network that connects them. NbIT analysis calculates the normalized coordination information (NCI), which describes how much information of the receiver system can be gained from information about the transmitter system (see (LeVine and Weinstein, 2014) and Method for details). Here we applied the NCI analysis to the trajectories from the StarD4-cholesterol complex simulations to quantify the information transmission between motifs that emerged from the RED analysis described above. The residue index in Method identifies the components of the structural motifs used in this analysis summarized in Table 2. In Supplementary Table S3 we report the structural motifs that were used as negative controls for which the NCI quantification is expected to reveal only background-level sharing of information with the cholesterol binding site.

The results from the NCI calculation show that the β1, β2β3, H4 head, and Ω1 loop motifs are strongly coordinated with the cholesterol binding site (S136, S147, W171, R92, Y117), while weaker coordination is found with the β9 and the β7β8 loop regions. Notably, the H4 head and Ω1 loop constitute the H4-Ω1 gate, and the β9 and the β7β8 loop constitute the β8-β9 corridor of the hydrophobic pocket. The stronger coordination of the CHL binding site by the H4-Ω1 gate than β8-β9 corridor is in accordance with the structural analysis in the previous section showing that the conformational changes and interaction changes at the H4-Ω1 gate always accompany CHL translocation to the Trp-binding state (Figures 3, 5E)

Because β1 and β2β3 are also strong coordinators of the binding site despite the long distance separating these motifs, we calculated the normalized mutual coordination information (NMCI) to detect the coordination channel for information transmission from β1 to the CHL binding site. The calculated NMCI shows that β2β3, H4, and β9, all share the information with the binding site, which identifies them as components of the coordination channel shown in Figure 6A. These results show how the coordination of cholesterol binding mode by β1 is mediated through the interaction of β1 with β2β3 and H4, which stabilizes the β9 and H4 in the gate opening conformation, as summarized in Figure 6B

The release pathway of CHL in the yeast sterol transport protein Osh4 (Singh et al., 2009) was shown to comprise step-wise interactions of the CHL hydroxyl with a sequence of hydrophilic residues inside the hydrophobic pocket. In StarD4 we observed a similar sequence of CHL occupying binding sites along the pocket, composed of the S136&S147-site (Figure 1D), a water-mediated S136&S147-site (Figure 1E), the W171-site (Figure 1F), and a water-mediated W171, R92 & Y117 site (Figure 1G). Given the coordination between CHL translocation in the binding site and specific changes in the conformation of the StarD4 protein identified from the detailed analysis of MD simulation trajectories we report, we probed the impact of impairing such a binding site by introducing a W171A mutation. We used a Free Energy Perturbation (FEP) protocol to introduce the mutation in representative conformations of the two most populated Macrostates (Figure 5): State 1 with CHL in a Ser-binding mode, and State 3 where CHL interacts with W171 and is also involved in the water bridged interactions with R92 and Y117. From each binding mode we launched 24 replicate runs of 2,720ns/each, resulting in an ensemble of simulations totaling 131μs for the W171A-StarD4 mutant.

In simulations starting from the Ser-binding state, we find that the W171A mutation eliminated the translocation of CHL away from the Ser-sites, as the ligand is either in the Ser-binding mode or in the water-bridged Ser-binding mode. In the simulations started from the Trp-binding-like state where the CHL is interacting with R92 and Y117 (through a water-bridge), the W171A mutation destabilizes this binding mode and the CHL tends to translocate from the Trp-binding-like state back to the Ser-binding site. The backward transitions are observed in 8 out of the 24 replicas, with the Ser-binding conformation contributing 21% of the population in the ensemble started from the Trp-binding-like mode (Supplementary Figure S9). In comparison, stable backward transitions were observed only 2 times in the 403 μs simulations of the WT StarD4 (Supplementary Figure S1). Thus, W171 turns out to be an important intermediate site that enables accessibility and stability for CHL binding to the downstream binding sites around W171, R92 and Y117.

Notably, structural analysis suggests that the allosteric connection persists in the W171A StarD4-CHL complex and the conformation of StarD4 is compatible with the CHL in the Ser-binding mode. Thus, when CHL transitions back from the R92 Y117 water network to the Ser-binding site, the initially open Ω1-H4 gate closed, and the unraveled N-terminal of H4 folded back (Figure 7; Supplementary Figures S9B–D). These conformational changes are consistent with the allosterically connected changes in WT StarD4 where the Ω1-H4 gate opening and the unfolding in H4 are found to be dominant features in the RED analysis only in the Trp-binding mode (Figure 5E; Supplementary Figure S7). Results from the NbIT analysis of the mutant confirm that the allosteric effect remains strong between the CHL and binding site to the peripheral transmitter motifs around H4 (Supplementary Table S5).

The mutation was chosen based on the results from analysis with the RED algorithm which revealed that the transitions of CHL in the binding site are coupled to a shift of the K49-S208 interaction to a K49-T204 configuration, and the establishment of a V51-A207 contact (Supplementary Table S1). According to the allosteric model created with NbIT, β1 is an allosteric transmitter that coordinates the dynamics of bound CHL through its interaction with H4 and β2. The H4-β1 interaction involved in H4 repositioning is stabilized by both electrostatic and hydrophobic interactions. Thus, when the H4-β1gate is closed the K49 sidechain relocates near S208 and its hydrocarbon chain interacts with A207 (Figure 8A). When the H4-β1gate opens, the K49 sidechain shifts closer to T204, enabling V51 to form a contact with A207 and remove this hydrophobic residue from an unfavorable aqueous environment (Figures 8A, B). The NbIT coordination contribution analysis confirms that K49 is contributing the most to the allosteric interaction between β1 and the CHL binding site residues. Therefore, we mutated K49 expecting changes of the β1-H4 interaction mode that weaken the allosteric coordination.

Mutations K49A and K49W were introduced with the same FEP protocol as above into the StarD4-CHL complex in Microstates 1 (Ser-binding mode) and 3 (Trp-binding mode). Then, simulations were carried out for 24 replicas from each state for 2,320ns per replica of the K49A-StarD4 (111μs total), and 2,720ns per replica for K49W-StarD4 (131μs total). When started from the Ser-binding mode, the CHL transitioned to the Trp-binding site in 4 out of 24 trajectories for K49W-StarD4 (containing 7% of the population, Supplementary Figure S10A), and in only 1 out of 24 trajectories for K49A-StarD4 (containing 2% of the population, Supplementary Figure S10B). These results suggest that the dynamics of CHL are significantly slower in the mutants compared to the WT StarD4 in which 6 out of 12 trajectories sampled the transition within the first 2720ns. No backward transition from a Trp-binding state to a Ser-binding state was observed in the K49A and K49W mutant constructs. These results are consistent with decreased CHL dynamics in the K49W and in the K49A mutants.

To study the conformational changes and allosteric effect, adaptive sampling was carried out to expand the sampling along the reaction pathway of the K49A mutation for another 48 μs in 24 replicas (Supplementary Figure S10C). The impact of the mutation on the allosteric mechanism of conformational changes in StarD4 is evident in the change of contact probability between residues in β1 and H4 as shown in Table 3.

The resulting structures, compared in Figure 8, show how the K49W and K49A mutations changed the mode of interaction between β1 and H4. In the WT StarD4, K49 establishes both hydrophobic and hydrophilic contact with A207, S208 and T204, and its translocation opens space for the formation of a V51-A207 contact pair (Figures 8A, B). With the K49W mutation, the large W49 forms stable contacts with both S208 and T204 but hinders sterically the V51-A207 interaction (Figure 8C; Table 3). The K49A mutation, on the other hand, resulted in a decreased contact rate with both S208 and T204, eliminated the hydrophilic interactions, and destabilized the V51-A207 contact (Figure 8D;Table 3).

The interactions between β1 and H4 are weakened overall in both mutants but the interaction modes differ. The probability of contact is decreasing gradually from WT to K49W to K49A, and the electrostatic interaction is also gradually removed in the same order. Consistent with the weakening of interaction and the decrease of CHL dynamics, the NbIT results showed a decrease in coordination information (CI values) in the allosteric network between the β1, H4 and β9 motifs the K49W mutant (Supplementary Table S6), and a loss of coordination effect from β1 to the CHL site in the K49A mutant (Supplementary Table S7)–as expected from our results from the analysis of WT StarD4.

The structure of mouse apo-StarD4 was taken from the X-ray crystal structure (PDB ID: 1jss), which includes residues 24 to 222 of the protein (Romanowski et al., 2002). Lys223 and Ala224 were introduced using Modeller 9.23 software (Sali, 1995), resulting in the conformation of mouse StarD4 22–224 segment with N-terminus acetylation. StarD4 complexes with cholesterol (CHL) were obtained by docking the ligand in the hydrophobic pocket of StarD4 using Schrodinger’s Induced Fit Docking protocol (Sherman et al., 2006a; Sherman et al., 2006b). The common features shared by the top 15 docking results are: the interaction of the cholesterol hydroxyl head that forms 2 or 3 hydrogen bonds with the sidechains of Ser136 & Ser147 (Figure 1A) and with the backbone carbonyl of Cys148, and the sequestering of the cholesterol hydrocarbon tail inside the hydrophobic pocket, away from the water environment. This top docking mode was used as the initial structure for the MD simulations of the StarD4-CHL complex solvated in 0.15M K⁺Cl⁻ ionic aqueous solution (∼32,600 atoms).

The cholesterol-bound StarD4 system was equilibrated first in MD simulations with NAMD version 2.12 (Phillips et al., 2005) using a multi-stage protocol. In the first stage the backbone of StarD4 and the heavy atoms of the ligand were harmonically constrained and gradually released in three steps of 0.5ns each, with restraining force constants changing from 1 kcal/(mol·Å2) to 0.5, and 0.1, respectively. This stage was followed by 6ns unbiased MD simulation using NAMD. In the third (production) stage, the velocities of all the atoms were reset, and the system was simulated with openMM software (Eastman et al., 2017) in 12 independent replicates, each for 33.6μs, resulting in a cumulative simulation time of 403μs (>0.4 milliseconds). The OpenMM simulations were conducted in NPT ensemble (T = 310K, p = 1 atm) using a 4fs integration time-step and a Monte Carlo Membrane Barostat.

Mutations were introduced into the StarD4-CHL complex with the free-energy perturbation protocol in NAMD version 2.12 (Phillips et al., 2005). Accordingly, hybrid StarD4 models containing overlapped residues (e.g., overlapped K49 and A49 to introduce the K49A mutation) were constructed using representative conformations obtained in the most popular Ser-binding state 1 and Trp-binding state 3. Then, the mutations are introduced by gradual annihilation of the interactions between the WT residue and its surrounding, achieved by decreasing the coupling parameters from 1 to 0 in 10 steps, and concomitantly increasing the interactions between the mutant residue and its surrounding from 0 to 1 in 10 steps. The simulation was run for 5ns per window. After the last window, the annihilated original residues are removed from the hybrid StarD4 models, and the mutant StarD4 is equilibrated for another 10ns. For production runs, 24 replicas are made for each system and sampled in long unbiased simulation with openMM (Eastman et al., 2017).

Rare events in MD trajectory data were detected as described with our RED protocol (Plante and Weinstein, 2021) which utilizes a Non-negative Matrix Factorization (NMF) algorithm implemented in the Scikit-Learn python package (Fabian Pedregosa et al., 2011). The unsupervised machine learning algorithm learns to decompose a trajectory into a set of components that represent structural motifs that move simultaneously, and reports their appearance as dominant structural characteristics during a particular period in the course of the trajectory. NMF is a machine learning technique used to decompose high-dimensional non-negative data. It has been widely applied in computational biology research, where it has proven to be a powerful tool for tasks such as the molecular pattern discovery, class comparison and prediction, and functional characterization of genes as described in Plante and Weinstein (2021) and Lee and Seung (1999). Briefly summarized, the event detection method in the RED protocol uses the NMF algorithm to learn a sparse, parts-based representation of the data (Lee and Seung, 1999). For trajectory analyses the input array (I) is a contact matrix constructed as described below from the data in each frame of the trajectory. By construction, this array encodes the information about the protein structure as it evolves over the trajectory time. The NMF algorithm decomposes the I array into a “spatial” array and a “temporal” array, which together are responsible for the integrated detection of the temporally defined rare events that provide the mechanistic information. Thus, the columns of the “spatial” output matrix produced by the NMF analysis are termed “components”, as each of them represents a set of conformational features that change simultaneously even if they are in different structural motifs. As the determinant features of a conformational state evolve over trajectory time, the different components containing them dominate the structural characteristics of the molecule at different times. Thus, the corresponding temporal array represents the contribution of each component along the trajectory, and thus identifies the time period where a particular component dominates the structural characteristics of the protein (Plante and Weinstein, 2021). A simple illustration provided in the full description of the method (Plante and Weinstein, 2021) refers to the analysis of a conformational transition in the unfolding of a fully folded alpha helical segment. The analysis results in one component that is the folded structure, and a second one that is the unfolded structure. Each frame of the simulation trajectory will be a mixing of these archetypes, with the folded component dominating the trajectory frames at the time segments before the transition, and the unfolded component dominating those after the transition.

Building the time-evolution contact map for the input array ( I ): A contact map is constructed between every pair of residues for each frame in the trajectory, with contact marked as 1 if any atom of one residue is <3.5Å away from any atom in another residue. For the 201 residues of StarD4 this yields a tensor of 0 and 1 values with dimensions (t, 201,201), where t represents the number of frames. Subsequently, the tensor is rearranged to a matrix of size 201² columns and t rows, where 201² is the total number of residue pairs. The matrix columns, representing all the residue pairs, are further trimmed by excluding the contact pairs that remain invariant throughout the entire simulation, which yields a I matrix with n_pairs = 3,880 columns. Each column represents a contact pair that has altered its contact state at least once during the entire simulation, encapsulating protein dynamics information for further analysis.

To focus on the conformational changes around the time of cholesterol relocation between its binding modes in the hydrophobic pocket, we excerpted 6 trajectory segments from the 6 independent runs as indicated by the time intervals (traj0, from 6.5 to 14.5μs; traj1, 0.1–4.1μs; traj3, 0.1–8.1μs; traj4, 0.1–6.5μs; traj6, 0.1–4.1μs; traj7, 0.1–4.1μs; 34.4μs total). In order to obtain results that are generalizable among trajectories (which decreases the sensitivity to incidental structural fluctuations and benefits the reproducibility of RED analysis), the contact matrices of the 6 segments are inputted at the same time as the training data for the NMF algorithm. This input of multiple trajectories is constructed by concatenating the data matrices along the time-axis. This contact map is recorded every 0.8ns, which yields an I matrix with t = 43,000 rows. Smoothing is applied to the contact matrix across the time dimension by averaging the array in a sliding window of 20 nanoseconds (25 frames).

Details of the protocol: The algorithm decomposes the input contact matrix described above into multiple components, where each component has a structural state described by a “spatial” row array with length n_pairs (here n_pairs = 3,880), and a corresponding “temporal” column array of length t (here t = 43,000). The number of components (c) is a super-parameter chosen by the user according to the number of independent processes in the system and was set to be c = 5 to study generalizable features shared in independent runs. Thus, the ensemble of the “temporal” column arrays is a temporal matrix of size (t, c) in this study, and the ensemble of the “spatial” row arrays is a spatial matrix of size (c, n_pairs = 3,880), where the multiplication of the temporal matrix and the spatial matrix is fitted to resemble the time-evolution contact map data matrix. Each row in the c = 5 rows of components in the spatial matrix is an array that encodes the determinant conformational features that evolve concurrently over time. The corresponding column in the temporal matrix represents the contribution of these features in determining the conformational state of the protein at each time. It and identifies the trajectory times at which this particular set of features dominates the structural characteristics of the protein by showing it to have the highest “weight” in this period of trajectory time.

Interpretation of the results: As shown in Figures 2E–G, the spatial arrays scatter plots show that there are groups of pairs that form a horizonal line, suggesting that these pairs are constantly contributing similarly to every component. These groups of pairs that are salient in every component by being in contact throughout the StarD4 simulation, and contribute similarly to each component, are used for normalization of the spatial arrays. Thus, each spatial array is divided by a normalization factor that equals to the mean value of the constructive residue pairs contributing to this particular component, thereby placing the contact features of these constitutive residue pairs at 1 in Figures 2E–G; Supplementary Figure S2B To preserve the relative weight between components, each corresponding temporal array is multiplied by the same normalization factor derived from the spatial array of this particular component, so that the multiplication result of spatial and temporal matrix remains the same and resembles the input contact matrix.

This normalization allows for direct comparison between components. Subtracting the spatial array of comp A from comp B, as represented in (Supplementary Figure S3), cancels out the contributions of the constitutive residue pairs so that pairs with positive contributions are those newly established in the event corresponding to the A→B transition, whereas the negatively contributing pairs are those breaking contact during this A→B event (Supplementary Figure S3). These pairs are salient in the transition as they contribute to the structural differences between components and are termed structure differentiating contact pairs (SDCPs). The SDCPs that are more salient than 4 times the standard deviation, are summarized in Supplementary Table S1 and are grouped based on the motif movement they represent.

The MD trajectory frames were projected onto a space of 7 collective variables described in Supplementary Table S2, with the “time-structure based independent component analysis” (tICA) method. The tICA method identifies the slowest reaction coordinates of a system in the dynamic dataset (Molgedey and Schuster, 1994; Naritomi and Fuchigami, 2011; Perez-Hernandez et al., 2013) by solving for the generalized eigenvectors and eigenvalues of the time-lagged covariance matrix: C x τ V i = λ i C x 0 V i in which the C x τ is the time-lagged covariance matrix C x τ = X t X T t + τ t between collective variables, C x 0 = X t X T t t is the regular covariance matrix, X t is the data vector at time t , τ is the lag-time, V i is the i^th tICA eigenvector, and λ i is the i^th tICA eigenvalue. In this study, the stride of the trajectories is 0.4 ns/frame, and the lag-time τ for the covariance matrix construction is 10ns (25 frames), thus ignoring any fluctuations faster than 10ns which are not relevant to the slow processes in the system. The first 1μs in every independent run are discarded to equilibrate the systems.

1. Construction of the Markov State Model

Markov State Model (MSM) is a powerful tool to study the kinetics of equilibrium state-transition processes in protein function or activation pathways (Beauchamp et al., 2012; Kohlhoff et al., 2014; Shukla et al., 2014). The implementation and validation of MSMs are well documented (Noe and andFischer, 2008; Pande et al., 2010; Prinz et al., 2011). To construct MSMs, the 3D conformational space of the first three tIC vectors was discretized into microstates. The transition count matrixes between the microstates were then recorded in a transition count matrix, which is then symmetrized by averaging with its transpose matrix in order to satisfy detailed balance and local equilibrium (Beauchamp et al., 2011; Prinz et al., 2011). Finally, the transition probability matrixes (TPMs) were built by normalizing the transition count matrix by the transitions departing from the same microstate.

2. Selection of parameters for markov model construction

The choice of super parameters: the lag-time of the transition, and the discretization of the conformational space (i.e., number of microstates), are optimized to ensure the quality of MSM (Beauchamp et al., 2011; Prinz et al., 2011). The lag-time is optimized with implied timescale test: τ i = τ ′ ln ⁡ ⁡ λ i in which τ ′ is the lag-time used for building the transition matrix, λ i is the i^th eigenvector of the TPM, and τ i is the implied timescale corresponding to the i^th relaxation mode of the system. In this study, the lag-time of 120 ns defined as shown in Supplementary Figure S5A ensures the best available Markovian behavior with the minimal loss of data. The best discretization of the conformational space is chosen by comparing the resulted MSM with the generalized matrix Rayleigh quotient (GMRQ) method (Beauchamp et al., 2011; McGibbon and andPande, 2015), where the MSM is trained on a randomly picked training set comprised of half the set of trajectories, and then cross-validated on the test set which includes the other half of the trajectory set. The method tests how well the eigenvectors of the training set diagonalize the TPM of the test set, scored as: R = T r V T S T V V T S V − 1 in which T r denotes trace of the inner matrix product, V contains the first n eigenvectors of the training set, S is the diagonal matrix composed of equilibrium populations of the test set (i.e., a diagonal matrix composed by the first eigenvector of the test set TPM), and T is the test set TPM. The best MSM parameters yield the highest GMRQ score. In this study, the 200 k-means microstates were chosen so as to provide a good balance between the consistency shown by GMRQ result (Supplementary Figure S5B) and the coverage of the 3D conformational space.

3. Interpretation of MSM, and macrostate assignment

To extract information from the MSM, TPM is decomposed into its eigenvalues and eigenvectors. The largest eigenvalue is 1, and its corresponding first eigenvector represents the equilibrium state population of the system. The remaining eigenvectors represent the population flow between microstates. The second largest eigenvalue corresponds to the eigenvector representing the population flow contributing to the slowest relaxation process of the system. In turn, the next largest eigenvalues and eigenvector correspond to the next slowest relaxation process, and so on. Here, the 200 microstates were grouped into 9 macrostates according to their kinetics similarities learnt from MSM, using the using the Robust Perron Cluster Analysis (PCCA+) algorithm (Deuflhard and Weber, 2005).

The transition paths among the 9 macrostates were obtained from the TPT analysis (Berezhkovskii et al., 2009), which is based on the flux matrix between states whose components are defined as: J i j = π i 1 − P t a r g e t , i T i j P t a r g e t , j where π i is the equilibrium population of state I, T i j is the transition probability from state i to state j, and P t a r g e t , i is the probability of visiting the target state i as the first state during the transition before visiting other states. The most probable transition pathways between specific states was obtained by iteratively finding the next pathway with the highest flux in the flux matrix using a graph theory algorithm, the Dijkstra algorithm (Dijkstra, 1959), implemented in the MSM builder software package (Beauchamp et al., 2011).

1. NbIT analysis focused on multiple motifs in the StarD4-CHL complex

The information theory-based NbIT analysis was used to measure the coordination between distanced motifs and to quantify the allosteric effect (LeVine et al., 2014; LeVine and Weinstein, 2014). For the calculation of the coordination information we chose the following structural motifs: CHLsite: S136 S147 W171 R92 Y117; β1: res R46 V47 A48 K49 K50 V51 K52; β2β3: res R58 K59 P60 Y67 L68 Y69; β2β3loop: E62 E63 F64 N65; H4head: Q199 S200 A201 D203 T204 A207; Ω1loop: L124 N125 I126; β9: D192 R194 G195; β7β8loop-nearβ9: V162 R163 G164; β7β8loop-mid: T157 R158 P159 E160; β7β8loop-nearβ6: E153 W154 S155 E156; H4tail: R218 K219 G220 L221; β3tail: G73 V74 M75 D76; β8β9loop: S179 P180 S181 Q182. The allosteric effect was then quantified from the coordination information as described in ref 32.

Briefly, the first step in the calculation of the coordination information was the alignment of the trajectories to a reference structure–we used the StarD4 crystal structure (PDB ID: 1JSS). The alignment of snapshot (per 0.4ns) was carried out on the alpha carbons of residues in β3-9 and H4, excluding the motifs have been found to undergo conformational changes (70–76, 79–86, 103–108, 113–118, 131–140, 146–150, 170–175, 184–189, 211–220). With the aligned trajectory, the entropy in each motif was calculated analytically through the differential entropy: H X = 1 2 ln ⁡ ⁡ 2 π e C X in which C X is the covariance matrix describing all variables in X, and X is the displacement array describing the position of every atom in the motif by its distance to the ref conformation. Then, the Total Correlation (TC) describing the total amount of information share in a set of motifs was calculated as: T C X 1 , … , X N = ∑ i N H X i − H X 1 , … , X N .

And the Coordination Information (CI) describing the amount of information shared in the receptors that is also shared with another motif that works as the transmitter: C I X 1 , … , X N , X m = T C X 1 , … , X N − T C X 1 , … , X N   |   X m

Here T C X 1 , … , X N | X m is the conditional total correlation between X 1 , … , X N conditioning on X m : T C X 1 , … , X N   |   X m = ∑ i N H X i , X m − H X m − H X 1 , … , X N , X m − H X m .

The Coordination Information was then normalized by the Total Correlation of the receiver, to represent the portion of dynamics of the receptor that is allosterically coupled with the transmitter. C I ¯ X 1 , … , X N , X m = C I X 1 , … , X N , X m T C X 1 , … , X N ∗ 100 % .

In Table 2 and the Supplementary Figures S3, S5, S6, S7, the Total Correlation (TC) within each motif is displayed along the diagonal, while the Normalized Coordination Information is presented in the off-diagonal elements, with residues on the top (columns) acting as the Transmitter and residues on the left (rows) being the coordinated Receiver.

3. Coordination channel analysis based on mutual coordination information

To define the allosteric channels that mediate the coordination information (Figure 6A), we calculated the mutual coordination information that measures the amount of coordination information shared between the receiver and the transmitter motifs that is also shared with another structural element: M C I X 1 , … , X N , X m , X n = C I X 1 , … , X N , X m + C I X 1 , … , X N , X n − C I X 1 , … , X N , X m X n

Here X m X n is the union set constituted by the residues in the Transmitter X m and the channel X n . The mutual coordination information was then normalized by the Coordination Information between the receiver and the transmitter to obtain the NMCI values used in the determination of the allosteric pathway in Figure 6B N M C I X 1 , … , X N , X m , X n = M C I X 1 , … , X N , X m , X n C I X 1 , … , X N , X m .