Unperturbed hydrocarbon chains and liquid phase bilayer lipid chains: a computer simulation study

In this work, the properties of saturated and unsaturated fatty acid acyl chains 16:0, 18:0, 18:1(n-9)cis, 18:2(n-6)cis, 18:3(n-3)cis, 18:4(n-3)cis, 18:5(n-3)cis, 20:4(n-6)cis, 20:5(n-3)cis and 22:6(n-3)cis in a bilayer liquid crystalline state and similar hydrocarbon chains (with CH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}3 terminal groups instead of C=O groups) in the unperturbed state characterised by a lack of long-range interaction were investigated. The unperturbed hydrocarbon chains were modelled by Monte Carlo simulations at temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T = 303$$\end{document}T=303 K; sixteen fully hydrated homogeneous liquid crystalline phosphatidylcholine bilayers containing these chains were studied by molecular dynamics simulations at the same temperature. To eliminate effects of the simulation parameters, the molecular dynamics and Monte Carlo simulations were carried out using the same structural data and force field coefficients. From these computer simulations, the average distances between terminal carbon atoms of the chains (end-to-end distances) were calculated and compared. The trends in the end-to-end distances obtained for the unperturbed chains were found to be qualitatively similar to those obtained for the same lipid chains in the bilayers. So, for understanding of a number of processes in biological membranes (e.g., changes in fatty acid composition caused by environmental changes such as temperature and pressure), it is possible to use, at least as a first approximation, the relationships between the structure and properties for unperturbed or isolated hydrocarbon chains.


Introduction
Biological membranes are very complex heterogeneous systems composed of various molecules such as lipids, sterols, proteins, carbohydrates, etc. Lipid molecules contain different head groups and a wide variety of acyl chains of fatty acids (FAs; Cook and McMaster 2002;Nelson and Cox 2008;Mouritsen and Bagatolli 2016). The FAs are the fundamental building blocks of all lipids in living matter. The most abundant class of lipids in the biological membranes of animals and plants is phosphatidylcholine (PC). FA acyl chains of PC lipids usually contain 12-24 carbon atoms; the most common chain lengths fall between 14 (or 16) and 22. Most of the FA acyl chains are unsaturated, containing 1-6 double bonds of the cis configuration in different positions; the majority of the double bonds in the tails are methylene-interrupted (i.e., one methylene group is localized between each pair of double bonds; Cook and McMaster 2002;Nelson and Cox 2008;Mouritsen and Bagatolli 2016). It is most common to find chains with an even number of carbon atoms, whereas odd ones are found in rare cases.
Abstract In this work, the properties of saturated and unsaturated fatty acid acyl chains 16:0, 18:0, 18:1(n-9)cis, 18:2(n-6)cis, 18:3(n-3)cis, 18:4(n-3)cis, 18:5(n-3)cis, 20:4(n-6)cis, 20:5(n-3)cis and 22:6(n-3)cis in a bilayer liquid crystalline state and similar hydrocarbon chains (with CH 3 terminal groups instead of C=O groups) in the unperturbed state characterised by a lack of long-range interaction were investigated. The unperturbed hydrocarbon chains were modelled by Monte Carlo simulations at temperature T = 303 K; sixteen fully hydrated homogeneous liquid crystalline phosphatidylcholine bilayers containing these chains were studied by molecular dynamics simulations at the same temperature. To eliminate effects of the simulation parameters, the molecular dynamics and Monte Carlo simulations were carried out using the same structural data and force field coefficients. From these computer simulations, the average distances between terminal carbon atoms of the chains (end-to-end distances) were calculated and compared. The trends in the end-to-end distances obtained for the unperturbed chains were found to The FA chains are often denoted in accordance with 'n-minus ' nomenclature, as N:d(n − j)cis, where N refers to the total number of carbon atoms in the chain, d is the number of the methylene-interrupted double bonds and j is the position of the first double bond, counted from the methyl (CH 3 ) terminus of the chain, with the methyl carbon as number 1. For brevity, below the term (n − j)cis in the notation will be occasionally omitted.
PU FA chains have been linked to a great number of biochemical processes. One of the most notable, observations Table 1 Occurrence of unsaturated and polyunsaturated (PU) fatty acid (FA) chains of lipids in animal and plant membranes Fatty acid chain Occurrence (findings) Refs.
18:1(n-9)cis FA The most abundant monoenoic FA in plant and animal tissues 35-60% of the total FAs of peanut oil acylglycerols Carrin and Carelli (2010) 35-69% of the total FAs of peanut oil acylglycerols Köckritz and Martin (2008) 60% of the total FAs of the oil from 00-quality oilseed rape Wittkop et al. (2009)
Unfortunately, experimental data for different properties of such hydrocarbons or FA acyl chains are scarce or lacking. Computer simulation is nowadays one of the most powerful tools for studying the properties of different molecular systems (Leach 2001;Berendsen 2007;Gould et al. 2007;Landau and Binder 2009;Binder and Heermann 2010;Brooks et al. 2011;Satoh 2011) including lipid membranes, lipids and lipid chains (see, e.g., several reviews on MD simulations of lipid membranes: Bennett and Tieleman 2013;Rabinovich and Lyubartsev 2013;Pluhackova and Böckmann 2015;Baoukina and Tieleman 2016;Bunker et al. 2016;Kirsch and Böckmann 2016;Lyubartsev and Rabinovich 2016;Pasenkiewicz-Gierula et al. 2016;Pöyry and Vattulainen 2016, and other articles published in Special Issue 10 of BBA-Biomembranes, 2016, v.1858 because it allows one to obtain information on an atomistic level. On the other hand, computer simulations of lipid bilayer systems with all possible combinations of chains are still very time-consuming, and, therefore, a different approach to the task would be valuable. To obtain 'structure-property' relationships for different hydrocarbon chains which can be compared with each other, one should use uniform conditions: the same state at the same temperature. The 'unperturbed' state of chain molecules (Flory 1969) was used as the uniform state of hydrocarbon chains in this work. The exact definition of this state is presented below, in the '2.2 Monte Carlo simulations' section. It has been proposed (Flory 1969) that properties of chain molecules in this state correspond to the properties in the bulk amorphous state. Neutron scattering experiments were later carried out (Dettenmaier 1978;Yoon and Flory 1978) and the results substantiated this prediction. On the other hand, biological membranes in a physiological form exist in a liquid crystalline (fluid) state having a relatively high degree of disorder and dynamical behaviour; this state is vital for the proper functioning of membranes.
The aim of the present study was to compare properties of a set of hydrocarbon chains in the unperturbed state unaffected by long-range interactions and for comparison in the liquid crystalline state of lipid bilayers, to assess if these properties are similar to each other.
A structural scheme of the chains considered in the present work is as follows: where a, d, b are integers. The total number of carbons of the chain is N = a + 3d + b + 2. For clear visualization of the connection between structure of the chains and FA acyl chains, these hydrocarbons will be denoted below as alk − N:d(n − j)cis, i.e., similar to the 'n-minus' nomenclature for FAs.
In addition, 16 fully hydrated homogeneous liquid crystalline PC bilayers containing these chains (as FA acyls) were studied by molecular dynamics (MD) simulations at the same temperature (T = 303 K). The MD and MC simulations were carried out using the same force field, to eliminate the effect of the simulation parameters. Both techniques are described below.

Molecular dynamics simulations
Sixteen fully hydrated homogeneous phosphatidylcholine (PC) bilayers were studied by MD simulation in an NPTensemble at temperature T = 303 K and pressure P = 1 bar. The simulation software was the MDynaMix package (Lyubartsev and Laaksonen 2000). The bilayers studied were comprised of one of the PC molecules which contained a saturated sn-1 chain (16:0 or 18:0) and an unsaturated sn-2 chain: 1-palmitoyl-2-oleoyl-sn-glycero-3-PC [16:0/18:1(n-9)cis PC]; 1-stearoyl-2-oleoyl-sn-glycero-3-PC The choice of the lipid set was motivated by the following considerations: (1) the sixteen bilayers listed contain the most important, biologically meaningful types of PU lipids (see "Introduction" section); (2) An inspection of these bilayers under the same conditions allows one to study: 1. the 'double bond number dependence' ('d dependence') of an sn-2 chain of lipid properties over the whole possible range of d from 1 to 5 at fixed N = 18 The presence of the acyl chain 16:0 in the 8 bilayers 16:0/... PC and acyl chain 18:0 in the 8 bilayers 18:0/...PC is sufficient to study geometrical properties of acyls 16:0 and 18:0 by MD simulation. End-to-end distances (between carbons) calculated for each saturated acyl were averaged over eight corresponding PC bilayers; the influence of the position (sn-1 instead of sn-2) of these chains on their end-to-end distances was neglected. Experimental melting temperatures T c of the bilayers were checked to make sure that the temperature T = 303 K in the MD simulations was appropriate. Available published experimental temperatures T c for the lamellar gel to liquidcrystalline phase transition of five of eight studied mixedchain PCs with 16:0 chains in the sn-1 position are gathered in Tables 2 and 3, and six of eight simulated PCs with 18:0 sn-1-chains are gathered in Tables 4 and 5. To the authors' knowledge, no experimental investigations of the melting Table 2 Published temperatures T c of the experimental lamellar gel to liquid-crystalline phase transition of phosphatidylcholines 16:0/18:1(n-9)cis PC DSC differential scanning calorimetry, 2 H NMR deuterium nuclear magnetic resonance, SANS small-angle neutron scattering, Hydrol. meas. hydrolysis measure. T c is the temperature averaged over the gel to liquidcrystalline and liquid-crystalline to gel phase transition temperatures, i.e., the heating and cooling transition temperatures, in all cases when a hysteresis was observed. The T c values for each lipid are presented in order of increasing temperature; values from the lipid database (Koynova and Caffrey 1998) are also presented in the end of the experimental data list. a The existence of metastable forms for the pure 16:0/18:1(n-9)cis PC liposomes was detected in Ref. Lavialle and Levin (1980). The I 2940 /I 2885 peak height intensity ratio as index was used (interchain disorderorder parameter), where I 2940 and I 2885 represent the peak height intensities for 2940-and 2885-cm −1 transitions. b Aqueous 50 wt% ethylene glycol solution. c The I 1100 /I 1130 peak height intensity ratio as index was used (intramolecular gauche-trans isomerization parameter), where I 1100 and I 1130 represent the peak height intensities for 1100-and 1130-cm −1 transitions, respectively Hydrol. meas. Kamp et al. (1975) 270.65 ± 2.4 lipid database Koynova and Caffrey (1998) temperatures of 16:0/18:4(n-3)cis PC, 16:0/18:5(n-3)cis PC, 16:0/20:5(n-3)cis PC, 18:0/18:4(n-3)cis PC or 18:0/18:5(n-3)cis PC have been published.
It is seen that the temperature T = 303 K of the MD computer simulations is higher than experimental gel to liquid-crystalline phase transition temperatures T c Table 3 Published temperatures T c of the experimental lamellar gel to liquid-crystalline phase transition of mixed-chain phosphatidylcholines containing an sn-1 palmitoyl (16:0) chain and various sn-2 unsaturated fatty acid chains For abbreviations, see footnote in Table 2 Lipid Lipid database Koynova and Caffrey (1998) Davis et al. (1980Davis et al. ( , 1981 and Keough (1986)  Hydrol. meas. Kamp et al. (1975) 280.05 ± 2.9 Lipid database Koynova and Caffrey (1998) of all the bilayers in Tables 2, 3, 4 and 5. In spite of the fact that T c values of several lipids, 16:0/18:4(n-3)cis PC, 16:0/18:5(n-3)cis PC, 16:0/20:5(n-3)cis PC, 18:0/18:4(n-3)cis PC and 16:0/18:5(n-3)cis PC, are unknown, there is good reason to believe from the analysis of the noted Tables that missing values of T c are also less than T = 303 K. This temperature is acceptable also for MC simulations of hydrocarbon chains because the main phase transition temperature of octadecane (alk-18:0) is 301.2 K and that of hexadecane (alk-16:0) is 291.2 K (Dirand et al. 2002), and phase transition temperatures of unsaturated (alkene) chains are substantially lower than those of n-alkanes. A description of the MD simulations technique of PC bilayers was presented in a previous paper (Rabinovich and Lyubartsev 2014). The simulation boxes contained 128 PC molecules of one of 16 studied types per bilayer (64 lipids in each leaflet) and 30 H 2 O molecules per lipid that corresponds to a condition of full hydration (overall 3840 water molecules). The two hydrocarbon tails, the glycerol section and the head group of the lipid molecules were treated in accordance with their known chemical structure, all hydrogen atoms were explicitly included in the computations.
In the starting configuration, the lipids were set parallel to each other, organized in a regular manner in two layers, and water molecules were distributed outside the bilayer. The system was put into a rectangular periodic cell, with the Z axis parallel to the bilayer normal. The size of the box was varied during the simulations under a semianisotropic NPT-ensemble with two degrees of freedom: one in the Z direction and another in the XY direction, so that the box sizes in X and Y direction were equal at each time moment.  Barry et al. (1991) 269.35 ± 1.8 Lipid database Koynova and Caffrey (1998) To calculate the energy of the lipid molecules in the course of MD simulations, the CHARMM27 force field parameter set (Feller and MacKerell, Jr. 2000) with modifications described in a previous paper (Högberg et al. 2008) was used. It was demonstrated in a number of publications that the original CHARMM27 force field has some disagreements with the experiment (Benz et al. 2005), especially in the tensionless isothermal-isobaric (NPT) ensemble simulations (Hyvönen and Kovanen 2005;Sonne et al. 2007). CHARMM27 force field was recommended to apply only with a fixed surface area; otherwise the simulated bilayer tends to form a gel-like state (Feller et al. 2002;Koubi et al. 2003;Jensen et al. 2004;Siu et al. 2008). For bilayers composed of 14:0/14:0 PC lipids, modifications introduced previously (Högberg et al. 2008) provided perfect agreement with experimental data for the area per lipid, as well as with the X-ray structure factor and NMR order parameters. For lipids considered in this work, we used the same partial charges as in Högberg et al. 2008, the lipid head group including esters (these charges were recalculated on the bases of ab-initio computations), while for tails, we used charges adopted from the original CHARMM27 force field (Feller and MacKerell, Jr. 2000), with scaling of 1-4 electrostatic interactions by factor 0.83 (Högberg et al. 2008). All intramolecular bond and angle parameters, as well as Lennard-Jones interactions, were also borrowed from the original CHARMM27 force field (Feller and MacKerell, Jr. 2000). Water molecules were described by the flexible SPC model (Toukan and Rahman 1985). Use of this water model, in connection with the modified CHARMM27 force field, was verified previously (Högberg et al. 2008).
The double time step algorithm (Tuckerman et al. 1992) was used to treat separately fast forces (covalent bonds, angles, torsions, collision Lennard-Jones forces within 5 Å distance) with time step 0.25 fs, and longer range forces with time step 2.5 fs. The long-range electrostatic interactions were treated by the Ewald summation method (Allen and Tildesley 1987). The reciprocal part of the Ewald sum was cut at the condition that the remaining terms do not contribute more than 0.0001 of the total value. The parameter of the Ewald sum was set to = 2.6r cut , and cutoff distance r cut = 14 Å was optimized for computational performance according to Ref. (Fincham 1994). The dispersion correction from the Lennard-Jones interactions outside the cut-off distance was included in the pressure (Allen and Tildesley 1987).
The systems were firstly simulated 1 ns under constant volume and then 1 ns under constant pressure and isotropic cell fluctuations. The obtained configurations were considered as starting points for longer simulations with independent cell fluctuations in Z and XY directions. The time reversible Nose-Hoover constant-temperature-constantpressure algorithm (Martyna et al. 1996) was implemented, with the thermostat and barostat relaxation time set to 30 fs and 1 ps, respectively.
All the PC bilayers were simulated for a total of 100 ns. One of the most fundamental properties of a lipid bilayer and one of the most common ways to determine whether the bilayer system has reached equilibrium is area per lipid A pl . When the area per lipid reaches a stable value, other structural properties (density distributions, NMR order parameters) do not show noticeable trends either. In the present work, area per lipid A pl was calculated as the crosssectional area of the simulation box divided by the number of lipids per monolayer. The time evolution of the A pl values of the each of 16 PC bilayers are shown in Figs. 3a, b.
From the observation of the different time evolution traces and calculation of block averages, we concluded that 20 ns of equilibration time is enough for all the bilayer systems considered. Therefore, the first 20 ns of the simulations (from 100 ns) were disregarded in further analysis. Atomic coordinates were saved each 1 ps in the trajectories.
The average areas A pl calculated as a result of our simulations are presented in Table 6.
It is possible to compare our data with the available published (experimental) data: in Table 7, available experimental average lipid areas for bilayers formed by PC lipids with one fully saturated and one unsaturated chain are collected. It is seen that our results are in good agreement with the experimentally deduced data.

Monte Carlo simulations
According to concepts developed by Flory (1969), the interpretation of the spatial configuration of a linear chain molecule dispersed in a dilute solution can be resolved into two parts: short-range and long-range interactions. (1) The short-range interactions of the chain are determined by the Table 6 Average areas per lipid, A pl , and relative fluctuations of the areas obtained for mixed-chain liquid-crystalline phase unsaturated phosphatidylcholine bilayers by MD simulations of the present work; T = 303 K a Statistical error Δ for 20-100 ns is evaluated from the variance of 10-ns block averages   Separovich and Gawrisch (1996) Gawrisch and Holte (1996) 18:0/18:3(n-3)cis PC 303 0.666 2 H and 31 P NMR Separovich and Gawrisch (1996) Gawrisch and Holte (1996) 18:0/20:4(n-6)cis PC 303 0.706 2 H and 31 P NMR Separovich and Gawrisch (1996) Gawrisch and Holte (1996) 18:0/20:5(n-3)cis PC 303 0.691 2 H and 31 P NMR Separovich and Gawrisch (1996) Gawrisch and Holte (1996) 18:0/22:6(n-3)cis PC 303 0.682 ± 0.004 X-ray diffraction Eldho et al. (2003) 303 0.692 ± 0.009 2 H NMR and X-ray Köenig et al. (1997b) 303 0.716 2 H and 31 P NMR Separovich and Gawrisch (1996) Gawrisch and Holte (1996) geometrical parameters (bond lengths and bond angles), together with the potentials affecting rotation about bonds, including the effects of steric interactions between atoms and groups which are near neighbors in sequence along the chain; in other words, the short-range effects are determined by interactions between groups separated by only a few bonds.
(2) The long-range interactions are dominated by interactions involving pairs of atoms and groups which are remote in the chain sequence, though near to one another in space when involved in mutual interactions; to put it differently, the longrange interactions are determined by interactions between pairs which are separated by many bonds (Flory 1969).
work. From the mathematical point of view, to reach this state, the long-range interactions should be excluded. Thus, the MC simulations were performed here for unperturbed hydrocarbon chains, in which only intramolecular interactions between near neighbours along the chain were included. Let U be the conformational energy of a chain in the unperturbed state, i.e., the short-range interaction energy. The equilibrium properties of all chains were calculated using the classical flexible model (Gö and Scheraga 1976). The average value of the physically observable ⟨H⟩ Θ of an unperturbed molecule in the canonical ensemble for this model is given by where N is the number of carbons of the main chain; k B -the Boltzmann constant; T-temperature; and 1 , 2 , … , N−1 are torsion angles of the main chain. Here, in (1) and below in the MC procedure, we keep bond lengths and bond angles equal to their equilibrium values which correspond to the chosen force field parameters; the torsions for all double bonds arranged cis were also fixed at the equilibrium value.
Assume we generate angles 1 , 2 , … , N−1 with probability density p( 1 , 2 , … , N−1 ). Then, an assessment H of the value ⟨H⟩ Θ by the MC method (Gould et al. 2007;Landau and Binder 2009;Binder and Heermann 2010;Satoh 2011) is Here, is the sample size and is the number of the current conformation. The value H from expression (2) converges to the value ⟨H⟩ Θ from expression (1).
In previous work (Zhurkin and Rabinovich 2015), an important sampling technique was developed for the efficient generation of chain conformations, with continuous variation of all single C-C bond torsions within [0, 2 ] range. The conformations were generated using the probability density The long-range interactions introduce alterations (perturbation) in the chain configuration obtained when only the short-range interactions are considered. It is important to note that the long-range effect depends not only on the actual volume of the chain group (fragment, unit) but also on its interaction with the solvent; it is reasonable, therefore, to discuss the effective covolume. The covolume for the chain group can be enhanced by use of a 'good' solvent for the chain. It may also be diminished by choice of a 'poor' one barely capable of dissolving the chain. Through judicious selection of solvent and temperature, the finite volume of the group can be compensated exactly by the mutual attractions between chain groups when immersed in the poor solvent (Flory 1969). This state was called the 'Theta (Θ) point'; as this takes place, the perturbation of the chain configuration must vanish and the chain become unperturbed (Flory 1969). Further, according to the prediction by Flory, in the bulk amorphous state, perturbation of the chain configuration must vanish. Neutron scattering experiments were later carried out (Dettenmaier 1978;Yoon and Flory 1978), and the results substantiated this prediction.
The 'unperturbed' state of chain molecules was used as the uniform state of the different hydrocarbon chains in this To calculate the conformational energy U of a hydrocarbon chain in the unperturbed state (the short-range interactions energy), a scheme of interdependence of each of three torsions along the chain was taken into account in our work. The energy U was calculated as the sum of energies U m ( , +1 , +2 ) of N f structural units (e.g., N f = N − 3 for a saturated chain with N carbons): where , +1 , +2 are torsion angles, and m is the structural unit type. The units reproduced precisely the structure of various chain fragments.
The energy U is arranged [according to (4)] in such a way that it is possible to calculate energy U m of any m unit at the preliminary step, before MC simulation of the chain. To calculate the energies of all units and chains as a whole, the same force field parameters as in the MD simulations were used (CHARMM27; Feller and MacKerell, Jr. 2000) with modifications performed in the paper (Högberg et al. 2008).
To construct a chain of the above-mentioned structure and then calculate the energy U according to expression (4), several of the 16 structural units presented in Fig. 4 should be properly combined. Three variable torsions in each unit in Fig. 4 are marked by red arrows. It is seen from expression (4) that the structural units in a chain should be connected in such a way that each two consecutive (neighboring) units have two mutually variable torsions.
When calculating U m of each structural unit, the torsion energy, non-bonded energy and electrostatic energy from the force field were considered; bond lengths and bond angles were fixed at equilibrium values. Therefore U according to (4) only approximately; nevertheless, it is quite sufficient for the next step of the algorithm. An exact calculation of the energy U of each conformation is carried out after it has been generated-see below: it will be denoted as U units in the final assessment H of the value ⟨H⟩ Θ .
It is possible to demonstrate here an 'interdependence phenomenon' of three consecutive variable torsions along the chain by the example of structural unit 9 from Fig. 4. Energy U 9 ( 1 , 2 , 3 ) of unit 9 was calculated (tabulated) and presented in Fig. 5 in the form of six two-dimensional energy maps containing equienergy contours, i.e., lines connecting the points of equal energies.
Equienergy contours for a pair of torsions can, of course, be demonstrated at any fixed value of the third torsion angle (not only at 60 • and 120 • ). It is seen that map (a) of Fig. 5 is significantly different from map (b) of this figure; a large difference is observed for maps (c) and (d), maps (e) and (f), i.e., dependence of energy U 9 on torsion angles 1 , 2 , 3 is pronounced (in other words, interdependence exists between the three torsions).
To generate the density (3) for a given chain, a special technique was developed. Let N f = N − 3 for simplicity; then, we can substitute U from (3) by expression (4) and rewrite expression (3) in a different form: The energy U m ( , +1 , +2 ) of each structural unit m was tabulated with a step of 1 • . Then exp[−U m ( , +1 , +2 )∕k B T] values under given T were calculated, and integrals were computed numerically. Then, the configurational space of torsion angles , +1 , +2 of each chain's unit m, where 0 ≤ ≤ 2 , 0 ≤ +1 ≤ 2 , 0 ≤ +2 ≤ 2 (i.e., a 'cube'), was divided numerically into 100 3 = bond length energy and bond angle energy were constant. In the calculations, multipliers of 1/2 and 1/3 were used for some energy items to exclude the possibility of double (or triple) summation of any energy items in the final expression (4). The multiplier 1/2 was used for the energy items which are strictly dependent on two variable torsions, and 1/3 was used for items strictly dependent on one variable torsion. The energy items that are strictly dependent on the three variable torsions of the considered structural unit were accounted without multipliers (i.e., a multiplier was equal to 1). It should be recorded that, at this step, it is possible to calculate the short-range energy Fig. 4 Sixteen structural units; to construct a linear hydrocarbon (n-alkane or n-alkene) chain of the structure, e.g., as in Fig. 2, that is typical for the biomembrane phospholipid chain structure (Fig. 1), and calculate the energy U according to expression (4), several of the presented units should be properly combined. Three variable torsions in each unit are marked by red arrows. The number at the bottom right of the unit is the unit's number 1 3 Fig. 5 Six two-dimensional energy maps of structural unit 9 from Fig. 4. The value 0 • of any torsions (angles 1 , 2 and 3 ) corresponds to the eclipsed conformation. The numbers near equienergetic contours are energies (kJ/mol). The energy of each map is measured from the global energy minimum of structural unit 9 1,000,000 parallelepipeds ('states') in such a way that they all have equal Boltzmann realization probabilities under given T. As a result, boundaries between the 'states' along three directions (torsions) in angle units were calculated; to do that, a special mathematical algorithm was developed. The idea of the algorithm is as follows. At first, values of boundaries along axis were calculated using a recurrent relationship, to divide the 'cube' onto 100 quadratic (2 × 2 ) strata ('layers') having equal realization probabilities and, as a consequence, different widths. Then, each stratum, using a recurrent relationship, was numerically divided along the second axis, +1 , onto 100 'columns' ('rods') in such a way that all 'columns' have equal realization probabilities (and, because of this, different sizes); the values of boundaries between the 'columns' were obtained. Finally, each 'column' was similarly divided along the third axis, +2 , onto 100 equiprobable parallelepipeds ('states') and, hence, all sizes (edge lengths in angle units) of the parallelepipeds are not equal to each other. It is evident that with this method of splitting, the boundaries of parallelepipeds of each structural unit m gather in the areas where energy minima, i.e. the number of states (parallelepipeds) in the energy minima is much more than those around the maxima.
The calculated boundaries between 1,000,000 equiprobable parallelepipeds for each chain's unit m are then used in MC simulations of the different hydrocarbon chains.
is the short-range interactions energy of the generated -th chain conformation, properly calculated using all terms of the force field: all bond length and angle energies, all torsion (dihedral) angle energies, out-ofplane energies, Urey-Bradley terms, and non-bonded and electrostatic interactions energies for such pairs of atoms which are included into the sequence of structural units of the given chain. Since all values of the torsions in the generated conformation are already known, the short-range energy U units is calculated correctly [in contrast to the energy U calculated approximately from expression (4)]. Thus, in the assessment (6), the probability of generation of each chain conformation and probability of its realization are calculated, and hence we obtain H To calculate average characteristics, approximately 10 12 conformations of each chain were generated in the present work.

Results and discussion
Average distances between terminal carbon atoms of the chains (end-to-end distances) considered in the unperturbed state, ⟨h⟩ Θ , and those of PC lipid chains in liquid crystalline bilayers, ⟨h⟩ bil , were calculated. The obtained data are presented in Table 8. Table 8 shows that ⟨h⟩ Θ values are somewhat less than ⟨h⟩ bil . This is because only one hydrocarbon chain terminus in the lipid molecule is free in space; the other one is chemically linked to the head group. Due to interactions of all lipid molecules with their neighbors and water molecules the lipids' head groups are arranged in the vicinity to each other. Therefore, possibilities of rotations around several C-C bonds adjoining the head groups are more restricted than those for the C-C bonds of the opposite end of the chain, so the chain region near the head groups is more stretched, in contrast to the unperturbed chain in which both ends are free.
To compare the results quantitatively, the relative (in percentage) difference between values ⟨h⟩ bil and ⟨h⟩ Θ for each chain was calculated: The values of are also presented in Table 8. The calculations showed that the relative difference between both states increases as the number of carbons and/or number of double bonds in the chain increase; is approximately equal to 8-10% (or ∼9%) for saturated 16:0 and 18:0 chains; to 10-15% for unsaturated chains with N = 18, d = 1-5; to 16-17% for PU chains with N = 20, d = 4 -5 and maximum ∼19% for 22:6(n-3)cis chain.
In this connection, a remark should be made. It is possible to divide interactions in the bilayers into three parts: intramolecular short-range, intramolecular long-range interactions (Flory 1969) of the chains, and intermolecular interactions between the chains and the neighboring chains (and PC head groups). Unperturbed hydrocarbon chain properties are fully defined by the short-range interaction energy (Flory 1969). Therefore, the value of can be considered as an assessment of the influence of the longrange interactions inside the chain and interactions with the neighboring chains and PC head groups of the lipid bilayer on the distance ⟨h⟩ bil , compared with the influence of only short-range interactions inside the chain on the ⟨h⟩ Θ .
While conformations of the unperturbed chains are not the same as those in liquid crystalline bilayers, the relative difference between the end-to-end distances ⟨h⟩ bil and ⟨h⟩ Θ of the considered typical acyls was found to be comparatively moderate; it is approximately equal to 9-19%. Therefore, such properties of the listed FA chains of phospholipids in bilayers, as ⟨h⟩ bil are significantly determined by the short-range interactions in the chains (indeed, it is determined approximately by 81-91%). There is good reason to believe that other geometrical properties of these chains are also determined mainly by their short-range interactions. We mention that the considered saturated and unsaturated chains are widespread, typical constituents of phospholipids, so the obtained result seems to be valuable for biomembranes. In other words, the value seems to be of the same order of magnitude for most of different membrane hydrocarbon chains with methylene-interrupted cis double bonds.
It should be also pointed out, that a relationship between ⟨h⟩ Θ and the number d of double bonds for the unperturbed hydrocarbon chains with constant number N of carbon atoms is the same as ⟨h⟩ bil for acyl chains of PC molecules in bilayers: both values decrease as d increases (Table 8). In other words, the trends in changes of ⟨h⟩ Θ and ⟨h⟩ bil in the order 18:0 → 18:1 → 18:2 → 18:3 → 18:4 → 18:5 are the same; the trends in the order 20:4 → 20:5 are also the same.
This decrease is obviously caused mainly by extension of the chain segment with double bonds, in which energy minima are wide and correspond to various collapsed chain conformations, and by shortening of chain's saturated segments in which the main energy minimum corresponds to the extended chain conformations (Flory 1969).
It is possible to consider also the trends of both discussed values of ⟨h⟩ in the order 18:4 → 20:4 and the trends in the order 18:5 → 20:5, i.e., relationships between ⟨h⟩ Θ , ⟨h⟩ bil and number N of carbon atoms under constant number d of double bonds. Table 8 shows that they are similar to each other, respectively. The discussed trends are shown schematically in Fig. 6. Of course, more rigorous treatment requires consideration also of the value of the parameter j in all trends, i.e., the third structure parameter of the chain in the expression N:d(n − j)cis that means the positions of double bonds along the chain. In particular, j differs for the Table 8 Average end-to-end distances, ⟨h⟩ Θ for unperturbed hydrocarbon chains and ⟨h⟩ bil for the acyl chains in liquid crystalline phosphatidylcholine (PC) bilayers obtained by computer simulations of the present work; T = 303 K a stat. error Δ 1 evaluated from the variance of ∼ 10 12 conformations b stat. error Δ 2 for 20-100 ns evaluated from the variance of 10 ns block averages  To our knowledge, a quantitative assessment of the difference between properties of the two chain states (in particular, the difference in ⟨h⟩ values) was obtained for the first time in the present work, while a qualitative similarity of chain properties in the two states was discussed in the literature before. For instance, it has been concluded (Rabinovich et al. 2003) that the bond-order parameters and orientation distribution characteristics of the chains in the lipid monolayer and bilayer 'liquid' regions, as found in experiments and in MD computer simulation models (Rabinovich et al. 1999a(Rabinovich et al. , b, 2000Rabinovich and Balabaev 2001), are qualitatively similar to the intramolecular order parameters and the intramolecular bond orientation distributions in single unperturbed (Flory 1969) unsaturated hydrocarbon chains previously investigated with MC simulations Ripatti 1999, 2000). Therefore, the behavior of the acyl chains in the liquid region of lipid bilayers (somewhat remote from the membrane-water interface) is dominated by the intramolecular short-range interactions. The long-range interactions of the segments of the lipids in this region of the bilayer and the interactions with the bilayer-water interface may be considered as a disturbance: the intermolecular interactions are largely used to orient the lipid molecules in the direction of the membrane normal.
Thus, from the two above-mentioned facts (about a comparatively moderate quantitative difference in the ⟨h⟩ bil and ⟨h⟩ Θ values of the chains and the similarity of their trends), a common conclusion can be made: to treat and analyse a number of processes in biological membranes (e.g., changes in FA composition caused by the environmental changes such as temperature and pressure), it is possible to use, at least as the first approximation, the relationships between structure and properties obtained for the unperturbed hydrocarbon chains. This seems not unreasonable: biomembranes are known to contain a wide variety of FA chains; the available 'bilayer' relations between their structure and property are incomplete and insufficient for the analysis, whereas the properties of the unperturbed chains and corresponding 'structure-property' relationships have been already studied for tens of variants (see, e.g., Zhurkin and Rabinovich 2015).

Conclusions
The average characteristics of hydrocarbon chains calculated in the unperturbed state (which is fully defined by short-range interaction energies) make it possible to estimate the influence of additional energy components on the state of these chains, if they are under other conditions or located in other systems. For the same temperature and force field parameters, it is acceptable to use any characteristic as a criterion. The average end-to-end distance of the chains was chosen as such a criterion in the computer simulations carried out in this work. The relationships between structure and the average end-to-end distances obtained for the considered unperturbed chains were shown to be qualitatively similar to those of lipid chains in bilayers. Such data for the majority of possible lipid acyls in bilayers (as a rule, it is several tens of chains or more) are yet unknown because MD simulations of various lipid bilayers are very time-consuming. On this basis, it is reasonable to investigate the unperturbed hydrocarbon chains instead. As a first approximation of the desired 'structure-property' relationships for the lipid chains in bilayers, the corresponding relationships for the unperturbed chains can be used.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// Fig. 6 Average end-to-end distances ⟨h⟩ Θ , [nm], for unperturbed hydrocarbon chains obtained by Monte Carlo (MC) simulations, and ⟨h⟩ bil , [nm], for the same acyl chains in liquid crystalline phosphatidylcholine (PC) bilayers obtained by molecular dynamics (MD) simulations (triangles for the marked chains in 18:0/... PC bilayers, circles for the marked chains in 16:0/... PC bilayers). The ranges for saturated acyl chains 16:0 and 18:0 are obtained for eight appropriate mixed-chain PC bilayers. Computer simulations of the present work, T = 303 K. Arrows show qualitatively trends in ⟨h⟩ Θ and ⟨h⟩ bil . To compare the obtained trends, the same names used here both for acyls (ordinate axis) and hydrocarbon chains (abscissa axis) creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.