
==== Front
J Phys Chem B
J Phys Chem B
jp
jpcbfk
The Journal of Physical Chemistry. B
1520-6106
1520-5207
American Chemical Society

39240094
10.1021/acs.jpcb.4c03835
Article
Membrane Stabilization of Helical Previtamin D Conformers as Possible Enhancement of Vitamin D Photoproduction
Smith Adam C. †⊥
Plazola Matthew †⊥
Hudson Phillip S. ‡§#
https://orcid.org/0000-0002-0640-0297
Tapavicza Enrico *†∥
† Department of Chemistry and Biochemistry, California State University Long Beach, 1250 Bellflower Boulevard, Long Beach, California 90840, United States
‡ Laboratory of Computational Biology, National Institutes of Health, National Heart, Lung and Blood Institute, 12 South Drive, Rm 3053, Bethesda, Maryland 20892-5690, United States
§ Department of Chemistry, University of South Florida, 4202 East Fowler Avenue, CHE205, Tampa, Florida 33620-5250, United States
∥ Faculty of Chemistry and Pharmacy, Institute of Physical and Theoretical Chemistry, University of Regensburg, Universitätsstraße 31, 93040 Regensburg, Germany
* Email: enrico.tapavicza@csulb.edu.
06 09 2024
19 09 2024
128 37 89568965
08 06 2024
30 08 2024
29 08 2024
© 2024 The Authors. Published by American Chemical Society
2024
The Authors
https://creativecommons.org/licenses/by-nc-nd/4.0/ Permits non-commercial access and re-use, provided that author attribution and integrity are maintained; but does not permit creation of adaptations or other derivative works (https://creativecommons.org/licenses/by-nc-nd/4.0/).

Photoinduced vitamin D formation occurs 10–15-fold faster in phospholipid bilayers (PLB) than in isotropic solution. It has been hypothesized that amphipatic interactions of the PLB with the rotationally flexible previtamin D (Pre) stabilize its helical conformers, enhancing thermal intramolecular [1,7]-hydrogen transfer, forming vitamin D. To test this hypothesis, we carried out molecular dynamics (MD) simulations of Pre in a PLB composed of dipalmitoylphosphatidylcholine (DPPC). We designed a classical force field capable of accurately describing the equilibrium composition of Pre conformers. Using adaptive biasing force MD simulations, we determined the free energy of Pre conformers in isotropic environments (hexane and gas-phase) and in the anisotropic environment of a DPPC PLB. We find a total increase of 25.5% of the population of both helical conformers (+20.5% g+Zg+ and +5% g–Zg−) in DPPC compared to hexane. In view of ab initio simulations, showing that hydrogen transfer occurs in both helical conformers, our study strongly suggests the validity of the initial hypothesis. Regarding the amphipatic interactions of Pre with the PLB, we find that, similar to cholesterol (Chol) and 7-dehydrocholesterol (7-DHC), Pre entertains hydrogen bonds mainly to the carbonyl groups of DPPC and, to a lesser extent, with phosphate oxygen atoms and rarely to water molecules at the interface. We further report order parameters of the Pre/DPPC system, which are slightly smaller than those for Chol/DPPC and 7-DHC/DPPC, but larger than for pure DPPC. This indicates a loss in membrane viscosity upon photochemical ring-opening of 7-DHC to form Pre.

National Science Foundation 10.13039/100000001 CHE-1464946 National Institute of General Medical Sciences 10.13039/100000057 UL1GM11897G- 02 National Institute of General Medical Sciences 10.13039/100000057 TL4GM118980 National Institute of General Medical Sciences 10.13039/100000057 RL5GM118978 National Institute of General Medical Sciences 10.13039/100000057 R15GM126524 National Institute of General Medical Sciences 10.13039/100000057 R01GM129519 National Institute of General Medical Sciences 10.13039/100000057 1 R16GM149410-01 National Heart, Lung, and Blood Institute 10.13039/100000050 NA document-id-old-9jp4c03835
document-id-new-14jp4c03835
ccc-price
==== Body
pmcIntroduction

Vitamin D is an important prohormone responsible for various regulatory processes in vertebrates.1 The main source of vitamin D for the majority of the world population is natural, sun-mediated photosynthesis, which occurs in the membrane of epidermal skin cells. Due to the lack of sun exposure in modern indoor lifestyles, vitamin D deficiency is still a major public health concern, particularly for dark-pigmented individuals living in northern latitudes, where UVB light is limited in winter.2 Vitamin D has recently gained renewed interest as a possible mitigating factor in the severity of infectious diseases.3

Besides its impact on public health, natural vitamin D photosynthesis serves as a textbook paradigm for both a biological example of an electrocyclic photochemical reaction (Figure 1) and a photochemical reaction that occurs within a biological cell membrane.4 Factors such as involvement of nuclear tunneling effects in the thermal intramolecular [1,7]-hydrogen shift (Figure 1) converting previtamin D (Pre) to vitamin D add further complexity in understanding this system.5 While the membrane has been found to have subtle effects on the excited-state ring-opening dynamics of provitamin D,6 it plays a crucial regulatory role in natural vitamin D formation, affecting both photochemical and thermal reactions involved in the formation pathway. On one hand, the membrane and the specific spectral distribution of sunlight prevent vitamin D overproduction upon prolonged irradiation.7−9 This effect, due to the conformationally controlled photochemistry of vitamin D seco-steroids, has been studied extensively both experimentally8 and by quantum mechanics-based molecular dynamics (MD) simulations.9−16 On the other hand, the membrane also has the ability to strongly increase the rate of vitamin D production up to 10–15-fold compared to isotropic solution.17−19 This speed up is possibly due to enhancing the rate of the last formation step, the thermal sigmatropic hydrogen transfer reaction, which is only possible if carbons C19 and C9 are sufficiently close to each other. As possible cause for this enhancement, Holick and co-workers18,20 suggested that the amphipathic interactions between the membrane and Pre shift the conformational equilibrium of previtamin D toward its so-called helical rotational isomers g+Zg+ and g–Zg–, depicted in Figure 2.20 The increased hydrogen transfer rate from C19 to C9 in g+Zg+ and g–Zg– seems plausible from steric considerations, as hydrogen donor C19 and acceptor C9 are in close proximity in these conformers.

Figure 1 Formation of vitamin D by photoinduced electrocyclic ring-opening of provitamin D (Pro) followed by [1,7]-sigmatropic hydrogen shift of previtamin D. Atom numbers are indicated. The second step involves a hydrogen transfer from C19 to C9. Dihedral angels ϕ1 and ϕ2 are defined by atoms C10–C5–C6–C7 and C6–C7–C8–C9, respectively.

Figure 2 Density functional theory (DFT) potential energy (kcal/mol) surface as a function of the dihedral angles ψ1 and ψ2. ψ1 and ψ2 are related to ϕ1 and ϕ2 by a shift of approximately 180 deg (see Supporting Information, Figure S2). Ground state equilibrium structures of the conformers and their location in dihedral space are indicated together with the (C19)H–C9 distance in Å (red). The hexatriene units of each conformer are color highlighted. In the labels for the conformers g stands for gauche and t stands for trans, indicating the single-bond conformation of the central double bond in Z conformation. + and – indicated the sign of the associated dihedral angle.

However, to date, no evidence for this mechanism has been provided. In addition, this hypothesis was challenged by Meana-Pañeda and Fernández-Ramos,5 who, suggested an alternative explanation: a stabilization of the formed vitamin D by the membrane. While the two explanations are not mutually exclusive, our simulations confirm the first hypothesis, showing a statistically significant increase in the helical conformer populations by the phospholipid bilayers (PLB) compared to gas-phase and isotropic hexane solutions.

To investigate the effect of the anisotropic environment of the PLB on the conformer equilibrium of Pre, we carried out MD simulations of Pre in a PLB composed of dipalmitoylphosphatidylcholine (DPPC) and in the isotropic environments of hexane solution and the gas-phase. To gauge the effect of the Pre on the membrane properties, we also studied the analogous systems, composed of cholesterol (Chol) in DPPC, 7-dehydrocholesterol (7-DHC) in DPPC, and pure DPPC for comparison. While the composition of epidermal skin cells is complex, DPPC is considered to be a valid model for epidermal skin cell membrane and has been used to mimic mammalian cell membrane.21,22 Regardless of the validity of using DPPC to model epidermal skin cell membranes, the accelerating effect of vitamin D production has been observed equally in epidermal skin cells and in DPPC model systems.19,20

Since an accurate description of the dihedral angle conformation of Pre is crucial for our study, we adopted the CMAP method23−25 and parametrized an empirical correction to the torsional potential associated with dihedral angles ϕ1 and ϕ2. Then we computed the free energy landscape in the dihedral space using adaptive biasing force (ABF) MD. Furthermore, to assess which Pre conformers are prone to undergo thermal hydrogen transfer, we also carried out constrained ab initio molecular dynamics (AIMD) simulations with C19–C9 distances constrained to specific values.

Methods and Computational Details

All ab initio calculations were carried out with TURBOMOLE V7.226,27 and employed the PBE approximation to the exchange–correlation functional28 in combination with the D3 dispersion correction29 and the def2-SVP basis set.30 The PBE-D3 approach has been found to yield reasonable structures for medium-sized organic molecules, with errors below 0.5 kcal/mol with respect to reference structures obtained by CCSD(T).31 Typical average errors on reaction barriers due to the small def2-SVP basis set have been found to be about 0.12 kcal/mol with respect to the large aug-cc-pVTZ basis set.32 In view of the approximate nature of the force field (FF) calculations used in MD, we consider PBE-D3/def2-SVP to be an appropriate level for the AIMD simulations and the CMAP parametrization, balancing computational efficiency with accuracy.

Classical MD simulations were carried out using CHARMM33 and NAMD.34

Constrained Ab Initio Molecular Dynamics

To test the ability of the conformers to undergo [1,7]-hydrogen transfer, we investigated g+Zg+, g–Zg–, and g–Zt+, all exhibiting (C19)H–C9 distances below 3.5 Å (Figure 2) using ground-state AIMD based on DFT, employing the PBE-D3/def2-SVP level. To accelerate the occurrence of the [1,7]-hydrogen transfer, we constrained the C19–C9 distances to values ranging from 2 to 3.8 Å. AIMD was carried out with a Nosé-Hover thermostat35,36 with target temperatures ranging from 300 to 1200 K. High-temperature simulations are used to accelerate the observation of the hydrogen transfer. Purpose of the AIMD simulations was simply to test if a given distance allows hydrogen transfer to occur and were not meant to mimic realistic conditions in an experimental setting. A time step of 50 au was used to propagate the nuclear degrees of freedom for a total of 12 ps. Additional details are given in the Supporting Information, Section 1.

Parameterization of the CMAP

Due to its computational efficiency, classical MD employing an empirical FF is the method of choice to study the relatively large system of Pre in a PLB. However, due to the conjugation effects of the π-orbitals of the conjugated double-bonds in the triene unit,37 the potential energy associated with the two dihedral angles ϕ1 and ϕ2 of the central triene moiety of the molecule cannot be described as a sum of two independent torsional potentials, as it is usually done in most empirical FFs. Consequently, employing the additive CHARMM generalized force field (CGenFF)38 to Pre leads to an erroneous description of the conformer composition (see Supporting Information, Section 2). To address this deficiency, we first derived an accurate empirical FF for Pre based on DFT. Specifically, we adopted the CMAP method,23−25 originally developed to accurately describe backbone dihedral angles in proteins. The CMAP method is applied as a correction term to an additive FF, accounting for the mutual dependency of the potential energy terms associated with two interacting dihedral angle conformations. Using DFT results as reference, our parametrized CMAP together with CGenFF38 leads to an accurate description of the potential energy surface in the ϕ1/ϕ2 space (Figure 2), reproducing the existence of six major rotamers previously found by semiempirical QM39 and various ab initio approaches.5,10,40,41 To obtain the CMAP, used together with CGenFF in classical MD, we scanned the two-dimensional dihedral space in steps of 10° and optimized the structures at the PBE-D3/def2-SVP level. For the optimized structures, we computed the CGenFF energy. The CMAP is obtained as the difference between the total DFT energy and the CGenFF energy. For full results and the parameter file, see the Supporting Information, Section 3. When the CMAP is supplied with the CGenFF parameters to NAMD, the CGenFF forces are corrected by the CMAP using a multidimensional spline fitting of the two-dimensional grid data.23,24

System Assembly

The Pre/hexane system was constructed using the solvation function of VMD.42

We constructed PLB systems of DPPC with the steroid molecules Pre, 7-DHC, and Chol embedded in the hydrophobic region. Membranes consisted of 8 steroid molecules and 138 DPPC molecules yielding a steroid concentration of 5.48 mol % in the membrane—a concentration low enough that no significant increase of the order parameter compared to pure DPPC can be expected.43 Prior to membrane assembly, Pre, 7-DHC, and Chol were arranged into octamer clusters, featuring four steroid molecules per leaflet, arranged to minimize steroid–steroid interactions during equilibration, as shown in Figure S4, in the Supporting Information. Steroids were oriented with OH-groups pointing toward the water–DPPC interface and alkyl tails pointing to the inside of the membrane. Lipid bilayers were assembled using CHARMM-GUI membrane builder and the steroid clusters were incorporated using the replacement method.44,45 Equivalent membranes were also constructed with pure DPPC (146 molecules). Membrane systems contained approximately 6200 water molecules.

While initial structures and FF parameters for Chol were provided by the CHARMM-GUI, initial structures and FF parameters for Pre and 7-DHC were obtained on the basis of structure data format files from the PubChem database.46 Topology and parameter files for Pre and 7-DHC were generated by CHARMM-GUI using CGenFF during membrane assembly.47 Parameter files for DPPC lipid molecules were provided by CHARMM-GUI, employing the CHARMM C36 FF.48 The TIP-3P model was used for water molecules.49

Constant Pressure/Temperature Molecular Dynamics

We used the default input generated by the CHARMM-GUI membrane builder for equilibration and production, as previously described.45 Constant pressure/temperature (CPT) MD was carried out at 1 atm and 323.15 K (50 °C) using CHARMM.33 Additional simulation details are given in the Supporting Information, Section 4. To obtain statistics, six separate simulations of each system with varied initial conditions were propagated using CPT dynamics for 15 ns using a time step of 2 fs. To obtain the order parameter of the membrane systems, we used the order parameter tool from the MEMBPLUGIN50 extension of VMD,42 as previously described.51 The segmental order parameter SCD (eq 1) reflects the viscosity of a bilayer by measuring the degree of isotropic motion in the lipid tails, Sn-1 and Sn-2, depicted in Figure 5. Values of SCD that are close to 0.5 indicate that the lipids are packed tightly together and motion of the tails is inhibited. Values closer to 0 are related to a greater amount of motion in the tails, according to1

where θi is the angle between the C–H bond of the ith carbon atom and the bilayer normal, brackets indicate this value is averaged over the ensemble.

Density profiles and pair distribution functions were generated using VMD.42

Adaptive Force Biasing Molecular Dynamics of Previtamin D

ABF simulations were started from the previously equilibrated systems at two temperatures 300 and 323 K. For each system (gas-phase, hexane, DPPC), at each of the two temperatures, we generated 24 independent initial conditions for ABF simulations by CPT MD of varying simulation times.

Using the Colvars module52 of the NAMD package34 collective variables ϕ1 and ϕ2 were defined according to Figure 1. In the ABF method, the instantaneous force along the collective variables is recorded for bins of 5 × 5°. After having collected the forces for 100 MD steps in each bin, the time average of the force along the collective variables is computed and opposed by a biasing force. This leads to a flattening of the potential energy surface along the collective variables, allowing for increased sampling in the space of the collective variables. After completion of the ABF simulation, the integration of the biasing force along the collective variables yields the free energy profile.53 Total simulation time for each of the 24 ABF simulations of each system amounted to 40 ns.

For each of the 24 independent simulations, we determined the free energy and the corresponding density distributions. From the average free energy landscape, we determined the highest contour that surrounded one particular free energy well. These contours were used to integrate the total population of each conformer. Conformer populations of the 24 independent simulations were used to compute the mean and standard deviation.

Results and Discussion

Constrained Ab Initio Molecular Dynamics

We investigated the conformers that exhibit (C19)H–C9 distances smaller than 3.5 Å (Figure 2) with respect to their ability to undergo hydrogen transfer. The AIMD simulations show (see Supporting Information, Tables S1–S6) that both helical conformers g+Zg+ and g–Zg–, but also g–Zt+ are able to undergo hydrogen transfer if the C19–C9 distance is constrained to distances of 2.6 Å and below. To the best of our knowledge, the importance of the g–Zt+ conformer with regard to the hydrogen transfer reaction has not been thoroughly discussed in previous literature.5,20 However, as we will demonstrate, due to its low population, the contribution of g–Zt+ conformer to the enhancement of the hydrogen transfer can be neglected, and the main contribution arises from the increase of the helical conformer populations of the helical conformers.

Constant Pressure/Temperature Molecular Dynamics

Analyzing the trajectories of the CPT MD simulations of the steroid/DPPC systems, we find that the density profile across the bilayer of Pre/DPPC is qualitatively very similar to that of Chol/DPPC (Figure 3), which we find consistent with previous simulations.43 At more thorough inspection, we notice that the water density is more penetrating toward the bilayer center in the Pre/DPPC system than in the Chol/DPPC system. The water density profile of pure DPPC is more similar to Pre/DPPC than to Chol/DPPC and exhibits overall the most water-penetrating characteristic. This is consistent with previous studies that found the free energy barrier of water permeation to be smaller in pure DPPC than in Chol/DPPC.54,55 The density of individual functional groups between Pre/DPPC and Chol/DPPC (Figure 4) also exhibits an overall qualitative agreement between the two systems. The only noticeable difference lies in the density profile of the OH group of Pre, which has its maximum slightly closer to the bilayer center than the density distribution of the OH group of Chol, although we cannot rule out that this difference is due to a sampling error.

Figure 3 Electron density profile perpendicular to the bilayer of Pre/DPPC (solid), Chol/DPPC (dotted), and pure DPPC (dashed).

Figure 4 Electron density profile perpendicular to the bilayer of Pre/DPPC (solid), Chol/DPPC (dotted).

Regarding the average Sn-1/Sn-2 order parameter SCD of DPPC and of the steroid/DPPC systems (Figure 5), we note that within our simulations (solid lines) Chol/DPPC has the highest order parameter, followed by 7-DHC/DPPC. Consistent with our result, the membrane ordering effect of 7-DHC has previously been found to be slightly smaller than the one of Chol.56 In another study their effects were found to be practically indistinguishable.57

Figure 5 Upper panel: Lewis structure of DPPC. Carbon atom numbers are indicated for the tails Sn-1 and Sn-2 tails. Lower panel: order parameter averaged over Sn-1 and Sn-2 chains as a function of the carbon atom index indicated above. Results from our simulations are shown by circles and solid lines. Experimental data and previous simulations are indicated by the crosses and dotted lines. Experimental NMR data for DPPC and Chol/DPPC are taken from refs (58 and 63), respectively. Previous data from MD simulations is taken from ref (62).

The Pre/DPPC systems (with CMAP and without CMAP) show similar order parameters, both lower than Chol/DPPC and 7-DHC/DPPC, but higher than pure DPPC. The CMAP correction, however, slightly increases the order parameter in the Pre/DPPC system. This might be caused by increasing the amount of helical conformations, facilitating the packing in DPPC due to their flat, Chol-like structure. Interestingly, the deviation in the order parameter between the simulations with CMAP and without CMAP for Pre/DPPC is stronger toward the inside of the bilayer (carbon index larger than 7). This is possibly due to the orthogonal dihedral angle conformations of Pre (see Supporting Information, Section 2) if the additive FF is not corrected, disturbing the packing of DPPC.

NMR-derived order parameters58,59 for Chol/DPPC and pure DPPC show the same trend of increasing order parameter by addition of Chol to the membrane, an effect that is commonly known.43,60,61 In both cases, however, experimental order parameters are lower than their counterparts in our simulations. The spread of the NMR-derived order parameters between Chol/DPPC and DPPC, is comparable to the spread in our simulations. Previous simulations at 323 K by Smondyrev and Berkowitz of pure DPPC using an AMBER-like united atom FF predicted lower order parameter than our simulations and are in good agreement to the experimental results.62

Interestingly, we notice that Pre is relatively flexible in the membrane, adopting a range of different orientations, as seen in the simulation snapshot (Figure 6) and as quantified by the distribution of cos(θ)-values (red, Figure 7), where θ is the angle between the z-axis and the vector formed by atoms C9 and C16. Besides, the orientations parallel and antiparallel to the z-axis, indicated by values +1 and −1, respectively, Pre can also adopt a horizontal orientation, as indicated by the values close to zero. With respect to this property, Pre behaves similar to Chol, which can also adopt various orientations in the DPPC bilayer (blue, Figure 7).

Figure 6 Snapshot from a MD simulation. In the center, a Pre molecule in g+Zg+ conformation is highlighted by the thick stick model, with the hexatriene unit highlighted in red. The shortest (C19)H–C9 distance is indicated (3.13 Å). In addition, a hydrogen bond from the hydroxy group of Pre to a water molecule at the polar bilayer/water interface is indicated (1.75 Å). Noncarbon atoms of the polar DPPC head groups are indicated by the ball models (blue: nitrogen, red oxygen, orange: phosphorus). DPPC carbon tails are represented by the thin gray line model. Additional Pre molecules dissolved in the PLB are indicated by the thin stick model.

Figure 7 Distribution of cos(θ)-values for Chol (blue) and Pre (red). θ is defined as the angle between the z-axis (bilayer normal) and the vector formed by atoms C9 and C16.

To investigate the hydrogen bonding partners of the OH group of Pre, we extracted the oxygen–oxygen radial distribution function (RDF) of the Pre–OH oxygen with other oxygen atoms in the system (Figure 8). An O–O distance of approximately 2.7 Å indicates a shared H-bond between two oxygen atoms.64 We find that the Pre–OH oxygen most likely forms H-bonds with carbonyl oxygens O22 and O32 of chains Sn-1 and Sn-2, respectively. The second most probable bonding partner are the free oxygen atoms of the phosphate groups, O13 and O14. Less likely, with similar probability, we find H-bonds between water oxygen (OH2) and phosphate ester oxygen O11. With this H-bonding trend, Pre exhibits the same qualitative behavior as previously found for Chol/DPPC.43

Figure 8 RDFs between the oxygen of the hydroxyl group of Pre and different oxygen atoms of DPPC and water. The RDFs have been computed with a distance bin size of 0.05 Å. For a discussion of the long-range behavior of the RDFs, please see the Supporting Information, Section 4.2.

Adaptive Biasing Force Molecular Dynamics of Previtamin D

To test the hypothesis of an increase in the helical conformer populations in the presence of the PLB, we will compare the populations of the different Pre rotamers, depicted in Figure 2 in the anisotropic environment of the PLB to the conformer populations in isotropic environment using MD.

Using CGenFF plus the derived CMAP, we applied the ABF method53,65 to sample the free energy profile in the dihedral space of Pre via thermodynamic integration. This was done in isotropic environments (gas-phase and hexane solution) and in the anisotropic environment of DPPC. For each system, we carried out 24 independent ABF simulations of 40 ns, each with distinct initial conditions. The average free energy landscapes in the ϕ1/ϕ2 space (first column, Figure 9) appear to be qualitatively similar for all systems, reflecting the underlying DFT potential energy surface (Figure 2) to some extent. However, the resulting free energy profiles differ slightly in magnitude and extent of their energetic wells. These differences become more obvious when inspecting the density distributions (second column, Figure 9) resulting from the free energy profile: we recognize an increased density for the g+Zg+ and g–Zg– conformers in DPPC (panel H, red and green contours, respectively) compared to the gas-phase (panel C) and hexane solution (panel E). Integrating the density within each contour, we observe a population of 56.8% for the g+Zg+ conformer in DPPC (panel I), compared to 33.9% (panel C) and 36.3% (panel F) in the gas-phase and hexane, respectively. Also for the g–Zg– conformer, the population in DPPC is increased (6.86%) compared to the gas-phase and hexane solution (1.28 and 1.82%, respectively). However, the standard deviations of the conformer populations over the 24 individual simulations are in the same range as the population differences in anisotropic and isotropic environments (Figure 9, third column; Supporting Information, Table S7). To test for statistical significance of the population differences, we applied a two-sample independent t-test66 (for complete results see Supporting Information, Table S9). The two-sample t-test shows that the relevant population differences are statistically significant: within a 95%—confidence interval (CI), the g+Zg+ population is at least 7.5% larger in DPPC than in hexane and at most 33.5% larger. For g–Zg–, the 95%—CI predicts the g–Zg– population to be between 1.0 and 9.1% larger in DPPC than in hexane. Similar, but slightly larger, increases are observed for the comparison of the DPPC populations to the gas-phase populations. In contrast, the population differences between the gas-phase and the hexane solution are statistically insignificant. Additionally, we find a shallower energy minimum of t+Zg+ and a lower barrier to g+Zg+ in DPPC, compared to the simulations in isotropic environment.

Figure 9 Left column: Helmholtz free energy (kcal/mol) of Pre as a function of ϕ1 and ϕ2 as obtained from ABF simulations at 300 K. Middle column: probability distribution (1/deg2) in the ϕ1/ϕ2-space, obtained from free energy surface. Right column: conformer populations, obtained by integration of the dashed contour lines of the probability distribution. Error bars indicate the standard deviation over 24 independent ABF simulations of each 40 ns. Upper (A–C), middle (D–F), and lower (G–I) rows refer to simulations of Pre in the gas-phase, in hexane solution, and in DPPC, respectively. Color code for middle and right column: red: g+Zg+, green: g–Zg–, gray: t+Zg–, cyan: t–Zg+, magenta: t–Zg–, yellow: g–Zt+.

The same trends are also observed in the results for simulations at 323 K (see Supporting Information, Figure S5, Tables S8 and S10). According to the t-test, the differences between the simulations at 300 and 323 K are not statistically significant.

Regarding the population of the g–Zt+ conformer, which, according to our AIMD simulations, can undergo [1,7]-hydrogen transfer, we observed a reduction from less than 1% in isotropic environments to zero in the anisotropic environment, within the accuracy of our calculations. The contribution of this conformer to the increased vitamin D formation rate can therefore be neglected.

Conclusion

In summary, we conclude that Pre incorporates into the DPPC membrane with structural similarity to Chol and similar hydrogen bonding pattern.

Incorporation of Pre into DPPC leads to a slightly lower order parameter than Chol/DPPC and 7-DHC/DPPC, which can possibly be attributed to the higher degree of flexibility due to the open-ring structure of a seco-steroid compared to the closed ring structure of steroids.

The anisotropic environment of the DPPC bilayer in combination with the amphipatic character of Pre leads to an average increase of about 25.5% of both the helical conformers compared to hexane solution. The dominating increase of about 20.5% is found for g+Zg+ conformer, but also g–Zg– contributes about 5% to the increase in helical conformers. This refines the previous suggestion,67 in which the g–Zg– conformer is merely a transient species that is populated right after provitamin D ring-opening, but quickly diffuses into other conformers. Our simulations here show that it persists in the equilibrium with a population of about 7% in DPPC.

However, our simulations are not able to quantitatively predict the 10–15-fold increase of the vitamin D production rate. Previously, sophisticated ab initio calculations, including nuclear tunneling effects,5 predicted hydrogen transfer rates for helical conformers also to be sensitive for other conformational parameters, such as the position of the OH-group (axial vs equatorial) and the ring conformations of Pre (boot, seat, twisted boat). For each helical conformer, 16 distinct subconformers, characterized by these additional geometric parameters, have been identified, exhibiting different H-transfer rates. Since the collective variables in our ABF simulations only considered dihedral parameters ϕ1 and ϕ2 and did not consider these additional conformational parameters, we are not able to assess the relative free energies with respect to these conformational parameters, preventing quantitative prediction of the resulting hydrogen transfer rates in different media. Nevertheless, in view of the AIMD results (Supporting Information, Section 1), showing that hydrogen transfer is possible in g+Zg+ and g–Zg–, our study strongly suggests that the amphipatic interactions of Pre with the PLB (Figure 6) will accelerate the vitamin D formation by stabilizing the helical conformers, as hypothesized by Holick and co-workers in 1993.18

Supporting Information Available

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpcb.4c03835.Complete results from AIMD, results from MD simulations without CMAP, details of the CMAP parametrization, the full CMAP parameter file, details and additional results of the CPMD dynamics, additional details of the ABF simulations, results of the ABF simulations at 323 K, results of the two-sample independent t-test (PDF)

Supplementary Material

jp4c03835_si_001.pdf

Author Present Address

# Medicine Design, Pfizer Inc., Cambridge Massachusetts 02139, United States

Author Contributions

⊥ A.C.S. and M.P. contributed equally to this work.

The authors declare no competing financial interest.

Acknowledgments

We want to thank Stefan Boresch for his help in generating the CMAP corrections and Wataru Shinoda for useful discussions. Research reported in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health (NIH) under award numbers R15GM126524, 1 R16GM149410-01, UL1GM11897G-02, TL4GM118980, RL5GM118978. P.S.H. acknowledges funding support from the Intramural Research Program of the NIH, NHLBI. P.S.H. also acknowledges support from H. Lee Woodcock via the National Science Foundation (CHE-1464946) and NIGMS of the National Institutes of Health (R01GM129519). The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH. We also acknowledge technical support from the Division of Information Technology of CSULB.
==== Refs
References

Bikle D. ; Christakos S. New aspects of vitamin D metabolism and action—Addressing the skin as source and target. Nat. Rev. Endocrinol. 2020, 16 , 234–252. 10.1038/s41574-019-0312-5.32029884
Clemens T. ; Henderson S. ; Adams J. ; Holick M. Increased skin pigment reduces the capacity of skin to synthesise vitamin D3. Lancet 1982, 319 , 74–76. 10.1016/S0140-6736(82)90214-8.
Martineau A. R. ; Forouhi N. G. Vitamin D for COVID-19: a case to answer?. Lancet Diabetes Endocrinol. 2020, 8 , 735–736. 10.1016/S2213-8587(20)30268-0.32758429
Holick M. F. ; MacLaughlin J. A. ; Clark M. B. ; Holick S. A. ; Potts J. T. ; Anderson R. R. ; Blank I. H. ; Parrish J. A. ; Elias P. Photosynthesis of previtamin D3 in human skin and the physiologic consequences. Science 1980, 210 , 203–205. 10.1126/science.6251551.6251551
Meana-Pañeda R. ; Fernández-Ramos A. Tunneling and conformational flexibility play critical roles in the isomerization mechanism of vitamin D. J. Am. Chem. Soc. 2012, 134 , 346–354. 10.1021/ja2077075.22118472
Sofferman D. L. ; Konar A. ; Spears K. G. ; Sension R. J. Ultrafast excited state dynamics of provitamin D3 and analogs in solution and in lipid bilayers. J. Chem. Phys. 2021, 154 , 094309 10.1063/5.0041375.33685160
Holick M. F. ; MacLaughlin J. A. ; Doppelt S. H. Regulation of cutaneous previtamin D3 photosynthesis in man: skin pigment is not an essential regulator. Science 1981, 211 , 590–593. 10.1126/science.6256855.6256855
MacLaughlin J. ; Anderson R. ; Holick M. Spectral character of sunlight modulates photosynthesis of previtamin D3 and its photoisomers in human skin. Science 1982, 216 , 1001–1003. 10.1126/science.6281884.6281884
Thompson T. ; Tapavicza E. First-Principles Prediction of Wavelength-Dependent Product Quantum Yields. J. Phys. Chem. Lett. 2018, 9 , 4758–4764. 10.1021/acs.jpclett.8b02048.30048134
Tapavicza E. ; Meyer A. M. ; Furche F. Unravelling the details of vitamin D photosynthesis by non-adiabatic molecular dynamics simulations. Phys. Chem. Chem. Phys. 2011, 13 , 20986 10.1039/c1cp21292c.22020179
Tapavicza E. ; Bellchambers G. D. ; Vincent J. C. ; Furche F. Ab initio non-adiabatic molecular dynamics. Phys. Chem. Chem. Phys. 2013, 15 , 18336–18348. 10.1039/c3cp51514a.24068257
Cisneros C. ; Thompson T. ; Baluyot N. ; Smith A. C. ; Tapavicza E. The role of tachysterol in vitamin D photosynthesis - a non-adiabatic molecular dynamics study. Phys. Chem. Chem. Phys. 2017, 19 , 5763–5777. 10.1039/C6CP08064B.28105477
Tapavicza E. ; Thompson T. ; Redd K. ; Kim D. Tuning the photoreactivity of Z-hexatriene photoswitches by substituents – a non-adiabatic molecular dynamics study. Phys. Chem. Chem. Phys. 2018, 20 , 24807–24820. 10.1039/C8CP05181J.30229769
Schalk O. ; Tapavicza E. Photochemistry. In ACS In Focus; American Chemical Society, 2021.
Tapavicza E. In Time-Dependent Density Functional Theory: Nonadiabatic Molecular Dynamics; Zhu C. , Ed.; Jenny Stanford Publishing, 2022; pp 141–197.
Tapavicza E. ; Reutershan T. ; Thompson T. Ab initio simulation of the ultrafast circular dichroism spectrum of provitamin D ring-opening. J. Phys. Chem. Lett. 2023, 14 , 5061–5068. 10.1021/acs.jpclett.3c00862.37227143
Yamamoto J. K. ; Borch R. F. Photoconversion of 7-dehydrocholesterol to vitamin D3 in synthetic phospholipid bilayers. Biochemistry 1985, 24 , 3338–3344. 10.1021/bi00334a039.2992580
Tian X. Q. ; Chen T. C. ; Matsuoka L. Y. ; Wortsman J. ; Holick M. F. Kinetic and thermodynamic studies of the conversion of previtamin D3 to vitamin D3 in human skin. J. Biol. Chem. 1993, 268 , 14888–14892. 10.1016/S0021-9258(18)82416-4.8392061
Holick M. F. ; Tian X. Q. ; Allen M. Evolutionary importance for the membrane enhancement of the production of vitamin D3 in the skin of poikilothermic animals. Proc. Natl. Acad. Sci. U.S.A. 1995, 92 , 3124–3126. 10.1073/pnas.92.8.3124.7724526
Tian X. Q. ; Holick M. F. A liposomal model that mimics the cutaneous production of vitamin D3. J. Biol. Chem. 1999, 274 , 4174–4179. 10.1074/jbc.274.7.4174.9933613
Van Meer G. ; Voelker D. R. ; Feigenson G. W. Membrane lipids: where they are and how they behave. Nat. Rev. Mol. Cell Biol. 2008, 9 , 112–124. 10.1038/nrm2330.18216768
Moniz T. ; Costa Lima S. A. ; Reis S. Human skin models: From healthy to disease-mimetic systems; characteristics and applications. Br. J. Pharmacol. 2020, 177 , 4314–4329. 10.1111/bph.15184.32608012
MacKerell A. D. ; Feig M. ; Brooks C. L. Improved treatment of the protein backbone in empirical force fields. J. Am. Chem. Soc. 2004, 126 , 698–699. 10.1021/ja036959e.14733527
Mackerell A. D. ; Feig M. ; Brooks C. L. Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 2004, 25 , 1400–1415. 10.1002/jcc.20065.15185334
Buck M. ; Bouguet-Bonnet S. ; Pastor R. W. ; MacKerell A. D. Importance of the CMAP correction to the CHARMM22 protein force field: dynamics of hen lysozyme. Biophys. J. 2006, 90 , L36–L38. 10.1529/biophysj.105.078154.16361340
Balasubramani S. G. ; Chen G. P. ; Coriani S. ; Diedenhofen M. ; Frank M. S. ; Franzke Y. J. ; Furche F. ; Grotjahn R. ; Harding M. E. ; Hättig C. ; et al. TURBOMOLE: Modular program suite for ab initio quantum-chemical and condensed-matter simulations. J. Chem. Phys. 2020, 152 , 184107 10.1063/5.0004635.32414256
Franzke Y. J. ; Holzer C. ; Andersen J. H. ; Begušić T. ; Bruder F. ; Coriani S. ; Della Sala F. ; Fabiano E. ; Fedotov D. A. ; Fürst S. ; et al. TURBOMOLE: Today and tomorrow. J. Chem. Theory Comput. 2023, 19 , 6859–6890. 10.1021/acs.jctc.3c00347.37382508
Perdew J. P. ; Burke K. ; Ernzerhof M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77 , 3865–3868. 10.1103/PhysRevLett.77.3865.10062328
Grimme S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27 , 1787–1799. 10.1002/jcc.20495.16955487
Schäfer A. ; Horn H. ; Ahlrichs R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. J. Chem. Phys. 1992, 97 , 2571–2577. 10.1063/1.463096.
Vuckovic S. ; Burke K. Quantifying and understanding errors in molecular geometries. J. Phys. Chem. Lett. 2020, 11 , 9957–9964. 10.1021/acs.jpclett.0c03034.33170683
Wang M. ; He X. ; Taylor M. ; Lorpaiboon W. ; Mun H. ; Ho J. Molecular geometries and vibrational contributions to reaction thermochemistry are surprisingly insensitive to the choice of basis sets. J. Chem. Theory Comput. 2023, 19 , 5036–5046. 10.1021/acs.jctc.3c00388.37463146
Brooks B. R. ; Brooks C. L. ; Mackerell A. D. ; Nilsson L. ; Petrella R. J. ; Roux B. ; Won Y. ; Archontis G. ; Bartels C. ; Boresch S. ; et al. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30 , 1545–1614. 10.1002/jcc.21287.19444816
Phillips J. C. ; Hardy D. J. ; Maia J. D. C. ; Stone J. E. ; Ribeiro J. V. ; Bernardi R. C. ; Buch R. ; Fiorin G. ; Hénin J. ; Jiang W. ; et al. Scalable molecular dynamics on CPU and GPU architectures with NAMD. J. Chem. Phys. 2020, 153 , 044130 10.1063/5.0014475.32752662
Nosé S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81 , 511–519. 10.1063/1.447334.
Hoover W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985, 31 , 1695–1697. 10.1103/PhysRevA.31.1695.
DuBay K. H. ; Hall M. L. ; Hughes T. F. ; Wu C. ; Reichman D. R. ; Friesner R. A. Accurate force field development for modeling conjugated polymers. J. Chem. Theory Comput. 2012, 8 , 4556–4569. 10.1021/ct300175w.26605615
Vanommeslaeghe K. ; Hatcher E. ; Acharya C. ; Kundu S. ; Zhong S. ; Shim J. ; Darian E. ; Guvench O. ; Lopes P. ; Vorobyov I. ; et al. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31 , 671–690. 10.1002/jcc.21367.19575467
Dmitrenko O. ; Frederick J. H. ; Reischl W. Previtamin D conformations and the wavelength-dependent photoconversions of previtamin D. J. Photochem. Photobiol., A 2001, 139 , 125–131. 10.1016/S1010-6030(00)00425-1.
Bayda M. ; Redwood C. E. ; Gupta S. ; Dmitrenko O. ; Saltiel J. Lumisterol to tachysterol photoisomerization in EPA glass at 77 K. A comparative study. J. Phys. Chem. A 2017, 121 , 2331–2342. 10.1021/acs.jpca.6b12843.28234492
Ferro-Costas D. ; Sánchez-Murcia P. A. ; Fernández-Ramos A. Unraveling the Catalytic Mechanism of β-Cyclodextrin in the Vitamin D Formation. J. Chem. Inf. Model. 2024, 64 , 3865–3873. 10.1021/acs.jcim.3c02049.38598310
Humphrey W. ; Dalke A. ; Schulten K. VMD – Visual Molecular Dynamics. J. Mol. Graphics 1996, 14 , 33–38. 10.1016/0263-7855(96)00018-5.
Tu K. ; Klein M. L. ; Tobias D. J. Constant-pressure molecular dynamics investigation of cholesterol effects in a dipalmitoylphosphatidylcholine bilayer. Biophys. J. 1998, 75 , 2147–2156. 10.1016/S0006-3495(98)77657-X.9788908
Jo S. ; Kim T. ; Iyer V. G. ; Im W. CHARMM-GUI: a web-based graphical user interface for CHARMM. J. Comput. Chem. 2008, 29 , 1859–1865. 10.1002/jcc.20945.18351591
Wu E. L. ; Cheng X. ; Jo S. ; Rui H. ; Song K. C. ; Dávila-Contreras E. M. ; Qi Y. ; Lee J. ; Monje-Galvan V. ; Venable R. M. ; et al. CHARMM-GUI membrane builder toward realistic biological membrane simulations. J. Comput. Chem. 2014, 35 , 1997–2004. 10.1002/jcc.23702.25130509
Bolton E. E. ; Chen J. ; Kim S. ; Han L. ; He S. ; Shi W. ; Simonyan V. ; Sun Y. ; Thiessen P. A. ; Wang J. ; et al. PubChem3D: a new resource for scientists. J. Cheminf. 2011, 3 , 32 10.1186/1758-2946-3-32.
Klauda J. B. ; Venable R. M. ; Freites J. A. ; O’Connor J. W. ; Tobias D. J. ; Mondragon-Ramirez C. ; Vorobyov I. ; MacKerell A. D. ; Pastor R. W. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. J. Phys. Chem. B 2010, 114 , 7830–7843. 10.1021/jp101759q.20496934
Pastor R. ; MacKerell A. Development of the CHARMM force field for lipids. J. Phys. Chem. Lett. 2011, 2 , 1526–1532. 10.1021/jz200167q.21760975
Jorgensen W. L. ; Chandrasekhar J. ; Madura J. D. ; Impey R. W. ; Klein M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79 , 926–935. 10.1063/1.445869.
Guixà-González R. ; Rodriguez-Espigares I. ; Ramírez-Anguita J. M. ; Carrió-Gaspar P. ; Martinez-Seara H. ; Giorgino T. ; Selent J. MEMBPLUGIN: studying membrane complexity in VMD. Bioinformatics 2014, 30 , 1478–1480. 10.1093/bioinformatics/btu037.24451625
Allen P. ; Smith A. C. ; Benedicto V. ; Abdulhasan A. ; Narayanaswami V. ; Tapavicza E. Molecular dynamics simulation of apolipoprotein E3 lipid nanodiscs. Biochim. Biophys. Acta, Biomembr. 2024, 1866 , 184230 10.1016/j.bbamem.2023.184230.37704040
Fiorin G. ; Klein M. L. ; Hénin J. Using collective variables to drive molecular dynamics simulations. Mol. Phys. 2013, 111 , 3345–3362. 10.1080/00268976.2013.813594.
Comer J. ; Gumbart J. C. ; Hénin J. ; Lelièvre T. ; Pohorille A. ; Chipot C. The adaptive biasing force method: Everything you always wanted to know but were afraid to ask. J. Phys. Chem. B 2015, 119 , 1129–1151. 10.1021/jp506633n.25247823
Saito H. ; Shinoda W. Cholesterol effect on water permeability through DPPC and PSM lipid bilayers: a molecular dynamics study. J. Chem. Phys. 2011, 115 , 15241–15250. 10.1021/jp201611p.
Issack B. B. ; Peslherbe G. H. Effects of cholesterol on the thermodynamics and kinetics of passive transport of water through lipid membranes. J. Phys. Chem. B 2015, 119 , 9391–9400. 10.1021/jp510497r.25679811
Róg T. ; Vattulainen I. ; Jansen M. ; Ikonen E. ; Karttunen M. Comparison of cholesterol and its direct precursors along the biosynthetic pathway: effects of cholesterol, desmosterol and 7-dehydrocholesterol on saturated and unsaturated lipid bilayers. J. Chem. Phys. 2008, 129 , 154508 10.1063/1.2996296.19045210
Liu Y. ; Chipot C. ; Shao X. ; Cai W. The effects of 7-dehydrocholesterol on the structural properties of membranes. Phys. Biol. 2011, 8 , 056005 10.1088/1478-3975/8/5/056005.21865621
Seelig A. ; Seelig J. Dynamic structure of fatty acyl chains in a phospholipid bilayer measured by deuterium magnetic resonance. Biochemistry 1974, 13 , 4839–4845. 10.1021/bi00720a024.4371820
Seelig J. ; Niederberger W. Deuterium-labeled lipids as structural probes in liquid crystalline bilayers. Deuterium magnetic resonance study. J. Am. Chem. Soc. 1974, 96 , 2069–2072. 10.1021/ja00814a014.
Berkowitz M. L. Detailed molecular dynamics simulations of model biological membranes containing cholesterol. Biochim. Biophys. Acta, Biomembr. 2009, 1788 , 86–96. 10.1016/j.bbamem.2008.09.009.
Róg T. ; Pasenkiewicz-Gierula M. ; Vattulainen I. ; Karttunen M. Ordering effects of cholesterol and its analogues. Biochim. Biophys. Acta, Biomembr. 2009, 1788 , 97–121. 10.1016/j.bbamem.2008.08.022.
Smondyrev A. M. ; Berkowitz M. L. Molecular dynamics study of Sn-1 and Sn-2 chain conformations in dipalmitoylphosphatidylcholine membranes. J. Chem. Phys. 1999, 110 , 3981–3985. 10.1063/1.478278.
Hofsäß C. ; Lindahl E. ; Edholm O. Molecular dynamics simulations of phospholipid bilayers with cholesterol. Biophys. J. 2003, 84 , 2192–2206. 10.1016/S0006-3495(03)75025-5.12668428
Brown I. D. On the geometry of O–H···O hydrogen bonds. Acta Crystallogr., Sect. A: Found. Adv. 1976, 32 , 24–31. 10.1107/S0567739476000041.
Henin J. ; Fiorin G. ; Chipot C. ; Klein M. L. Exploring multidimensional free energy landscapes using time-dependent biases on collective variables. J. Chem. Theory Comput. 2010, 6 , 35–47. 10.1021/ct9004432.26614317
The MathWorks Inc. MATLAB Version: 9.13.0 (R2022b); The MathWorks Inc., 2023. https://www.mathworks.com/help/stats/ttest2.html.
Sofferman D. L. ; Konar A. ; Mastron J. N. ; Spears K. G. ; Cisneros C. ; Smith A. C. ; Tapavicza E. ; Sension R. J. Probing the Formation and Conformational Relaxation of Previtamin D3 and Analogues in Solution and in Lipid Bilayers. J. Phys. Chem. B 2021, 125 , 10085–10096. 10.1021/acs.jpcb.1c04376.34473504
