Breakage of Hydrophobic Contacts Limits the Rate of Passive Lipid Exchange Between Membranes

The maintenance of heterogeneous lipid compositions among cellular membranes is key to biological function. Yet, even the simplest process that could be responsible for maintaining proper lipid distributions, passive lipid exchange of individual molecules between membranes, has eluded a detailed understanding, due in part to inconsistencies between experimental findings and molecular simulations. We resolve these discrepancies by discovering the reaction coordinate for passive lipid exchange, which enables a complete biophysical characterization of the rate limiting step for lipid exchange. Our approach to identify the reaction coordinate capitalizes on our ability to harvest over 1,000 unbiased trajectories of lipid insertion, an elementary step of passive lipid transport, using all-atom and coarse-grained molecular dynamics simulations. We find that the reaction coordinate measures the formation and breakage of hydrophobic contacts between the membrane and exchanging lipid. Consistent with experiments, free energy profiles as a function of our reaction coordinate exhibit a substantial barrier for insertion. In contrast, lipid insertion was predicted to be a barrier-less process by previous computational studies, which incorrectly presumed the reaction coordinate to be the displacement of the exchanging lipid from the membrane. Utilizing our newfound knowledge of the reaction coordinate, we formulate an expression for the lipid exchange rate to enable a quantitative comparison with experiments. Overall, our results indicate that the breakage of hydrophobic contacts is rate limiting for passive lipid exchange and provide a foundation to understand the catalytic function of lipid transfer proteins. Graphical TOC Entry


INTRODUCTION
Lipid membranes are essential to the structural integrity of cells. They form fluid boundaries that organize and compartmentalize cellular functions within organelles. Each organelle requires a unique membrane composition for its proper function. 1,2 To maintain such varied compositions, lipids are heterogeneously trafficked between cellular membranes through vesicular or non-vesicular transport. Vesicular transport plays a major role in trafficking lipids and proteins between organelles in the secretory pathway. 3,4 Organelles that are not connected by vesicular transport machinery rely on non-vesicular mechanisms to receive and export lipids. Even for organelles in the secretory pathway, non-vesicular transport mechanisms provide an additional way to more rapidly exchange lipids, for example, to swiftly alter membrane compositions in response to environmental changes. 5,6 Despite its cellular importance, non-vesicular transport mechanisms in vivo have not been fully characterized. Non-vesicular transport predominantly involves the movement of individual lipids between membranes. Monomeric lipid transfer between membranes may occur passively, in which a lipid desorbs and freely diffuses to another membrane. 5 Alternatively, lipid transfer proteins may facilitate monomeric lipid exchange by enclosing lipids within their hydrophobic interiors dur-ing transport. 6,7 With half-times on the order of hours, passive lipid exchange is too slow to fully account for lipid transport in vivo, 5 indicating that lipid transfer proteins are likely catalysts of lipid exchange. Nevertheless, a mechanistic understanding of passive lipid exchange can inform our knowledge of how lipid transfer proteins increase the rate of lipid transport.
Molecular simulations are well suited to explore the microscopic dynamics and to identify the rate limiting step of passive lipid exchange. Quantifying the free energy barrier associated with the rate limiting step, however, requires knowledge of the reaction coordinate, which characterizes the collective motion of molecules that advances a transition. [8][9][10] Previous computational work [11][12][13][14][15] on lipid transport has presumed that a lipid's displacement normal to the membrane is the reaction coordinate ( Figure 1A) and has yielded results in conflict with experimental findings. [16][17][18][19] Here we show that the reaction coordinate for passive lipid exchange is indeed more subtle than a simple distance measurement. The reaction coordinate characterizes the creation (or disruption) of a locally hydrophobic environment around the incoming (or outgoing) lipid ( Figure  1B). This realization resolves qualitative (but not quantitative) discrepancies between simulation and experiment and suggests that the breakage of hydrophobic contacts between a lipid and membrane limits the rate of passive lipid transport.

B
✘ Incorrectly predicted to be a TS based on lipid's displacement ✔ Correctly predicted to be a TS based on hydrophobic contacts A Figure 1: A lipid's displacement normal to the membrane has conventionally been presumed to be the reaction coordinate for passive lipid exchange; however, the lipid's displacement does not reliably identify transition state (TS) configurations. Although the lipid's displacement (black arrow) in configuration A is consistent with values observed at the transition state, this configuration is not a transition state. By contrast, configuration B is correctly predicted to be a transition state based on the extent of hydrophobic contact between the lipid and membrane (highlighted in yellow and cyan).
Experimental and Computational Background. Numerous in vitro studies of passive monomeric lipid exchange between membranes have demonstrated that it is a first-order process and that the rate of lipid exchange strongly correlates with a lipid's solubility. [16][17][18][20][21][22][23][24][25][26] Based on these observations, the most widely accepted mechanism is characterized by aqueous diffusion. [16][17][18]20,[23][24][25][26][27][28][29][30] The diffusive mechanism involves lipid desorption, which is rate limiting, followed by diffusion through solvent and insertion into another membrane. Based on experimentally measured activation energies, calculated activation free energies for lipid exchange exceed free energies for transferring a lipid from water to a membrane, indicating that there is a barrier for lipid desorption and insertion. [16][17][18][19] At this barrier, the desorbing lipid is hypothesized to only have the terminal carbons of its tails left within the membrane. 16,27 Thus, the activation free energy required to form the transition state has been attributed to the creation of a cavity in the membrane due to partial removal of a lipid and another cavity in the solvent to accommodate that lipid. 16,31 However, because the experimental methods currently used to study lipid exchange have coarse temporal and spatial resolution, molecular features of the transient transition state can only be hypothesized. Molecular dynamics (MD) simulations offer an attractive means to test these hypotheses since the necessary time and length scales are accessible. Additionally, a complete free energy profile can be obtained from MD simulations. By calculating the free energy as a function of the reaction coordinate, which describes the system's dynamics during transitions between stable states, [8][9][10] free energies of activation can in principle be accurately quantified. If, however, the free energy is computed as a function of an order parameter that is not the reaction coordinate, then the apparent barrier generally underestimates the rate-determining free energy of activation. Thus, identifying the reaction coordinate and corresponding free energy profile for passive lipid exchange can provide fundamental insights into the physical processes and work required to maintain heterogeneous cell membrane compositions.
Most previous computational studies have focused on obtaining a full free energy profile for lipid desorption and insertion. [11][12][13][14][15] Traditionally, free energy profiles have been computed as a function of a lipid's displacement normal to a bilayer measured from the lipid's phosphate group to the bilayer's center-of-mass (COM) ( Figure 1A). [11][12][13][14][15] These free energy profiles lack a barrier for insertion, 11-15 seemingly in conflict with experimental results. [16][17][18][19] If the COM displacement is not the reaction coordinate for lipid exchange, the kinetically relevant barrier may not be resolved in these free energy profiles; instead, a barrier may exist along a different degree of freedom that is the reaction coordinate. Consistent with this idea, Vermaas and Tajkhorshid's MD study of lipid insertion indicated that the COM displacement is not sufficient to fully describe the microscopic dynamics of lipid insertion. They demonstrated that after the lipid associates with a bilayer, each tail of the lipid enters the bilayer successively to complete the insertion processes. The observation of splayed lipid intermediates during insertion, which are indistinguishable from other configurations based on the lipid's COM displacement, suggests that other degrees of freedom need to be considered to construct an accurate reaction coordinate for lipid exchange. 32 Our Approach. The discrepancy between experimental [16][17][18][19] and computational [11][12][13][14][15] reports about a barrier for lipid exchange, together with the evidence suggesting that the lipid's COM displacement is a poor reaction coordinate from a MD study of lipid insertion 32 prompt two questions: (1) What is the reaction coordinate for lipid exchange? (2) What, if any, activation free energy barrier impedes the process of lipid insertion? In this article, we aim to answer these questions using molecular simulation. In doing so, we identify the reaction coordinate for passive lipid exchange in all-atom and coarse-grained lipid models. Knowledge of the reaction coordinate allows us to elucidate key biophysical details of the transition state ensemble and properly assess the free energetic cost of lipid transport.
Rather than driving the system along a presumed reaction coordinate, we instead harvest natural, unbiased trajectories in which a lipid spontaneously inserts into a membrane. Statistical analysis of this ensemble of dynamical pathways reveals the key collective motions required for lipid transport. Firstly, we find that lipid insertion is a barrier crossing event and occurs via three different pathways, distinguished by various splayed lipid intermediates. Secondly, we find that the reaction coordinate characterizes the formation and breakage of hydrophobic lipid-membrane contacts. The lipid's displacement normal to the bilayer, even formulated to distinguish splayed configurations, is not the reaction coordinate and obfuscates the barrier for insertion. Consistent with previous experimental results, [16][17][18][19] free energy profiles as a function of our reaction coordinate display a barrier for insertion and yield insertion rates calculated from Kramers theory that agree with those obtained directly from simulations. Finally, using our newfound reaction coordinate, we formulate a Smoluchoski equation for the lipid exchange rate to directly compare our simulation results with experiments. Overall, our results demonstrate that the rate limiting step for passive lipid exchange is the breakage of hydrophobic contacts, suggesting that lipid transfer proteins may catalyze lipid transport in part by lowering the associated activation free energy.

Molecular Dynamics
Simulations. Spontaneous desorption of a lipid from a membrane is a very slow process occurring over minutes to hours, which is well beyond the timescale accessible in MD simulations. Most previous work has addressed this problem by introducing an external bias that allows rare configurations otherwise inaccessible in MD simulations to be readily sampled. While this approach generates configurations plausibly found along a lipid desorption trajectory, it can fail to reveal the natural, unbiased route of lipid desorption. Our approach instead exploits a fundamental statistical property of microscopic dynamics, namely its time reversibility. Natural desorption trajectories are simply the time reverse of spontaneous lipid insertion trajectories. As shown by Vermaas and Tajkhorshid, 32 the latter are straightforward to generate in an unbiased way since lipid insertion is a rapid process. Thus, to gain insights into the dynamics of lipid exchange, we harvested trajectories of lipid insertion into a bilayer of 128 lipids using both all-atom and coarse-grained MD simulations. All simulations were performed in an isothermal-isobaric (NPT) ensemble using GROMACS 5. 33 The temperature was maintained at 320 K, ensuring that the bilayers were in the liquid crystalline phase. The pressure was maintained at 1 bar using semi-isotropic pressure coupling to allow the z dimension, which is perpendicular to the bilayer, to fluctuate separately of x and y, ensuring tensionless bilayers.
Coarse-Grained Systems. A large number of lipid insertion trajectories are required to make conclusions with high statistical accuracy. To obtain 1,000 lipid insertion trajectories in a computationally tractable manner, we performed MD simulations using the coarse-grained MARTINI force field. 34 All coarse-grained simulations used dilauroylphosphatidylcholine (DLPC) to compare to all-atom simulations of dimyristoylphosphatidylcholine (DMPC). Because MARTINI maps roughly four heavy atoms to a single coarse-grained bead, DMPC, which has 14 carbon atoms per tail, is best represented by the MARTINI model for DLPC, which has 3 beads per tail.
First, a bilayer surrounded by 3 nm thick slabs of standard MARTINI water was built using INSANE. 35 Prior to the addition of a tagged lipid into the solvent, a 50 ns simulation was performed to allow the bilayer's structure to fully equilibrate, as monitored by the area per lipid ( Figure S1). Next, 1,000 replicate systems with a free lipid were built by inserting a tagged lipid at a random location in the solvent around the equilibrated bilayer such that the tagged lipid's COM and bilayer's COM are separated in z by at least 3.2 nm. Production runs of 1 µs were performed, and lipid insertion occurred in all replicates during this time. Reported times for MARTINI simulations are not scaled by a factor of 4, as done in other work to account for the erroneously fast solvent diffusion in MARTINI as compared to real water. 34 Complete simulation details are provided in the Supporting Information (SI).
All-Atom Systems. Additionally, we harvested 10 all-atom lipid insertion trajectories to compare with our results using MARTINI. The CHARMM36 force field 36 was used in combination with the CHARMM TIP3P water model 37 since it accurately reproduces many experimental observables, including the volume and area per lipid, bilayer thickness, lipid lateral diffusion coefficient, and neutron density profiles, for a liquid crystalline DMPC bilayer. 38,39 We used a simulation protocol similar to that described above for MARTINI.
First, a bilayer surrounded by 3 nm thick slabs of solvent was built using the CHARMM-GUI Membrane Builder. 40,41 A 50 ns production run was performed to allow the bilayer to fully equilibrate ( Figure S1). Next, 10 replicate systems with a tagged lipid in solution were built as for the MARTINI systems. Each replicate was simulated in increments of 100 ns until the tagged lipid inserted into the bilayer. Complete simulation details are provided in the SI.

Characterization of Transition Paths.
From the harvested lipid insertion trajectories, we identified transition paths, trajectory segments that connect "reactant" and "product" states A and B. In state A, the tagged lipid is fully surrounded by solvent; in state B, the tagged lipid is fully within the bilayer. States A and B are characterized by the displacements d lip , d sn1 , and d sn2 shown in Figure 2A. d lip is the displacement in z from the COM of the tagged lipid to the COM of the closest leaflet. Similarly, d sn1 is the displacement in z from the terminal carbon of the sn1 tail of the tagged lipid to the COM of the closest leaflet, and d sn2 is the analogous distance for the sn2 tail. Based  on distributions of these three distances obtained from coarse-grained and all-atom MD simulations ( Figure  S2), state A is defined by d lip > 24Å, which ensures that any configurations with the tagged lipid adsorbed onto the surface of the bilayer are not included. State B is defined by d sn1 < −3Å and d sn2 < −3Å, which ensures that the both tails of the tagged lipid are inserted into the bilayer. In addition to d lip , d sn1 , and d sn2 , we evaluated over 50 order parameters as putative reaction coordinates for lipid exchange. Hydrophobic contacts between the tagged lipid and bilayer are judged according to the distance d CC between a hydrophobic carbon of the tagged lipid and a hydrophobic carbon of the closest membrane leaflet, where i indexes the many pairs of such atoms. For MARTINI lipids, hydrophobic carbons include all tail beads. For CHARMM36 lipids, hydrophobic carbons include atoms C23 -C214 and C33 -C314.
More specifically, hydrophobic contacts are measured by the order parameters min(d CC ) and n CC ( Figure  2B). min(d CC ) = min i d (i) CC is the minimum distance between a hydrophobic carbon of the tagged lipid and a hydrophobic carbon of the closest leaflet. n CC is the number of close contacts between hydrophobic carbons of the tagged lipid and hydrophobic carbons of the closest leaflet. The ith pair of hydrophobic carbons was counted as a close contact if d These cutoff values encompass approximately two water solvation shells around a hydrophobic carbon of a lipid in solution and also two carbon solvation shells around a hydrophobic carbon of a lipid in a bilayer. We found smaller cutoffs to be insufficient for fully characterizing the tagged lipid's hydrophobic environment and, thus, also insufficient for constructing the reaction coordinate. Complete descriptions of all other order parameters are provided in the SI. The MDAnalysis Python library 42 was used to analyze all trajectories for each order parameter.
Free Energy Calculations. Based on the dynamics observed along transition paths, lipid insertion is a barrier-crossing process. To identify the physical origin of this barrier, we calculated free energy surfaces as a function of different order parameters. For the MARTINI system, we calculated a free energy surface as a function of d sn1 and d sn2 . We also calculated free energy surfaces as a function of min(d CC ) and n CC for both the MARTINI and CHARMM36 systems. To obtain these free energy surfaces, we performed umbrella sampling simulations 43 using the PLUMED 2 patch 44 for GROMACS and used the weighted histogram analysis method (WHAM). 45 Complete details about the free energy calculations are provided in the SI.
Committor Analysis. To identify the reaction coordinate for passive lipid exchange and determine if transition states are sampled in umbrella sampling simulations, we characterized configurations according to their tendency to proceed to state B using committor analysis. [8][9][10] The committor, p B , is the probability that a configuration will reach state B prior to state A when its momenta are chosen randomly from a Maxwell-Boltzmann distribution. By construction, the committor distinguishes transition states, which have p B = 0.5, from stable state A and B configurations, which have p B = 0 and 1, respectively. Thus, p B is the true reaction coordinate. Through committor analysis of configurations found along transition paths, we identified order parameters that are strongly correlated with the committor and, thus, can be used as approximate reaction coordinates. These order parameters have the advantage of being more physically descriptive and, thus, more easily interpreted than the committor. Henceforth we refer to order parameters that correlate strongly with p B as the reaction coordinate. We calculated commit-tor values for 98,094 MARTINI and 138 CHARMM36 configurations sampled along transition paths. Additionally, we calculated the committor for 500 MAR-TINI and 100 CHARMM36 configurations from each umbrella sampling simulation. For each MARTINI configuration, the outcome of 50 trajectories, each 3 ns long and initialized with random velocities sampled from a Maxwell-Boltzmann distribution, were used to calculate its committor value. For each CHARMM36 configuration, the outcome of 20 trajectories, each 12 ns long, were used to calculate its committor value.
Rate Calculations. To further assess how well the dynamics of lipid transport are captured by monitoring lipid-membrane hydrophobic contacts, we calculated the rate constant for lipid insertion, k ins , using Kramers theory 46 coupled with thermodynamic information about hydrophobic contacts. The resulting value was compared with the insertion rate calculated from the mean first passage time in MD simulations. Based on Kramers theory, where D is the diffusion coefficient along a reaction coordinate x for lipid insertion, x B is the value of x in state B, L is the maximal value of x in state A sampled in our simulations, U (x) is the effective interaction potential biasing the dynamics of x, and β = (k B T ) −1 is the inverse of Boltzmann's constant, k B , multiplied by temperature. We have taken x to be min(d CC ) since min(d CC ) describes the spatial motion of the lipid during insertion in addition to hydrophobic contact formation. For the same reason we set D to be the diffusion coefficient of a freely diffusing lipid in solution. Note that min(d CC ) alone is not the reaction coordinate for lipid exchange (a linear combination of min(d CC ) and n CC is the reaction coordinate), but it is more simply interpretable to utilize a single order parameter for these calculations. U (x) is taken to be the free energy profile along min(d CC ). Details about how U (x) is obtained from the free energy surface that depends jointly on min(d CC ) and n CC are provided in the SI. We evaluated Eq. 1 with numerical integration using the trapezoidal rule. The diffusion coefficient was calculated from an additional MD simulation of a single lipid solvated in a cubic box with the same area and simulation protocol as the bilayer systems, but with isotropic pressure coupling. D was calculated according to Einstein's relation from the mean squared displacement of the lipid's COM obtained from 1 µs and 100 ns trajectories of a MARTINI and CHARMM36 lipid, respectively.

RESULTS AND DISCUSSION
Lipid Insertion Is a Barrier Crossing Process that Occurs via Multiple Pathways.
We first investigated the dynamics of lipid insertion by harvesting 1,000 MARTINI DLPC and 10 CHARMM36 DMPC insertion trajectories from MD simulations. A lipid insertion event is classified by a transition from state A, in which the tagged lipid is fully in the solvent, to state B, in which the tagged lipid resides within the bilayer. State A is distinguished by d lip , which measures the COM displacement of the tagged lipid along the bilayer's normal; state B is distinguished by d sn1 and d sn2 , which measure the displacement of each tail of the tagged lipid along the bilayer's normal ( Figure  2A). Precise definitions are given above in the Methods section. The fact that state B configurations cannot be reliably identified by d lip alone suggests that the COM displacement is not the reaction coordinate for passive lipid exchange.
Free energy profiles as a function of the lipid's displacement along the bilayer normal obtained from previous computational studies [11][12][13][14][15] give the impression that lipid insertion is a barrier-less process. If that were true, insertion should occur immediately once the tagged lipid reaches the bilayer. However, as seen in snapshots from a MARTINI and a CHARMM36 trajectory shown in Figure 3 and Movies S1 and S2, the tagged lipid repeatedly arrives at the bilayer and adheres to its surface without inserting. Instead, it detaches from the bilayer's surface and returns to the solvent. In typical trajectories, many such unproductive encounters occur before the lipid inserts into the bilayer. Indeed, as seen in time traces of d lip in Figure 4A, many adsorption events commonly precede insertion. Similar adsorption events have been observed during simulations of 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) insertion into a bilayer. 32 When lipid insertion eventually occurs in these trajectories, it does so suddenly. Characteristic of barriercrossing dynamics, d lip , d sn1 , and d sn2 change sharply from values of state A to those of state B ( Figure  4A). The fact that the transition times are much faster than the inverse rate constant for lipid insertion, 1/k ins , points to a substantial free energy barrier for insertion. Free energy profiles as a function of the lipid's displacement simply do not resolve this barrier. [11][12][13][14][15] Thus, the barrier must exist along a different degree of freedom that captures other important features of the dynamics.
In fact, lipid insertion occurs via three different pathways which cannot be differentiated by d lip . Each pathway is characterized by a distinct lipid configuration, which is distinguished by d sn1 and d sn2 , near the bilayer's surface: (1) In the sliding pathway, the two tails enter the bilayer almost simultaneously as the tagged lipid slides into the bilayer (Figure 3, CHARMM36 trajectory and Figure 4A, magenta time traces). Near the bilayer's surface, both tails are a similar distance from the bilayer ( Figure 4B, magenta region). (2) In the sn1 splayed pathway, the sn1 tail enters the bilayer first ( Figure 4A, Table  S3 reports the frequency of each pathway in our simulations. Splayed intermediates have also been observed in a previous study of lipid insertion 32 and postulated as transition states for stalk formation during membrane fusion, 47-52 a key step in vesicular transport. The existence of distinct insertion pathways might suggest that the displacement of individual tails along the bilayer's normal could serve as reaction coordinates since they encapsulate dynamically relevant information that is not contained in the COM displacement. We demonstrate below, however, that the displacements of individual tails are not the reaction coordinates for lipid exchange. The Reaction Coordinate Characterizes Hydrophobic Contacts Between the Lipid and Membrane. Based on experiments, a cavity model has been proposed to describe the transi-tion state. 16,31 According to the cavity model, solvent is evacuated above the desorbing lipid and a void forms in the membrane below. Such a focus on cavities is reminiscent of modern theories of the hydrophobic effect, which characterize hydrophobicity in terms of the statistics of solvent density fluctuations. 53 Although no true cavities are observed at transition states sampled in our simulations ( Figure S3), the importance of hydrophobicity in lipid transport is evident from our analysis of over 50 order parameters. For example, the number of water molecules solvating the tagged lipid steadily decreases during insertion ( Figure S4). The density of hydrophobic molecular fragments below the tagged lipid gradually increases while the number of hydrophilic ones decreases ( Figure S3). Defects in the polar head group region of the bilayer that expose hydrophobic membrane patches 54,55 are also observed in the transition state ensemble ( Figures S5 and S6). Based on these results, we hypothesized that the reaction coordinate for passive lipid exchange monitors the formation and breakage of hydrophobic contacts between the tagged lipid and membrane ( Figure 1B).
To rigorously test this hypothesis, we employed com- mittor analysis. [8][9][10] The committor, p B , is the probability that a trajectory initiated at a given configuration will reach state B prior to state A when initial momenta are chosen randomly from a Maxwell-Boltzmann distribution. Configurations within stable states A and B have p B = 0 and 1, respectively. Transition states are equally likely to advance to states B and A, such that p B = 0.5. [8][9][10]56 Since p B directly measures the progress of a reaction, p B is the true reaction coordinate. 9,10,57 However, p B is a complicated function of the system's microscopic configuration -a very large set of variables in the case of biomolecular systems. The complete functional form for p B is practically unobtainable and would provide little physical insight into reaction mechanisms. Instead, it is more informative to identify an order parameter, q, that is strongly correlated with p B and closely approximates the true reaction coordinate. 57,58 Henceforth we refer to such order parameters as the reaction coordinate. As one assessment of correlation between an order parameter q and the true reaction coordinate p B , we compare the probability distribution of q from transition states to distributions of q from states A and B. The typical values of q at the transition state should differ from its values in states A and B such that these probability distributions do not overlap. Otherwise, q cannot reliably distinguish transition states from stable reactant and product states and, thus, poorly recapitulates the true reaction coordinate. Additionally, we compare the probability distribution of q from transition states to distributions from pre-and post-transition states. Because pre-transition states are intermediates between state A and the transition state and post-transition states are intermediates between the transition state and state B, it is often challenging to devise an order parameter q that distinguishes transition states from them, as p B naturally does. A good approximation to the true reaction coordinate should reliably distinguish transition states from pre-and post-transition states in addition to state A and B configurations.
A more stringent test of a putative reaction coordinate examines a histogram of p B values for configurations with a particular value of q. If q is the reaction coordinate, then a histogram of p B for configurations with a given value of q will be peaked at the corresponding value of p B . This histogram test is useful to determine if the mapping from q to p B is approximately one-toone, a requirement for q to be the reaction coordinate. Bimodal histograms of p B are clear indicators that q is not the reaction coordinate. By contrast, a histogram sharply peaked at p B = 0.5 for values of q characteristic of transition states clearly indicates that q accurately describes the true reaction coordinate. 9,10 Using the first criterion that q must reliably distinguish transition states from all other configurations to be the reaction coordinate, we assessed measures of hydrophobic lipid-membrane contacts as approximations to the true reaction coordinate. Specifically, we characterize hydrophobic contacts by the minimum distance between hydrophobic carbons of the tagged lipid and hydrophobic carbons of the bilayer, min(d CC ), and the number of close hydrophobic carbon-carbon contacts between the tagged lipid and bilayer, n CC ( Figure 2B). Precise definitions are given above in the Methods sec-tion. To determine if this criterion was satisfied, we compared the probability distributions of min(d CC ) and n CC from several ensembles: equilibrium configurations representative of (1) state A and (2) state B (as defined above in Methods); and three ensembles drawn from transition paths of (3) pre-transition state configurations identified by p B = 0, (4) transition states identified by p B ≈ 0.5 (specifically, 0.45 ≤ p B ≤ 0.55 for MARTINI and 0.4 ≤ p B ≤ 0.6 for CHARMM36 configurations), and (5) post-transition state configurations identified by p B = 1. Joint distributions of min(d CC ) and n CC in these five different ensembles are shown in Figure 5. Corresponding 1D probability distributions of min(d CC ) and n CC are shown in Figures S7 and S8. The distributions from MARTINI and CHARMM36 configurations exhibit similar features. In state A, min(d CC ) effectively measures the separation between the tagged lipid and the distant bilayer. Large values of min(d CC ) are therefore typical. No close hydrophobic contacts are formed since the tagged lipid is fully solvated. As the tagged lipid progresses from state A towards the transition state, it becomes a pre-transistion state configuration. Pre-transition state configurations are closer to the bilayer than state A configurations, resulting in decreased values of min(d CC ), but hydrophobic contacts between the tagged lipid and bilayer still scarcely exist. At the transition state, the tagged lipid has made only a few initial hydrophobic contacts with the bilayer. The transition state distribution is centered at values of min(d CC ) and n CC intermediate between those of states A and B. As the tagged lipid progress from the transition state to state B, it becomes a post-transition state configuration. Post-transition state configurations, which include splayed lipid configurations, have a substantial number of hydrophobic contacts between the tagged lipid and bilayer compared to transition states but on average half as many as state B configurations. Both the post-transition state and state B distributions of min(d CC ) are sharply peaked at a value consistent with the minimum of the Lennard-Jones potential for a carbon-carbon interaction (or a bead-bead interaction in the MARTINI model). In state B, a maximal number of hydrophobic contacts exist between the tagged lipid and adjacent lipids in the bilayer. Importantly, the transition state distribution overlaps negligibly with the other four distributions. A combination of min(d CC ) and n CC distinguishes transition states from not only stable state A and B configurations but also pre-and post-transition state configurations, and, therefore, may serve as the reaction coordinate for lipid exchange.
If the true reaction coordinate is well described by a combination of min(d CC ) and n CC -in other words if the formation and breakage of hydrophobic contacts are the essential processes required for lipid exchangethe free energy surface ∆F (min(d CC ), n CC ) should exhibit a barrier for insertion and desorption. These free energy surfaces for both MARTINI and CHARMM36 are shown in Figure 6A. Corresponding 1D free energy profiles along min(d CC ) and n CC are shown in Figure  S10. A deep free energy minimum exists at values of min(d CC ) and n CC characteristic of state B due to the formation of many favorable hydrophobic contacts between the tagged lipid and membrane lipids. At larger values of min(d CC ) characteristic of state A, the free energy surface flattens; when the lipid is far away from the membrane, the free energy is no longer sensitive to min(d CC ). The surface plateaus at larger free energy values for MARTINI compared to CHARMM36. This discrepancy is consistent with the fact that compared to atomistic models, the MARTINI model overestimates the free energy difference between a lipid in solution and in the bilayer. 34 In qualitative agreement with experimental findings, [16][17][18][19] there is a free energy barrier for lipid insertion ( Figure 6A, outlined in dashes). The activation free energy for the formation of hydrophobic contacts is roughly 5 k B T . Additionally, transition paths closely follow the minimum free energy path in the space of min(d CC ) and n CC (Figure S9), indicating that the dynamics of lipid exchange are well described in terms of hydrophobic contacts.
Most importantly, a combination of min(d CC ) and n CC can locate transition states with high fidelity. We find that a linear combination suffices for this purpose (r c = α 1 min(d CC ) + α 2 n CC + α 0 with coefficients α i determined using a maximum likelihood approach 59,60 as detailed in the SI). By construction, r c = 0 at the dividing surface between states A and B where transition states are located. The dashed lines in Figure 6A outline a narrow range around r c = 0, which roughly traces the ridgeline of ∆F (min(d CC ), n CC ) between states A and B. To definitively test if r c is the reaction coordinate, we performed committor analysis of configurations drawn from a narrow range around r c = 0. Figure  6B shows that the ensemble defined by r c ≈ 0 predominantly includes transition states for both MARTINI and CHARMM36. Thus, a measure of hydrophobic contacts between the tagged lipid and bilayer is the reaction coordinate for lipid exchange. Together, min(d CC ) and n CC capture many different aspects of the lipid's hydrophobic environment, underlying the ability of r c to precisely describe the process of lipid exchange. Hydrophobicity can be quantified in other ways as well. For example, we constructed a more complicated reaction coordinate, r , that is a linear combination of 48 order parameters excluding min(d CC ) and n CC , with coefficients again determined using a maximum likelihood approach 59,60 (Table S4). As fully described in the SI, each of these 48 order parameters measure different details about the lipid's environment, including the number of water molecules solvating the tagged lipid and the size of exposed hydrophobic membrane defects near the tagged lipid. This more complex reaction coordinate identifies transition states almost as faithfully as r c does ( Figure S11A), but r (M =48) c is difficult to interpret physically. The reaction coordinates r c and r (M =48) c are highly correlated ( Figure S11B), demonstrating that detailed information about the lipid's environment is well represented by a simple combination of min(d CC ) and n CC . Not only is a linear combination of min(d CC ) and n CC the most accurate reaction coordinate out of all tested (Table S5), but it has the advantage of providing physical insight into lipid exchange: The rate limiting step for desorption is the breakage of hydrophobic contacts between the lipid and membrane.

The Lipid's Displacement Normal to the Bilayer Is Not the Reaction Coordinate.
Previous computational studies presumed that the lipid's displacement normal to the bilayer is the reaction coordinate for lipid desorption and insertion. [11][12][13][14][15] However, d lip is not the reaction coordinate, a fact we confirm by performing committor analysis. The probability distributions of d lip from transition states are compared to the distributions from state A, state B, pre-transition state, and post-transition state configurations in Figure 7. The distributions from MARTINI and CHARMM36 configurations are quite similar. In state A, where the tagged lipid is fully solvated and away from the bilayer, the typical value of d lip is large. As the lipid enters the bilayer, it progresses from being in state A to being a pre-transition state configuration, a transition state, a post-transition state, and finally in state B; correspondingly, the centers of each distribution shift to smaller values of d lip . In state B, the tagged lipid is in register with the other lipids in the membrane such that the distribution is centered near zero. For both MARTINI and CHARMM36, distributions from transition states overlap significantly with those from the other ensembles, indicating that transition states cannot be reliably identified by d lip ( Figure  1A). Furthermore, histograms of p B for configurations sampled along transition paths with values of d lip typical of transition states are bimodal and lack a peak at 0.5 ( Figure S12), indicating that d lip is not strongly correlated with p B . Therefore, d lip is not the reaction coordinate for lipid exchange.
Based on the observation of multiple pathways for lipid insertion characterized by splayed intermediates (Figure 4), combinations of the tail displacements d sn1 and d sn2 might have been considered as pathway-specific reaction coordinates. For transitions via the sliding pathway, 1 2 d sn1 + 1 2 d sn2 is a plausible reaction coordinate since both tails enter the bilayer at approximately the same time. For transitions via the sn1 splayed and sn2 splayed pathways, d sn1 and d sn2 are plausible reaction coordinates for each pathway, respectively, since the lipid is committed to fully insert after the first tail has entered the bilayer. Probability distributions of these putative reaction coordinates from MARTINI transition states are compared to the distributions from state A, state B, pre-transition state, and post-transition state configurations in Figure 8 for each pathway individually. With only 10 CHARMM36 insertion trajecto-ries, there is insufficient data to reliably examine distributions for each pathway individually (data from all CHARMM36 transitions is shown in Figure S13). With limited overlap between the transition state distribution and distributions from all other ensembles ( Figure  8), 1 2 d sn1 + 1 2 d sn2 , d sn1 , and d sn2 appear to be potential reaction coordinates that are specific to each pathway. d sn1 cannot be used to accurately identify transition states along the sn2 splayed pathway and vice versa ( Figure S14). While this indicates that transition states have values of d sn1 and d sn2 distinct from configurations in the other ensembles, it does not guarantee that 1 2 d sn1 + 1 2 d sn2 , d sn1 , and d sn2 are strongly correlated with p B , as required for them to be the reaction coordinates for each lipid exchange pathway.
Furthermore, if 1 2 d sn1 + 1 2 d sn2 is the reaction coordinate for the sliding pathway, d sn1 for the sn1 splayed pathway, and d sn2 for the sn2 splayed pathway, then they should fully describe the dynamics of lipid exchange. In that case, the free energy surface ∆F (d sn1 , d sn2 ) should have a barrier for lipid insertion and desorption to be consistent with the observed barrier crossing dynamics ( Figure 4A). This free energy surface is shown in Figure 9A for the MARTINI model. Between states A (top right corner) and B (bottom left corner), slight plateaus in the free energy profile occur at combinations of d sn1 and d sn2 characteristic of splayed configurations (top left and bottom right corners). The free energy profile has a saddle point near d sn1 ≈ d sn2 ≈ 15Å with a barrier for insertion of approximately 1 k B T (outlined in dashes), which is only slightly larger than the statistical error in the calculated free energy ( Figure S15). This small barrier would hardly impede the dynamics of insertion. Indeed, it is significantly lower than the barrier for hydrophobic contact formation ( Figure 6A), demonstrating that decreases in the tail displacments do not limit the rate of lipid insertion.
Finally, to determine if the tail displacements are strongly correlated with the true reaction coordinate, we constructed a histogram of p B values for configurations found near the saddle point during umbrella sampling simulations. The histogram of p B for configurations that are predicted to be transition states based on the free energy surface is strongly bimodal ( Figure 9B). Almost all of the tested saddle point configurations are committed to state A or B, and transitions states are seldom sampled in our simulations. Thus, the lipid's displacement normal to the bilayer, even reformulated to account for the three different insertion pathways, is not the reaction coordinate for lipid exchange.
Calculation of Lipid Exchange Rate Enables Comparison to Experiment. Having identified the reaction coordinate, which measures hydrophobic contact formation and breakage between the lipid and membrane, we utilized this knowledge to calculate the kinetic parameters of lipid exchange. First, we calculated the lipid insertion rate, k ins , using Kramers theory (Eq. 1) as detailed in the Methods section and compared it to the rate obtained directly from our unbiased trajectories. Kramers theory provides an expression for the reaction rate of a process whose dynamics are diffusive and well described by an overdamped Langevin equation. 46 During insertion, the lipid is buffeted by solvent molecules and other membrane lipids while crossing the broad free energy barrier for hydrophobic contact formation ( Figure 6A), resulting in diffusive barrier crossing dynamics ( Figure S9) and making Kramers theory appropriate. The value of k ins calculated with Kramers theory is 7.0 µs −1 for MARTINI and 21.0 µs −1 CHARMM36; k ins obtained from the mean first passage time in the spontaneous insertion simulations is 8.0 µs −1 for MARTINI and 5.7 µs −1 for CHARMM36. The insertion rates calculated based on hydrophobic contact formation differ from those obtained from dynamical simulations by a factor of 4 or less. Such good agreement demonstrates that hydrophobic contacts describe the dynamics of lipid insertion, an elementary step of lipid exchange, with quantitative accuracy. Unfortunately, the timescale of insertion is too fast to measure with straightforward experimental methods.
To compare our results to experiments, we sought to calculate instead the much slower lipid exchange rate, k ex . Although we obtained k ins directly from simulations, k ex cannot be obtained in the same manner. Complete exchange events, which involve lipid desorption and diffusion to another membrane in addition to insertion, could not be simulated due to their prohibitively long waiting times. Given the good agreement between k ins calculated directly from simulation and from Kramers theory, we instead calculated k ex from a Smoluchowski equation that contains information about the free energetics of hydrophobic lipid-membrane contacts.
To experimentally probe passive lipid exchange, two populations of vesicles, one initially composed of labeled lipids and the other initially of unlabeled ones, are combined in solution. k ex is then determined by monitoring the rate of mixing of labeled lipids between these two populations of vesicles. 17,25,26,[61][62][63] We can obtain an expression for the rate of mixing by first considering how the concentration of labeled lipids in a single vesicle changes over time. The concentration of labeled lipids in each vesicle will adopt a steady state with the concentration in the solution surrounding that vesicle. At steady state, the flux of lipids over any closed surface around a vesicle is constant. Importantly, since lipids that desorb may rapidly return to their original vesicle given the small barrier for insertion compared to desorption ( Figure 6A), the flux contains contributions from lipids both leaving and entering a vesicle. Accord-ing to the Smoluchowski equation, the flux, J, of lipids radially into a vesicle is where n B is the number of labeled lipids in a vesicle, r is the radial distance from the center of a vesicle, D is the diffusion coefficient of a lipid in solution, U (r) is the effective interaction potential between a lipid and a vesicle, and ρ(r) is the concentration of labeled lipids at a distance r. As for the calculations using Kramers theory, we use min(d CC ) to describe interactions between a lipid and a vesicle in 3D space. U (r) is then the free energy profile as a function of min(d CC ) ( Figure S10). We set U (r) = 0 far away from a vesicle and shift the peak of the free energy barrier to r = 50 nm, a typical radius of large unilamellar vesicles (LUVs) used in experimental studies. 17,25,26,[61][62][63] Assuming that the concentration profile of labeled lipids reaches steady state much faster than n B is varying, the time dependence of the bulk concentration of labeled lipids,ρ, can be regarded as constant while calculating the steady-state profile ρ(r)/ρ. Assuming as well that equilibration within state B is fast, the concentration of labeled lipids within a vesicle obeys equilibrium statistics, where i indicates whether the vesicle was initially composed of labeled (i = 1) or unlabeled (i = 2) lipids, r B is a particular location within a vesicle, q B is the partition function for labeled lipids in a vesicle, and the integral is performed over the range of r designated as state B. With these two assumptions, Eq. 2 integrates to where r c is the distance beyond which U (r) is zero. Finally, the rate of mixing is obtained by considering how the difference between the number of labeled lipids in type 1 vesicles and in type 2 changes over time: Eq. 5 correctly reflects the fact that lipid exchange is a first order process. 16  DLPC, our calculated k ex is in agreement with experimental values. However, in the case of CHARMM36 DMPC, our calculated k ex differs from experimental values by six orders of magnitude. Some discrepancy between experiment and simulation/theory is to be expected because small errors in the free energy profile are magnified exponentially in k ex . Yet, a six-orderof-magnitude difference between theoretical and experimental rates corresponds to a significant underestimation of the free energy barrier of approximately 12 k B T for CHARMM36. This large discrepancy could be due to any number of the following reasons: (1) Sampling errors, which are within about 1 k B T throughout the free energy surface ( Figure S9), may influence our rate calculations.
(2) Inaccuracies in the force field may cause errors in the calculated free energy surface. Based on reported differences between permeabilities and partition coefficients for water, alkanes, and other small molecules obtained from simulations with CHARMM36 and from experiments, 64,65 we estimate that errors due to the force field are approximately 1 − 2 k B T . We also repeated our calculations using a different force field (Stockholm lipids, also known as Slipids 66,67 ) as detailed in the SI. The height of the free energy barrier is approximately 2 k B T higher for Slipids compared to CHARMM36 ( Figure S16), resulting in an estimate of k ex that still differs from experiment by about five orders of magnitude.
(3) The ionic strength of the solution may influence the free energy surface. Our simulations are performed with neat water, but experimental systems are conducted in buffer, which contains significant amounts of monovalent salts, such as NaCl or KCl. 17,61-63 The addition of NaCl or KCl could salt out the free lipids in solution, increasing the difference in free energy between state A and B. Additionally, monovalent cations bind to the carbonyl region of phosphatidylcholine membranes, causing the thickness of the bilayer and order of the tails to increase. [68][69][70] These structural changes are expected to increase the free energy barrier for hydrophobic contact formation.
(4) The assumptions we made to obtain Eq. 4 may not be consistent with what occurs in experi-mental systems, causing discrepancies between our calculated k ex and experimental measurements. In writing down a Smoluchowski equation (Eq. 2), we assumed that lipids exchange between stable, spherical vesicles of uniform size. This may be justified since vesicles composed of biological phospholipids are kinetically-trapped, metastable states, which by definition do not quickly relax into equilibrium lamellar structures. [71][72][73][74] But, we lack the requisite knowledge about the processes and associated timescales by which vesicles relax into equilibrium structures to confirm the validity of this assumption. Other relaxation processes, such as uncatalyzed vesicle fusion, which has an experimentally measured rate similar to lipid exchange, 75 could occur simultaneously. We note that lipid flip-flop, the process by which a lipid moves between leaflets of the same bilayer, 1,5 also has a rate measured under some experimental conditions that is similar to the exchange rate. 61,76 Unlike vesicle fusion, lipid flip-flop is accounted for in experimental determination of lipid exchange rates.
(5) Similar issues of random error and unjustified kinetic assumptions could plague the inference of exchange rates from laboratory measurements.

CONCLUSIONS
The breakage of hydrophobic contacts limits the rate of passive lipid transport. To reach this conclusion, we investigated the elementary steps of lipid exchange, lipid insertion into and desorption from a membrane, using molecular dynamics simulations of the coarse-grained MARTINI and all-atom CHARMM36 lipid models. Results from MARTINI and CHARMM36 provide consistent pictures; even with a coarse description of lipids and water, the MARTINI model captures the essential features of lipid exchange exhibited in an allatom model. We discovered that the reaction coordinate for passive lipid exchange measures the formation and breakage of hydrophobic lipid-membrane contacts, which gives rise to a free energy barrier for both lipid desorption and insertion. Thus, knowledge of the reaction coordinate resolves previous qualitative discrepancies between simulations, which predicted that there is no barrier for lipid insertion, and experiments, which indicated that there is a barrier. This barrier likely plays an important biological role: A barrier for lipid insertion ensures that membrane compositions, which are spateotemporally regulated to maintain cell homeostasis, 1,2 are not easily disrupted. We suspect that the breakage and formation of hydrophobic contacts may also give rise to a free energy barrier for transporting synthetic amphiphiles. 77 Additionally, knowledge of the reaction coordinate allowed us to formulate a Smoluchowski equation to model lipid exchange between vesicles, which occurs over time and length scales inaccessible in MD simulations, and calculate the rate of lipid exchange. Differences between our calculated lipid exchange rate and experimental measurements indicate that considerable quantitative discrepancies between simulation and experiment still exist. Future studies will be performed to better assess the sources of these discrepancies.
Finally, this knowledge provides a foundation to understand how catalysts of lipid exchange work at a molecular level. Lipid transfer proteins may efficiently extract lipids from membranes by lowering the activation free energy barrier for hydrophobic contact breakage. Interestingly, catalysts of vesicle fusion, the key step in vesicular lipid transport, may function in a similar way. For example, carbon nanotubes aid vesicle fusion by facilitating the formation of hydrophobic contacts between two vesicles, 78 and viral fusion peptides are thought to catalyze fusion by promoting hydrophobic lipid tail protrusions and contacts. 47 Thus, common physical properties of lipids, specifically their hydrophobicity, may be exploited in vivo to precisely control both non-vesicular and vesicular lipid transport.
Acknowledgement J.R.R. acknowledges the support of the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE 1752814. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562; specifically, it utilized XSEDE Comet at the San Diego Supercomputing Center through allocation CHE180038. We thank Prof. Baron Peters for providing us with maximum likelihood approach code. We thank Dr. Georg Menzl and Layne Frechette for valuable discussions and helpful comments.

Supporting Information Available
Complete descriptions of MD and umbrella sampling simulation protocols; descriptions of all order parameters used to analyze transition paths; description and results of maximum likelihood approach to identify a reaction coordinate; demonstration of equilibration of initial lipid bilayers; probability distributions used to determine stable state definitions; probability distributions of additional order parameters; statistical errors in computed free energy surfaces; minimum free energy paths along min(d CC ) and n CC ; 1D free energy profiles of min(d CC ) and n CC ; results from simulations with the Slipids force field; frequency of each insertion pathway; movies of lipid insertion trajectories.

Graphical TOC Entry
Hydrophobic Contact Free Energy