Membrane Permeability of Fatty Acyl Compounds ... - ACS Publications

Subscriber access provided by Gothenburg University Library

Article

Membrane Permeability of Fatty Acyl Compounds Studied via Molecular Simulation Josh V. Vermaas, Gregg T Beckham, and Michael F. Crowley J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b08233 • Publication Date (Web): 17 Oct 2017 Downloaded from http://pubs.acs.org on October 18, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Membrane Permeability of Fatty Acyl Compounds Studied via Molecular Simulation Josh V. Vermaas,† Gregg T. Beckham,∗,‡ and Michael F. Crowley∗,† †Biosciences Center, National Renewable Energy Laboratory, Golden, CO 80401 ‡National Bioenergy Center, National Renewable Energy Laboratory, Golden, CO 80401 E-mail: [email protected]; [email protected] Abstract Interest in fatty acid-derived products as fuel and chemical precursors has grown substantially. Microbes can be genetically engineered to produce fatty acid-derived products that are able to cross host membranes and can be extracted into an applied organic overlay. This process is thought to be passive, with a rate dependent on the chemistry of the crossing compound. The relationship between the chemical composition and the energetics and kinetics of product accumulation within the overlay is not well understood. Through biased and unbiased molecular simulation, we compute the membrane permeability coefficients from production to extraction for different fatty acyl products, including fatty acids, fatty alcohols, fatty aldehydes, alkanes, and alkenes. These simulations identify specific interactions that accelerate the transit of aldehydes across the membrane bilayer relative to other oxidized products, specifically the lack of hydrogen bonds to the surrounding membrane environment. However, since extraction from the outer membrane leaflet into the organic phase is found to be rate limiting for the entire process, we find that fatty alcohols and fatty aldehydes would both manifest similar fluxes into a dodecane overlay under equivalent conditions, outpacing the accumulation of acids or alkanes into the organic phase. Since aldehydes

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

are known to be highly reactive as well as toxic in high quantities, the findings suggest that indeed fatty alcohols are the optimal long-tail fatty acyl product for extraction.

2


Page 2 of 46

Page 3 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Over the past two decades, significant effort has been invested in engineering microbes to create oleochemicals. In oleaginous yeast, genetic engineering can guide metabolism towards the overproduction of fatty acids, 1,2 which can be used as a feedstock for a variety of fuels and chemicals, including diesel replacements, detergents, and bio-derived plastics. 1,3–5 However, once the product has been synthesized, it must be extracted from the yeast in some way. One option is for cells to be lysed after product accumulation. 6–8 This approach has two major cost drivers; the direct cost of the lysing reagents, as well as the indirect cost of regrowing the majority of the microbial biomass lost in the process. Alternatively, it has been found that organic overlays in the growth media can extract fatty acids or especially fatty alcohols in a non-destructive manner, significantly reducing cost for the end products by eliminating the need to regrow host cells. 2,9–11 What has not been shown, however, is whether the products being extracted from the cells are ideal in terms of maximizing the rate of product extraction. Enzymes exist for generating alkanes, 12 alkenes, 13,14 and aldehydes 15 in addition to alcohols from fatty acids, 10 and it has not been conclusively determined which of these products have the highest membrane transit (kt ) or dodecane adsorption (kd ) rates (Fig. 1). Indeed, depending on bilayer transit kinetics of the compound produced, dodecane overlays can either quickly accumulate the desired compound, increasing yield significantly 16,17 or only modestly increase productivity. 18 If fatty acid products accumulate within the yeast membranes, rather than being extrated, the high concentration of product would effectively inhibit further fatty acyl product biosynthesis through competitive inhibition, in turn reducing the yield. To maximize productivity, the product should be chosen such that it quickly and efficiently transfers across the bilayer and into the dodecane overlay. Some preliminary hints as to the properties of the optimal product come from the pharmaceutical industry, where the challenge is to deliver drugs as efficiently as possible to the cell interior, effectively reversing the process pursued here. For a compound to be internalized, 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 4 of 46

Page 5 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


it must (1) leave its carrier, which can be a complex organic lattice such as a liposphere 19 or organic nanosheet, 20,21 (2) enter the outer leaflet of the cell membrane, and (3) flip across the membrane bilayer, which is exactly the reverse of the process under study here (Fig. 1). It has been recognized that typically these compounds should be sufficiently lipophilic to permeate across the membrane, but also soluble enough to be distributable by the bloodstream. 22–24 The solubility and lipophilicity therefore need to be balanced for best results, 25 a pattern repeated here for fatty acyl product extraction. Evaluating the specific interactions that influence the kinetics of membrane transit or adsorption can be tricky through an experimental lens, as neither ensemble or single-molecule experiments alone can fully explore the interactions between the molecule and its environment responsible for the observed kinetics. By contrast, computational investigation of the process kinetics and their molecular origins is much simpler, and is why there has been so much success in investigating permeation of drug-like compounds through simulation. 26–28 A controlled environment is setup, and the inherent atomic resolution of atomistic molecular dynamics simulations permits specific causal interactions to be directly identified. 29 Even for slow processes, combining free energy profiles and local diffusion constants from biased simulation is enough to determine membrane permeability, 30,31 which in turn predicts potential fluxes given a set concentration of product. 32 Through a combination of equilibrium and biased atomistic molecular dynamics (MD) simulations, we determine the relative kinetic parameters for excreting the fatty acid products into dodecane or water. From our findings, we observe that extraction from the cell membrane to a dodecane overlay is rate-limiting for most of the products tested, and as a result both alcohols and aldehydes have similar computed membrane permeability coefficients. However, since aldehydes lack a hydrogen bond donor, and thereby preferentially segregate below the membrane headgroups, aldehydes transfer between membrane leaflets very quickly during equilibrium simulation. Fatty acids or alcohols, by contrast, transfer between leaflets significantly more slowly. If the acyl tail length of the product were shorter such that transfer rather than

5




Page 6 of 46

Page 7 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cell (Nin ), within the inner leaflet (Nil ), within the outer leaflet (Nol ), and outside the cell (Nout ) over time, and relate this to the flux (J, here expressed in molecules per unit time) between adjacent partitions within the cell. N˙ in = Jproduction − Jin→il N˙ il = Jin→il − Jil→ol

(1)

N˙ ol = Jil→ol − Jol→out N˙ out = Jol→out − Jaccumulation Assuming the cell is at a steady state, no further products can accumulate within any partition (N˙ = 0), and we see that equation 1 simplifies to equalizing the fluxes between partitions (Fig. 2). Flux in general is a product of the difference in concentration between adjacent partitions (∆C), the surface area of the interface between the system partitions (A), and the permeability coefficient between those partitions (P):

J = P × A × ∆C

(2)

Since the geometry of a cell does not change depending on the fatty acyl product, and the partition coefficient is determined by properties of the membrane and the species itself, the concentration difference between individual partitions responds to keep up with the production rate. As a result, the final concentration within the dodecane overlay at steady state depends on the permeability coefficients for individual steps. These unknown permeabilities are what will be computed through molecular simulation. As described in recent work by Lee et al. 30 , bilayer crossing kinetic parameters can be efficiently computed given only the free energy profile and the position-dependent local diffusivity using the inhomogeneous solubility-diffusion model based on the work of Diamond and Katz 33 and first adapted for simulation use by Marrink and Berendsen 31 , which states

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 46

that the permeation resistance, R, is given by:

R=

Z

ξu ξl

exp (∆G (ξ) β) dξ D (ξ)

(3)

Where the free energy change (∆G) and diffusivity (D) are known along the reaction coordinate ξ between two points, which in the original development of the theory is simply the position along the membrane normal axis. This permeation resistance is inversely related to the permeability P (P = R−1 ), which is equivalent to the permeability coefficient in Eq. 2 for each component of the egress process. The apparent permeability for the entire process can also be modeled by summing together the resistivities from each step of the excretion −1 process (Rtotal = Rin→il + Ril→ol + Rol→out = Papparent ), similar to Ohm’s Law as applied

to a series circuit. Expressed in permeability form (analogous to conductance in electrical circuits), this is: Papparent =

1 1 Pin→il

+

1 Pil→ol

+

1 Pol→out

(4)

Methods To explore the phenomenon of bilayer crossing and extraction of fatty acyl products, there are two classes of simulation carried out. Unbiased equilibrium MD simulations were sufficient to observe membrane leaflet transfer of some products, and to arrive at relaxed membrane states suitable as starting points for the biased simulations required for free energy profile calculations. Since some processes of interest, in particular transfer of the product to the simulated dodecane overlay, occur on timescales well beyond what is achievable with current simulation techniques, the equilibrium simulations were supplemented with replica exchange umbrella sampling simulations 34 (REUS, also called Hamiltonian replica exchange 35 or bias exchange umbrella sampling 36 ) to sample along the transitions represented by bilayer flipping and bilayer extraction (Fig. 1). The aggregate equilibrium and non-equilibrium simulation time exceeds 37.5 µs. The details of setup, simulation, and analysis are presented in the 8


Page 9 of 46




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

having a 30:14:12:9:4:1 molecular composition of phosphatidyl choline (PC), phosphatidyl ethanolamine (PE), phosphatidyl inositol (PI), ergosterol, phosphatidyl serine (PS), and phosphatidic acid (PA) headgroups and sterols. This ratio was determined by examining different literature sources of experimentally determined headgroup compositions of the prototypical yeast Saccharomyces cerevisiae, 38–40 and specifically selecting for strains where composition is known for all major lipid classes as well as ergosterol. There are many strains of Saccharomyces cerevisiae to choose from, each with their own distinct lipid composition, however the W303-1A strain composition was selected as a compromise between different compositions observed. 38 While the membrane composition is derived from non-oleaginous yeast species, available mass spectrometry data for a variety of yeast species, including oleaginous yeast, suggest that the headgroup composition used in the model is reasonable in its proportions, but critically lack information as to sterol concentrations present in the membrane to complete the model. 41 Mass spectrometry data informed the choice of lipid tails. Most lipid tails seen from yeast cells have 18 carbons and one unsaturation (C18:1), which are approximately as abundant as saturated 16 carbon tails (C16:0) and 18 carbon tails with no or two unsaturation sites (C18:0, C18:2). 42,43 55% of our lipids are thus of the palmitoyl oleyl class (C16:0-C18:1), 22% of our lipids are of the stearyl-linoleyl class (C18:0-C18:2), and 23% of our lipids are of the dioleyl class (C18:1-C18:1). This ratio has the correct composition in terms of tail components, but likely does not reflect the true tail component ratios, which have many more possible combinations. The membrane was assembled with CHARMM-GUI, 44–46 and each leaflet had 9 additional fatty acyl products inserted into it, for a total of 18 product molecules per system (green molecules in Fig. 3). This corresponds to fatty acyl products making up 5% of the membrane by weight, an estimate that balances likely concentrations during fatty acyl overproduction with product toxicity observed at higher concentrations of product. 5,47,48 In addition, the additional copies permit sampling of product-membrane interactions from different micro-environments within the membrane arising from the initial distribution of lipids. The fatty acyl products considered were palmi-

10


Page 10 of 46

Page 11 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tate (fatty acid), cetyl alcohol (fatty alcohol), cetyl aldehyde (fatty aldehyde), pentadecane (fatty alkane), and 1-pentadecene (fatty alkene). Dodecane overlays were constructed prior to solvation and ionization by the Solvate and Autoionize plugins of VMD. 49 In the case of a 10 % dodecane overlay, 71 dodecane molecules were placed randomly in solution through the insert-molecules tool from the GROMACS simulation suite, 50,51 which quickly aggregated into a small sphere after only a few nanoseconds of simulation (Fig. 3B). For the organic phase, 2000 molecules were placed in a box as above, but then equilibrated to the correct density. The equilibrated box was then placed near the membrane, and elements were removed if they came within 5 ˚ A of membrane components. Ionization brought the total aqueous salt concentration to 150 mM NaCl, which has a surplus of up to 48 Na+ ions due to the net charge of the anionic membrane and the fatty acyl product. Even in the case of the organic phase, which has the least water present, there are at least 12 water molecules per leaflet component to solvate the headgroup layer. All systems have an initial total size of approximately 65 ˚ A×65 ˚ A×100 ˚ A.

Simulation Details There are common simulation conditions for both the equilibrium and REUS simulations, which are listed here for completeness. All simulations were performed with NAMD 2.11, 52 using the CHARMM36 force field for lipids 53 and sterols. 54 Both dodecane and the fatty acyl products used lipid parameters to describe their motion, using the CHARMM general force field (CGenFF) 55 to supplement the lipid force field for chemistries unique to the products (e.g. the aldehyde). The TIP3P water model 56 was used in conjunction with a 12 ˚ A nonbonded cutoff, consistent with the parameterization methodology for CHARMM force fields. 55 Long-range electrostatics was treated using particle mesh Ewald 57,58 with a 1.2 ˚ A grid. Anisotropic pressure coupling was maintained via the Langevin piston method 59,60 to a pressure of 1 atm, with periodic cell growth along the membrane normal axis decoupled from growth along the membrane parallel axes. A Langevin thermostat using γ=1 ps−1 maintained 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the system temperature 61 at either 300 K (equilibration) or 310 K (biased simulation). Each simulation timestep was 2 fs, enabled by using the SETTLE algorithm 62 to fix bond lengths to hydrogen. Equilibrium Simulation Unbiased simulation was carried out to equilibrate the membranes constructed above in preparation for biased bilayer extraction or crossing simulations. These unbiased simulations were each at least 100 ns long, with separate simulations for each combination of fatty acyl product and extraction environment. At this timescale, the membrane arrives at its preferred geometry, as demonstrated by the quick relaxation of the periodic box length parallel to the membrane plane (Fig. S1), although this is well short of the timescales required to see significant mixing or clustering. 63,64 It was noted that for aldehydes, alkanes, and alkenes, membrane crossing happened spontaneously, as did the fatty alcohol at a very slow rate. These four equilibrium product simulations were extended to 500 ns to better sample these spontaneous membrane crossing events. Another equilibrium simulation was carried out with the alcohol atomic charges altered to match those of the aldehyde to probe the effect of hydrogen bonds on controlling spontaneous membrane crossing. Bilayer Crossing Simulations For the cases where the membrane crossing was not spontaneous (the acid), or happened very slowly (the alcohol), the transition between leaflets was sampled through biased simulation. Since the protonation state of carboxylic acids has been previously shown to have a strong effect on membrane crossing, 65,66 additional simulations were carried out for palmitic acid in addition to palmitate. The reaction coordinate to track the progress of the crossing was chosen to be the position of the C1 atom, the first carbon in the acyl chain, as implemented in the collective variables (colvars) module of NAMD. 52,67 This consistent choice of reaction coordinate enables direct comparison between the free energy profiles produced here from

12


Page 12 of 46

Page 13 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


biased simulation and profiles inferred from population distributions observed in equilibrium simulation. The state at 100 ns from the equilibrium simulations was used to seed the starting point for the steered molecular dynamics (SMD) simulations, which in turn provided the starting configurations for the REUS simulations. In 18 separate SMD simulations, a single copy of the fatty acyl product in the simulation was pulled across the bilayer over 10 ns using a force ˚−2 . Due to the thickness of biological membranes of approximately constant of 2 kcal mol−1 A 40 ˚ A, 68,69 the tested range of C1 positions ranged from -25 ˚ A to 25 ˚ A along the membrane normal, with the zero-point located at the membrane midplane. This 50 ˚ A range for the reaction coordinate was subdivided into 72 windows, as window spacings between 0.5 and 1˚ A have worked well for previous studies of bilayer translocation. 30,70–72 Since part of the advantage of REUS simulations comes from their diverse sampling of different configurations, we enhance this by drawing 4 initial states for the REUS simulations from each SMD trajectory. Because each SMD trajectory is pulling a different copy of the fatty acyl product across the membrane which starts in a different position within the membrane, our final profile will be minimally biased by the local microenvironment around a single product molecule. The 4 states contributed by each SMD trajectory are chosen such that they are the closest states to the center of 4 randomly assigned windows. The REUS simulations themselves were conducted with an umbrella force constant of ˚−2 relative to the center of the umbrella for a given replica. Exchanges be4 kcal mol−1 A tween adjacent replicas were attempted every picosecond between alternating replica pairs, yielding a 2 ps overall exchange rate, a rate consistent with reported best practices suggesting frequent exchange attempts. 73 The combination of these parameters yielded a 20% exchange acceptance ratio, which has been argued to be optimal for replica exchange simulations. 74 Each replica was simulated for 20 ns, where independent free energy profiles for the first 10 ns and last 10 ns were typically very close to one another, indicating convergence (Fig. S2).

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 46

Bilayer Extraction Simulations In no case was extraction of the fatty acyl product into either of the dodecane phases or aqueous solution spontaneous on the timescale simulated. Thus, all 5 fatty acyl products required biased simulation. However, unlike the bilayer crossing presented above, pulling on the C1 carbon alone would result in a highly unphysical reaction coordinate, as the acyl tail within membrane lipids is known to come to the membrane surface well before the desorption process is complete. 75,76 Using a contact-based reaction coordinate in this case is most appropriate, as the system itself can choose how to reduce or increase the number of contacts between system components as the desorption progresses. Specifically, this is implemented as a “coordnum” (coordination number) collective variable within the colvars module 67 of NAMD 52 with the following form: 4 XX 1 − |xi − xj | /8˚ A C (g1 , g2 ) = 10 ˚ 1 − |x − x | /8 A i j i∈g1 j∈g2

(5)

This formalism calculates the number of weighted contacts between two groups (g1 and g2 ). Based on the groups, the weighted contact number goes to zero if the closest approach between the two groups is around 10 ˚ A (Fig. S3). When the closest approach distance is larger, there will be very few pairs of atoms that contribute significantly to the overall sum. If the point of closest approach for an atom pair is further away, the small number of atom pairs that will contribute to the overall sum do not contribute very much, and the total contact number remains small. It should be noted that the exponents within Eq. 5 are in principle arbitrary, however this particular combination both captures a close contact (distance less than 2.5 ˚ A) and falls off to small values at an appropriate distance from the membrane. Depending on the simulation conditions, the calculated contacts and the identities of g1 and g2 could change. For instance, for extraction into solution (Fig. 3A), one group consisted of the last 15 acyl carbons from each product (C2-C16 for the acid, aldehyde, and alcohol,

14


Page 15 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


C1-C15 for the alkane and alkene), while the other was the even-numbered acyl tail carbons from the membrane lipids to reduce the performance impact of computing the collective variable. Under these selected atoms and this particular formulation of the collective variable, the observed contact number at equilibrium is approximately 600, representing a fully membrane embedded state. With this formulation, a fully extracted product molecule floating in solution would have a contact number of 0. To capture the full range of possible contact numbers along the extraction reaction, our reaction coordinate range was [0,625], and was split into 54 equally spaced replicas. The same SMD preparation was used to setup each replica with a 0.01 kcal mol−1 contact−2 force constant applied to extract the fatty acyl product. Due to alternative atomic configurations that result in zero product-membrane contacts but a create a void within the membrane (Fig. S4), harmonic restraints along the membrane normal axis maintained the bilayer thickness and disfavored void formation. These restraints ˚−2 , and were applied to only the bilayer phosphorus had a force constant of 0.1 kcal mol−1 A atoms. Each of the 18 product molecule extraction simulations contributed 3 states to the 54 total replicas, which were simulated for 20 ns using the same 0.01 kcal mol−1 contact−2 force constant to harmonically restrain each replica to its respective center. Exchanges between replicas were attempted again every picosecond between alternating replica pairs, for a 2 ps overall exchange attempt rate. At this force constant, 30% of exchanges were accepted. PMF convergence and replica tracking in each simulation condition are reported in the Supporting Information (Figs. S5 and S6). When dodecane is included, an additional contact number is also computed, where one group is again the last 15 acyl carbons of one product molecule, but the second group is a subset of carbon atoms of the dodecane phase, specifically C1, C2, C4, C6, C7, C9, C11, and C12. Taking this symmetric subset of atoms within dodecane reduces the performance impact of computing the reaction coordinate, while not introducing a bias in the reaction coordinate. The total reaction coordinate is then represented by a combination of the two contact numbers, Ctotal = Cmembrane − Cdodecane . In this case, a positive value for Ctotal

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 46

would reflect a product that is making its way out of the membrane, whereas a negative value reflects a product that is extracted into the dodecane phase. In all but one case, the tested range for this total reaction coordinate is [-625,625], which was covered by 108 replicas during REUS simulation, and spans the extracted and membrane embedded states. For the fatty aldehyde extraction into the organic phase, an extra restraint was applied to ensure that the aldehyde group was not in the middle of the membrane. The restraint is described fully in the Supporting Information, along with a comparison of results relative to the original unrestrained set of simulations (Figs. S7 and S8). The SMD setup for each window was similar to the procedure for solution extraction simulations, using the same force constants and constraints, except that each SMD trajectory contributed 6 states for each pulled product molecule to cover the larger replica set. Analysis Analysis of the simulation trajectories was primarily carried out using a combination of builtin and purpose built VMD 49 scripts and features, leveraging the numpy and scipy libraries 77 to handle and manipulate trajectory data. Individual dataplots were constructed with the matplotlib library. 78 The inhomogeneous solubility diffusion relationship (Eq. 3) was first established for a reaction coordinate along the membrane normal axis, which has a natural unit of length, the width of the membrane. However, Eq. 3 also holds for other reaction coordinates, such as contact number, which does not have a natural unit of length, and would thus cause the general flux equation (Eq. 2) to have units with little experimental meaning. Thus, we instead report a modified quantity P ∗ when discussing membrane extraction reaction coordinates, where: P∗ =

∆d P ∆ξ

(6)

Here, ∆ξ is the extraction process described in reaction coordinate space, and ∆d is the mean change in C1 position between the two end states. Replacing P with P ∗ in Eq. 2 yields 16


Page 17 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fluxes that can be directly compared for these cases, and importantly also added together to determine an apparent permeability coefficient for the complete process. The free energy profiles themselves are computed using Bayesian reweighting combined with Gibbs sampling, 36 based on the work of Habeck 79 and Bartels 80 , and recently synthesized by Ferguson 81 , which reformulates the traditional weighted histogram analysis method (WHAM) 82 in probabilistic terms. The local diffusivity is estimated using the formalism put forth by Hummer 83 , directly from the variance and the reaction coordinate autocorrelation function at different restraint positions (Fig. S9). Details of the calculation are provided as Supporting Information. Due to the potential for rapid exchange between windows, the integrated autocorrelation time required to compute the diffusivity is estimated by extrapolating the approximately exponential decay of the reaction coordinate autocorrelation observed on short timescales to long time.

Results and Discussion Given the discussion of the underlying theory, it is important to identify the expected rate limiting steps. For typical fatty acyl products (fatty acids and alcohols), the presence of hydrophilic moieties is expected to slow the rate of bilayer transfer, as membrane crossing free energy profiles for amino acid analogs show large barriers to the translocation of polar species. 70 By contrast, for more hydrophobic components, it is expected that the extraction into solution or an organic phase would be rate limiting to the overall process, as crossing itself is spontaneous on comparatively short lengthscales for compounds like cholestrol. 28 To determine the permeability for both of these steps, we present first the results from bilayer crossing simulations, including their equilibrium distribution and the effect of hydrogen bonds in addition to the permeability computed from biased simulation. Thereafter, we analyze and present the results pertaining to bilayer extraction.

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Bilayer Crossing of Fatty Acyl Products From the equilibration simulations, there was one striking and unexpected feature. Whereas the fatty acid does not exchange between the two membrane leaflets present in the simulation, and fatty alcohols exchange only three times in our simulations (Figs. S10 and S11), the fatty aldehyde spontaneously and frequently exchanges between membrane leaflets on the relatively short sub-microsecond timescales studied here (Fig. S12). This aldehyde behavior is closer to what is observed for the alkane and alkene, which exchange rapidly between leaflets (Figs. S13 and S14). The overall distribution averaged over the full trajectory length is reported in Fig. 4, which highlights that the probability distribution for the aldehyde between the two peaks in either leaflet is connected, unlike that of the fatty acid or the fatty alcohol. However, unlike the alkane or alkene, the aldehyde products do not lie in the middle of the bilayer for significant periods of time, and instead quickly transit from one side of the bilayer to the other, with a minimal amount of time spent in the center of the bilayer (Fig. 4). Indeed, the alignment of the aldehyde product relative to the membrane normal is quite similar to an acid or alcohol (Figs. 4 and S15), and may even be “upside down”, with its oxygen within the membrane interior. The physical basis for this change in transferring between leaflets lies in the hydrogen bonding of the different products with the membrane environment (Fig. 5). The carboxylic acid carries a charge, and therefore readily interacts with water and proton donors in the membrane headgroup region. This is in part also why the carboxylic acid is distributed furthest from the membrane center (Fig. 4). Alcohols can act as hydrogen bonding donors, and thereby can interact with both headgroup phosphate groups as well as the exposed carbonyl groups of the lipid tails, in addition to acting as both a donor and acceptor to water molecules near the membrane interface. An aldehyde, in contrast to the acid, is charge neutral at physiological pH, and can only act as an acceptor to hydrogen bonds, unlike the alcohol, which can also donate a proton. There are few hydrogen bond donors at the membrane interface, and many potential acceptors, such as the acyl tail carbonyls or the 18


Page 18 of 46

Page 19 of 46





Page 20 of 46

Page 21 of 46




1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

where required (Fig. 6). The free energy profile follows trends previously seen in other systems, where a deprotonated carboxylate, in this case palmitate, encounters a high barrier to bilayer crossing due to solvating waters crowding around the naked charge as it crosses the membrane, creating a water defect within the membrane as the water molecules follow the charge across a membrane. 84,85 This is ameliorated considerably by forming a carboxylic acid at the membrane interface, which lowers the observed barrier to membrane transit by over 12 kcal mol−1 , to a significantly smaller 5 kcal mol−1 . Thus, even though the palmitic acid is a small fraction of the palmitate at neutral pH, the acids will still make up the majority of the population of crossing fatty acids. 65,66 The barrier to bilayer crossing for the alcohol is approximately double the barrier experienced for the aldehyde, indicative of precisely how unfavorable it is to break the previously described hydrogen-bond-like interactions between a fatty alcohol and the headgroups and waters at the membrane interface. The free energy profiles in Fig. 6 also put the non-symmetric probability distributions from Fig. 4 into their appropriate context, indicating that even in the worst case of probability asymmetry, that of the aldehyde, the free energy profile asymmetry does not exceed 0.5 kcal mol−1 . This bounds the systematic sampling error in this case, which represents the deviation from the long time behavior of the product distribution between the two leaflets being symmetric, given the symmetric membrane composition of the two leaflets. When transitions are more frequent, as is the case for the alkane or alkene, this systematic error is reduced. The error bars reported for the biased transitions from Fig. 6 are purely statistical, and do not account for errors due to limited sampling. Again, these systematic errors can be estimated from the asymmetry of the free energy profiles, and is approximately 1 kcal mol−1 in the case of the worst performer, the palmitate. Based on comparison of the first and second half of each REUS trajectory set (Fig. S2), significantly longer simulations would be required to completely eliminate this systematic error. The shape of the local diffusion profile determined from the simulations (Fig. 6) reflects local bilayer structure. In regions where the membrane is more disordered, such as near the

22


Page 22 of 46

Page 23 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


membrane center, 46,86 the local diffusion constant increases. This is not a universal feature of local diffusion constants within a bilayer, 30,87 although it has been measured before for some small compounds. 88,89 The rationale offered typically is that the center of the bilayer is the region of lowest density, and creates an environment where small molecules face little resistance from fewer atomic collisions. Despite our fatty acyl products being relatively large, the measured reaction coordinate only considers the C1 atom, and so effectively acts like a much smaller molecule. When |z| > 15 ˚ A, indicating the C1 atom is leaving the membrane core, the diffusion constant also increases as C1 is gradually exposed to the solvent environment, where small water molecules buffet the atom rather than being shielded within the membrane core. Table 1: Permeability coefficients (P) for bilayer crossing and bilayer extraction depending on the fatty acyl product. Each permeability is computed using the relationships associated with Eq. 3, with the upper and lower bounds of the integral taken to be the locations of the local free energy minima that delineate the process. The first column shows the permeability coefficient for leaflet-to-leaflet transfer, in the customary units of cm s−1 . The other columns are for the permeability coefficient equivalent for product extraction from the bilayer, which has been converted into customary units (cm s−1 ) from the natural units of the reaction coordinate used (contact number ns−1 ). The different reported permeability coefficients for extraction correspond to the conditions shown in Fig. 3; extraction to the aqueous phase (aq), a dodecane blob that is 10% of the total solvent volume (blob), and a dodecane phase above the membrane (phase). The equivalent extraction permeability coefficients in natural units are given in Table S1. Product Palmitate Palmitic Acid Cetyl Alcohol Cetyl Aldehyde Pentadecane 1-Pentadecene

Bilayer Crossing log10 P cm s -11.0 -1.6 -1.2 0.2 1.4 1.1

∗ log10 Paq

-1.7 – -3.2 -2.3 -5.3 -4.4

cm s

Extraction ∗ cm log10 Pblob s -1.6 – -2.4 -1.4 -3.9 -3.1

∗ log10 Pphase -1.6 – -1.6 -0.9 -2.3 -2.4

cm s

Together, the free energy profile and the local diffusion constant can be combined with the inhomogeneous solubility diffusion model (Eq. 3) to arrive at estimates for the permeability coefficient (Table 1). From the perspective of the bilayer crossing rates, the combination of 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the free energy and the local diffusion coefficients reinforces earlier findings of spontaneous crossing; species that were observed to cross spontaneously have the largest permeability coefficient, while those with the smallest permeability coefficients exhibit significant barriers to crossing. The most noteworthy aspect is the quantification of the bilayer transit that the aldehyde provides over other non-reduced species. Fatty aldehydes transit approximately 25 times faster than do fatty alcohols, about 60 times faster than even protonated fatty acids. By comparing the number of crossing events in Figs. S11 and S12, we see that this estimate is consistent with what we see in equilibrium simulations, where there are 3 alcohol crossing events in 500 ns, and between 60 and 70 (depending on how a crossing is defined) in the case of the aldehyde. There is the potential for accelerating the membrane crossing step by emphasizing aldehyde production.

Bilayer Extraction of Fatty Acyl Products The end result of our study into the extraction kinetics has already been presented as part of Table 1, and is a useful point at which to begin our comparisons. As the end point for extraction becomes more hydrophobic (aq→blob→phase), the extraction equivalent of a permeability coefficient increases. Interestingly, the permeation coefficients indicate that the slowest step of the process tends to be the extraction from the bilayer into solution, since the extraction P when an organic phase is present is lower than the permeability coefficient for bilayer crossing. This has implications for the overall apparent permeability (Eq. 4), where Papparent ≅ Pol→out if Pol→out

Membrane Permeability of Fatty Acyl Compounds ... - ACS Publications

Recommend Documents