1. Introduction

Cholesterol (CHOL) is an important component of the membranes of mammalian cells. There is evidence that cholesterol affects membrane structure ^[15], including physiological membranes ^[19]. Ergosterol (ERG) is the fungi and yeast equivalent, or analog, of CHOL. The chemical structures of CHOL and ERG are also strikingly similar, which fact has led us to expect that the structural arrangement of CHOL and ERG in physiological membranes may be quite similar.

Physiological membranes consist of a bilayer composed of two sheets of phospholipid molecules. A specific phospholipid is identified by the polar headgroup and by the length of the acyl chains. The polar headgroup is hydrophilic and the hydrocarbon tails are hydrophobic. The bilayer then has two hydrophilic surfaces and a hydrophobic interior resulting in a resistance of the bilayer to the passage of ions. The electrical resistance of the phospholipid bilayer membrane is of the order of $50-100$ ${\rm{G}}\Omega$.

Our experiments in the Turner Laboratory have dealt almost exclusively with channels formed by the antifungal drug nystatin (NYS) in the presence of ERG or CHOL. These channels are of a putative barrel-like structure consisting of $6-8$ NYS molecules ^[6], which we identify as ERG/NYS and CHOL/NYS channels. When NYS is applied to a living cell the channels formed in the cell membrane disrupt the physiology of the membrane resulting in cell death.

In our apparatus two chambers (CIS and TRANS) containing KCl solutions were separated by a hole of $50-150$ $\mu {\rm{m}}$ diameter, determined by the choice of cup used for the TRANS chamber. We form an artificial planar phospholipid bilayer over this hole isolating the chambers from one another and impose a concentration gradient of 410:150 mmol KCl and an electrical potential between the chambers. Artificial vesicles (AVs) with phospholipid membranes containing ERG (or CHOL), and NYS are porous because of the ERG/NYS or CHOL/NYS channels that form in the vesicle membrane. If introduced into the CIS chamber of our apparatus these porous AVs fuse with the planar bilayer separating the chambers as a result of the KCl concentration gradient between the chambers. The fusion of such an AV then transports a very small patch of membrane containing ERG/NYS or CHOL/NYS channels onto the much larger planar phospholipid bilayer separating the chambers. We can then study the dynamics of the channel activity in the planar bilayer reflected in the variations in planar bilayer current.

Another approach is to form the isolating planar bilayer between the chambers of phospholipids containing a known mol fraction of ERG or CHOL. Adding a small concentration of NYS to the CIS chamber then results in the formation of ERG/NYS or CHOL/NYS channels in the isolating planar bilayer.

Either case these approaches allowed us to study the structures of ERG and CHOL in the AV membranes through the channel dynamics reflected in measurements of the time dependence of planar bilayer currents.

There has been considerable interest in the interaction of the sterols ERG and CHOL with lipids in bilayers and in monolayers. The goal of much of this work has been to understand lateral organization of these sterols in biological membranes ^[9,21,22,23]. There is also considerable evidence regarding the importance of functional sterol microdomains in biological membranes ^[4,5,17,29]. The greater hydrophobicities of the sterols requires them to be shielded from the water by the polar headgroups of the phospholipids ^[4]. The idea that the Lipid headgroups form an umbrella shielding the sterol from the water on either side of the phospholipid bilayer, which has been referred to as the "umbrella effect", was originated by Huang and Feigenson ^[14].

There are three models of the interaction of sterols with lipids that have emerged from these investigations. These were discussed by Ali, et al. ^[1]. (1) A condensed complex model has resulted from thermodynamic studies carried out by McConnell and coworkers. These include lipid monolayers ^[21] and bilayers ^[2]. The condensed complex model has been used to understand the effects of depletion of CHOL in plasma membranes ^[9,19,31]. (2) A superlattice (SL) model was first proposed by Somerharju, et al. to interpret kinks in fluorescence intensity data taken on pyrene-labeled phosphatidylcholine (Pyr-PC) in liposomes ^[26,30]. In the model the acyl chains of the lipids were arranged in a hexagonal lattice. It was proposed that the pyrene-containing acyl chains took positions in the lattice that maximized their separation. In 1992 Tang and Chong studied the spectrum of Pyr-PC in multilamellar vesicles (MLVs) and discovered dips as well as kinks ^[28]. They proposed an extended SL model to account for these observations. The structures predicted by the extended SL model were obtained from Monte Carlo (MC) simulations conducted on sterols interspersed in a hexagonal lattice ^{[12,13,14,20,27]}. (3) The umbrella model in which the polar headgroups of the phospholipids shield the hydrophobic sterols from water by forming something of an umbrella over the sterol ^{[3,8,13,14,20]}.

These models are not completely distinct. Huang and Feigenson introduced the umbrella model in the context of the results they found from the interactions present in their Hamiltonian, which formed the basis of their MC simulations. And Radhakrishnan and McConnel indicate similarities between their thermodynamic condensed complex model and a general superlattice model ^[21]. However, Ali, et al. ^[1] conducted experiments in which large MLVs with membranes containing CHOL were immersed in a bath of cholesterol oxidase (COD) in which they measured the initial reaction rates and related these, through transition state theory (see e.g. ^[11]), to the chemical potential of the CHOL. Their data indicated support for a SL model and an understanding of the dynamics in terms of the umbrella model, but brought the condensed complex model into question for phospholipid bilayers.

Chong measured the density of dehydroergosterol (DHE), a naturally occurring fluorescent sterol ^[4]. He found that a plot of the DHE fluorescence vs. the mol fraction of DHE shows a number of intensity drops (^[4] Figure 1), which he proposed resulted from a difference in the vertical position of the DHE in the membrane in the irregular regions and the regular, or superlattice (SL) regions of the DHE. He reasoned that the depth of the DHE in the SL regions is less than that of the DHE in the irregular regions resulting in a higher dielectric constant for the DHE in the SL regions. The fluorescence of DHE decreases with increasing dielectric constant of the surrounding medium ^[24] Therefore the dips in the DHE fluorescence are associated with an increase in the regularity of the DHE in the phospholipid bilayer. These SL regions were in agreement with the model of Tang and Chong ^[28].

Figure 1. The forms of the superlattices that will fill a hexagonal lattice for ERG mol fractions of (a) 0.40 and (b) 0.50 Each open circle represents a lipid acyl chain and each filled circle an ERG.

DownLoad: Full-Size ImgPowerPoint

Liu, Sugár, and Chong ^[17] studied the dependence of the fluorescence of DHE in mixtures of DHE/cholesterol and DHE/ergosterol in phosphatidylcholine (PC) bilayers on total mol fraction of sterol present in the bilayer. The total sterol mol fractions were those of DHE + ERG or DHE + CHOL, for DHE mol fraction $0.05$ and $0.005$. They observed fluorescence intensity dips whenever the total mol fraction of sterol is close to the critical mole fractions for sterols in SLs ^[4]. The differences in the locations of the peaks were slight for ERG or CHOL, which was not important for the studies conducted by Liu, et al. The presence of SLs was clear. From the results we present here, however, the subtle differences in the location of the fluorescence peaks when CHOL replaced ERG may indicate the rather substantial differences in dynamics we discovered.

In a previous publication we showed that ERG/NYS channels form at the boundaries and not in the interiors of ERG domains ^[12]. This conclusion was based on AV fusion experiments as described above. Our statistical mechanical analysis predicted the linear decay schemes we observed for these bilayer currents and the linear relationship of peak bilayer current to AV diameter we had found in previous experiments. We were also able to subsequently supplement these experiments with MC simulations based on a modification of the Hamiltonian developed by Huang ^[13]. Our MC simulations based on this Hamiltonian resulted in SLs for the specific mol fractions of sterols predicted by Liu, et al. ^[17] and provided a detailed understanding of ERG/NYS channel dynamics when ERG was present in other than the specific mol fraction required for a superlattice.

Complex intermolecular interactions in a system result from changes in the electrochemical potential (see e.g. ^[11]), which is proportional to the entropy of the system. Entropy is a thermodynamic property that cannot be understood in terms of small numbers of molecules, particularly for systems with multibody interactions, for which the original analysis of Ludwig Boltzmann fails (see e.g. ^[11]). The results of MC simulations involving the Metropolis algorithm are driven toward a Boltzmann equilibrium. MC simulations, however, involve only the Hamiltonian in the physical description, and are based on a Markovian description of the system dynamics. Any accurate description of the physics of interaction of sterols at the boundary of the domains rich in ERG should, of course, be based on the thermodynamics of nonequilibrium at the boundary of the domain. In a MC simulation we have only the Hamiltonian to represent any and all interactions. The only way that the Hamiltonian can be modified to model boundary interactions is by adding an energy at the boundary.

For the investigations we describe here the Hamiltonian was modified to include a boundary binding energy (BBE) to model the interaction among ERG molecules in a bilayer based on the results of ^[25]. Sintes and Baumgärtner considered the interaction between two proteins in a lipid bilayer. The proteins were modeled by cylinders, which extended across the bilayer. The sterols in which we are interested do not extend across the total bilayer. However, the study by Sintes and Baumgä rtner reveals aspects of the energy of interaction between membrane molecules distinct from the lipids which are of interest to us. Specifically Sintes and Baumgärtner found two basic regions of lipid mediated attraction between the cylindrical molecules dependent on distance between them. These two regions were separated by a repulsion. In our Hamiltonian we introduced an attractive force between sterol molecules on the boundaries of SL domains. This reflects the results of Sintes and Baumgärtner, although we made no attempt to model the detailed spatial dependence they discovered except to require that our attractive force appears only when the sterol environment satisfies conditions that locate it on the boundary.

With this Hamiltonian, including a BBE, we found, for mol fractions of the sterol lying between those at which a SL fills a hexagonal lattice, that a mosaic of SL domains forms. We chose the hexagonal lattice to model a bilayer of closely packed phospholipid acyl chains. For the case in which the sterol is ERG we shall refer to these mosaics as Ergosterol Superlattice Domains (ESLDs).

Our MC simulations reproduced the ERG density measurements as functions of ERG mol fractions obtained by Chong ^[4] from DHE fluorescence. These simulations also provided a detailed description of the experimentally observed decay scheme for the channel currents. Our AV fusion experiments often revealed decay schemes consisting of multiple linear steps. Our MC simulations revealed that these steps result from differences in the rate at which sterols are released from the SLs corresponding to different mol fractions, which are present in the ESLDs. We can then claim that we have substantial evidence that the structure of ERG in bilayers is a SL or a mosaic of SLs.

However, our expectation that this program of bilayer experiments supplemented by MC simulations would extend to the study of AVs in which CHOL replaced ERG was dashed with our first AV fusion experiments. The dynamics of the nystatin channels supported by cholesterol (CHOL/NYS) simply bore no resemblance whatsoever to the decay scheme characteristic of the ERG/NYS channels. The procedure we had previously used to investigate the structure of ESLDs in bilayers could not then be used to investigate the structure of CHOL in bilayers.

Our primary interest was in whether or not the structure of the CHOL in the phospholipid bilayer took on a form similar to the SL domains we had identified in ERG. It seemed futile to attempt anything other than an experimental approach to the question. It also seemed prudent to try any experimental program first on ERG/NYS channels, since we knew considerable about them, and then extend the program to CHOL/NYS channels.

Based on the work of Coutinho, Prieto, and coworkers ^[6,7] on the interaction of nystatin with membranes and our results that the ERG/NYS channels form on ESLD boundaries, we reasoned that initial formation rates of ERG/NYS channels may reflect the internal ERG structure in the ESLDs. We prepared artificial bilayers of POPE or POPC containing the ERG mol fraction of interest. Then we stirred small quantities of NYS into the CIS chamber to give concentrations of $8.1-32.4\ \rm{nmol}/\rm{ml}$. The experiments produced initial exponential channel growth rates for which the rate constants varied in the same fashion as the ERG density for a SL structure. However, an identical set of experiments failed for CHOL/NYS channels because of the rapid fluctuation of the initial channels.

During our investigations, however, we had noticed the regular appearance of what seemed to be fluctuations of single or small numbers of channels appearing as spikes on our records of the bilayer current. These spikes were of the order of a $p\rm{A}$ and had durations of the order of $0.1\ \rm{s}$. We finally turned to a study of the frequency of these fluctuations.

For our studies we formed bilayers with measured mol fractions of CHOL in POPE or POPC. Then in the CIS chamber we added NYS to produce densities of either $16.2$ or $32.4\ \rm{nmol}/\rm{ml}$. On the first $400\ \rm{s}$ of the bilayer current record we then conducted a Fast Fourier Transform (FFT). To obtain an accurate evaluation of the frequency of the lowest FFT peak we conducted a tenth order orthogonal polynomial fit to the FFT. Our studies encompassed the CHOL mol fraction range of $0.16-0.44$ in order to obtain the frequency of the fluctuation spikes as a function of the mol fraction of CHOL. The fluctuation spikes that had attracted our attention were comparatively large and of the lowest frequency appearing in the current record.

The functional dependence of our FFT data very closely resembled the fluorescence data from Liu, et al. ^[17] for the region covered by the fluorescence data. (their Figure 2D). Our data indicate also a rather well defined, although blunt peak at a mol fraction of $0.333$. This is observed in the DHE fluorescence data of Chong ^[4] and appears in Figure 4 in Liu, et al. for DHE/CHOL/DMPC with a sample incubation time of 28 days.

Figure 2. Results of MC Simulation for ERG mol fraction = 0.43.

DownLoad: Full-Size ImgPowerPoint

Figure 3. Percent of the ERG molecules in each of the ESLDs vs ERG mol fraction in bilayer.

DownLoad: Full-Size ImgPowerPoint

Figure 4. Total percent ERG molecules in ESLDs regardless of the identity of the ESLD plotted against the overall ERG mol fraction of the bilayer.

DownLoad: Full-Size ImgPowerPoint

Liu, et al. conclude that their data support idea that sterols are regularly distributed because of the presence of a bulky and rigid tetracyclic ring ^[4]. This is the basis for the SL structure, as we pointed out above. These data are also consistent with the hexagonal SL, although, because we have proposed no MC model for our CHOL experiments, the form of any model lattice does not enter here. We then conclude that our FFT presentation of the presumed fluctuation spikes indicates the presence of a SL structure in CHOL in artificial planar phospholipid (POPC) bilayers.

2. Materials and Methods

2.1. Artificial vesicle preparation

The artificial vesicles (AVs) were prepared using a mixture ($10\ \rm{mg}/% \rm{ml}$) of synthetic lipids, 1-Palmitoyl-2-Oleoyl-sn-Glycero-3-Phosphoethanolamine (POPE), 1-Palmitoyl-2-Oleoyl-sn-Glycero-3-Phosphocholine (POPC), 1-Palmitoyl-2-Oleoyl-sn-Glycero-3-[Phospho-L-Serine] (POPS), (Avanti Polar-Lipids, Inc., Alabaster, AL) and ERG (Sigma Chemical Co., St. Louis, MO). For the initial AV fusion experiments the lipids were in a 2:1:1 weight ratio and the ERG fixed by the mol fraction desired. In the latter experiments POPC and POPE were used alone in the phospholipid bilayers. Quantities were measured using Hamilton syringes (Hamilton Company, Reno, NV). NYS (Sigma Chemical Co., St. Louis, MO) dissolved in dry methanol ($2.5\ \rm{mg}/\rm{ml}$) was added ($55\ \mu \text{l}$). The solution was dried under nitrogen and $0.25\ \rm{ml}$ of $150\ \rm{mmol}$ KCl ($8\ \rm{mmol}$ HEPES, pH $7.2$) was added. The solution was then vortexed for fifteen minutes and bath sonicated for two minutes then freeze-thaw-sonicate cycle once again. Then it was capped with nitrogen and incubated until use. Finally, the AVs were incubated at room temperature under nitrogen. Before use the AVs were sonicated for $3-4$ seconds.

2.2. Bilayer preparation

Lipid Bilayers were initially prepared using synthetic lipids 1-Palmitoyl-2-Oleoyl-sn-Glycero-3-Phosphoetanolamine (POPE) and 1-Palmitoyl-2-Oleoyl-sn-Glycero-3-Phosphocholine (POPC) (Avanti Polar-Lipids, Inc., Alabaster, AL). In the initial AV fusion experiments a $% 10\ \mu \rm{l}$ suspension of POPE:POPC in a weight ratio of 7:3 was prepared using glass pipettes and dried under nitrogen gas. The dried lipids were resuspended in $50\ \mu \rm{l}$ decane and vortexed for $20$ seconds. The lipid in decane suspension was brushed onto a $50-150$ $\mu {\rm{m}}$ hole to form the planar bilayer using a micropipette. In the later experiments only POPC or POPE were used alone.

2.3. Experiments with NYS stirred in CIS

Bilayer components were prepared using a mixture ($10\ \rm{mg}/\rm{ml}$) of synthetic lipids, 1-Palmitoyl-2-Oleoyl-sn-Glycero-3-Phosphoethanolamine (POPE) or 1-Palmitoyl-2-Oleoyl-sn-Glycero-3-Phosphoch-oline (POPC), (Avanti Polar-Lipids, Inc., Alabaster, AL) and ERG or CHOL (Sigma Chemical Co., St. Louis, MO). Quantities were measured using Hamilton syringes (Hamilton Company, Reno, NV). The solution was dried under nitrogen gas. The dried lipids were resuspended in $50\ \mu \rm{l}$ decane and vortexed for $20$ seconds. The lipid in decane suspension was brushed onto a $50-150$ $\mu {\rm{m}}$ hole to form the planar bilayer using a micropipette. NYS (Sigma Chemical Co., St. Louis, MO) dissolved in dry methanol was added to the CIS. Concentrations of NYS in the CIS ranged from $8.1-32.4\ \rm{nmol}/\rm{ml}$.

2.4. Apparatus

All laboratory experiments were conducted using a planar artificial bilayer. In the experiments we imposed an osmotic gradient of 410:150 $\rm{mmol}$ KCl (CIS:TRANS) and an electrical potential across the bilayer. A combination of the concentration gradient and the electrical potential produced a bilayer current in the picoampere ($p\rm{A}$) range.

The stir bar was driven at $4\ \rm{Hz}$ and the data for the bilayer current were collected at $100\ \rm{Hz}$.

2.5. Monte Carlo Simulations

The general Monte Carlo approach is based on a time independent solution to a Markovian master equation. We used a Metropolis algorithm ^[18] in which all sterol moves resulting in a lower system energy are accepted and moves resulting in higher system energy have an acceptance ratio of $A=\exp \left(-\Delta \varepsilon /kT\right)$, where $% \Delta \varepsilon$ is the energy increase that would result if the move were accepted. The result of a MC step is then a statistical change in the system state, which is biased toward a local Gibbs/Boltzmann thermodynamic equilibrium.

Our MC simulations were all conducted on a hexagonal lattice in which the sites were phospholipid acyl chains or were occupied by a sterol molecule. These lattices had $L=38-98$ sites on a side for the simulations of AVs and lattices of up to $L=240$ sites on a side for the empty lattice in which the AV was placed to model a fusion experiment. We employed periodic boundary conditions on the lattices.

We previously showed that in a MC simulation of the diffusion of sterols from an ESLD a MC step is a representation of the actual time step in the diffusion process ^[12]. This allowed us to interpret our MC simulations of the dispersion of fusion remnants as representations of the corresponding time record. We make no attempt to transform MC steps into real time intervals (see e.g. ^[25]).

We wrote the computer programs for all the MC simulations reported here in Maple code and conducted simulations on a PC with an Athlon 64 3000+ microprocessor and 1.0 GB of RAM. We cross checked our simulations, which resulted in complete SLs, against those of Huang and Feigenson ^[14].

For the initial condition for the simulation of a AV we positioned the number of sterols required to produce the desired sterol mol fraction randomly at sites on the hexagonal lattice. For the initial condition for a MC simulation of a AV fusion we placed the prepared vesicle in the center of the larger lattice.

2.6. Hamiltonian

The Hamiltonian for our MC simulations was based on that presented by Huang ^[13]. To this Hamiltonian we added a term $H_{\text{bound}}$, which represented the binding energy on the boundaries of the ESLDs. Our total Hamiltonian was then

The contributions $H_{\text{sterol}}$ and $H_{\text{chain}}$ to the Hamiltonian appearing in Eq. 1 are described in ^[13]. We chose $H_{\text{bound}}$ to represent an attraction of the sterols in the SL to those on the boundary of the SL. For example for the case of a mosaic made up of SLs with sterol mol fractions of $0.40$ and $% 0.50$ we defined $H_{\text{bound }}$ as

The occupation of the lattice sites is carried by the occupation variables $% L_{\text{si}}$, $L_{\text{x}_{\text{1}}\text{i}}$, and $L_{\text{x}_{\text{2}% }\text{i}}$, which are

and

The occupation variable $L_{\text{x}}$ ensures that if a sterol molecule is on the boundary of both a $0.50$- and a $0.40$- mol fraction SL domain it doesn't experience a binding energy from either SL domain. To do this we defined $L_{\text{x}}$ as

The energies $\Delta E_{\text{B50}}$ and $\Delta E_{\text{B40}}$ are the Boundary Binding Energies (BBE) associated with the SL domains of sterol mol fraction $0.50$ and $0.40$ respectively. The contribution of $H_{\text{ bound}}$ to the Hamiltonian is favorable if the energies $\Delta E_{\text{B50 }}$ and $\Delta E_{\text{B40}}$ are both taken to be $< 0$. The contribution of $H_{\text{bound}}$ then models the cohesive force of Sintes and Baumgä rtner between two cylindrical molecules in a bilayer ^[25]. We note that since the channels, and consequently the NYS molecules, are on the SL domain boundaries, the environment of the sterols on the boundaries of ESLDs differs from that within the SL domains.

The terms $T_{\text{x}}$ are the energy scaling factors for the BBE. These factors must be chosen such that $H_{\text{bound}}$ vanishes if the sterol is not on the boundary of an SL domain. We define a sterol as being on the boundary of a SL domain if there are $3$, $4$ or $5$ sterols in the correct next-nearest-neighbor (NNN) positions. To accomplish this $\left(T_{\text{0 }}, T_{\text{1}}, T_{\text{2}}, T_{\text{6}}\right) =\left(0, 0, 0, 0\right)$. We then chose the remaining values of $T_{\text{x}}$ equal to corresponding scaling values in $H_{\text{sterol}}$ and $H_{\text{chain}}$. Specifically $\left(T_{\text{3}}, T_{\text{4}}, T_{\text{5}}\right) =\left(6, 10, 15\right)$ (see ^[13]). In Eq. 2, however, the index $i$ is a sterol while in $H_{\text{chain}}$ it is an acyl chain.

From the unit cells for the sterol mol fractions $0.50$ and $0.40$ the sterols in the NNN positions determine the identity of the SL. The summations over the indices $x_{\text{1}}$ and $x_{\text{2}}$ in Eq. 2 are over the six possibly occupied NNN sites for SLs with sterol mol fractions $0.50$ and $0.40$ respectively. This prevents the routine from biasing the BBE toward either ESLD.

In our MC simulations we varied all energies over limited ranges. We accepted combinations of energies which produced simulations of fusion events in fundamental agreement with observations. This was true for essentially all energy combinations that produced mosaics of SLs in the model fusion remnant.

There is generally an increase in energy at the boundary of a domain rich in sterols. This results from the necessity of covering the otherwise exposed hydrophobic portions of the phospholipids at the boundary (see e.g., Figures 1 and 8 of ^[16]). As we pointed out in the introduction, however, the Hamiltonian must model an interaction at the boundary that depends on entropy as well as energy. The attraction discovered by Sintes and Baumgärten is representative of this.

In each MC simulation we recorded the number of ERG molecules within and on the boundary of each ESLD. This allowed us to compare our predictions with fluorescence data for ERG density in structured regions and to calculate the electric current passing through the bilayer during the simulation.

3. Results

3.1. ERG/NYS channels in bilayers

In Figure 1 we have the results of MC simulations for ERG mol fractions of (a) $0.40$ and (b) $0.50$ for which Liu, et al. predict SLs (see Table 1 of ^[17]).

And in Figure 2 is an ESLD with ERG mol fraction of $0.43$. Results of MC Simulation for ERG mol fraction = $0.43$.

In Figure 2 each open circle represents a lipid acyl chain and each filled circle an ERG molecule. Except for some vacancies and individual molecules, the phospholipid membrane is filled with a mosaic of ESLDs with SLs corresponding to those in Figure 1.

Each MC simulation began with a specific mol fraction of ERG distributed randomly on the hexagonal bilayer using a (pseudo) random number generator. At each MC step we recorded the number of ERG molecules within and on the boundary of each type of ESLD on the bilayer. This provided us with an indication of when the final state was reached. From these data we could also calculate the relative number of ERG molecules in each ESLD by averaging over the final steps, which we accepted as a theoretical prediction of the relative distribution of the ERG molecules at a particular overall mol fraction.

In Figure 3 we present the results.

From these data we calculate the total percent of the ERG molecules in ESLDs. We present the result in Figure 4.

If we invert our data in Figure 4 and plot them along with Chong's data from Figure 1 of ^[4] we have a comparison of the ERG in ESLDs as predicted by our MC simulation and as measured experimentally by fluorescence of DHE. We present this comparison in Figure 5.

Figure 5. Dehydroergosterol (DHE) in regular regions of ERG and ERG in ESLDs from MC simulation as functions of the ERG Mol Fraction in the membrane.

DownLoad: Full-Size ImgPowerPoint

To model the physics of the decay of these ERG/NYS channels we placed the ESLD mosaic, such as that in Figure 2, into the center of an open phospholipid membrane as we illustrate in Figure 6.

Figure 6. Model of a fused AV in a phospholipid bilayer.

DownLoad: Full-Size ImgPowerPoint

We then conducted MC simulation with this as the initial state obtaining results such as that shown in Figure 7, which is a snapshot after $500$ MC steps.

Figure 7. Dispersion of the ESLDs in Figure 6 after 500 steps into MC simulation. Note the ERG diffusing into the phospholipid bilayer and the changing form of the ESLDs.

DownLoad: Full-Size ImgPowerPoint

In Figure 7 we can see the ERG molecules diffusing into the phospholipid bilayer and the rearrangement of the ESLDs with ERG mol fraction of $0.50$. After $2600$ MC steps we have the remnant shown in Figure 8.

Figure 8. Dispersion of the ESLDs in Figure 6 after 2600 steps into MC simulation. Note that only ESLDs with ERG mol fraction 0.40 remain.

DownLoad: Full-Size ImgPowerPoint

We note in Figure 8 that only ESLDs with ERG mol fraction $0.40$ remain on the phospholipid bilayer.

Since the ERG/NYS channels lie on the perimeters of the ESLDs ^[12], a calculation of the perimeters of each of the ESLDs at each MC step will produce a theoretical plot of the bilayer current resulting from the dispersing mosaic. We present an example of such a plot in Figure 9, comparing it with actual experimental data.

Figure 9. Bilayer current vs time from experiment (left panel) and from MC simulation (right panel). MC simulation contains plots of the total perimeters of the ERG domains with mol fractions 0.50 (red), 0.40 (blue) and total (black). Data and MC simulation with mol fraction 0.43.

DownLoad: Full-Size ImgPowerPoint

3.2. CHOL/NYS channels in bilayers

As we pointed out, our initial experiments on the CHOL/NYS system yielded results from AV fusion experiments that were radically different from those from the ERG/NYS system. In Figure 10 we provide the results of one of our earliest AV fusion experiments.

Figure 10. AV fusion experiment in a phospholipid bilayer. AVs contained CHOL of mol fraction of 0.27.

DownLoad: Full-Size ImgPowerPoint

In this experiment the form of the bilayer currents observed differ dramatically in both magnitude and form from that in Figure 9, left panel, which is characteristic of the bilayer currents we observed at all mol fractions of ERG in AV membranes.

In Figure 11 we present details of the channel dynamics appearing in Figure 10.

Figure 11. AV fusion experiment in a phospholipid bilayer. AVs contained CHOL of mol fraction 0.27. Insert is detail of one of the spikes from an equivalent experiment at a later date.

DownLoad: Full-Size ImgPowerPoint

In the insert in Figure 11, taken from a later experiment, we provide details of a current spike, such as those that we observed whenever the bilayer current differed from zero. These spikes in the bilayer current are regular with frequency increasing as the total bilayer current increases. We also identified channel closings appearing as negative spikes, as are also evident in Figure 11. The duration of the closings may indicate that multiple channels were involved in these events, while the event in the insert in Figure 11 is possibly from the opening of a single channel.

As we described above, these initial experiments with CHOL/NYS channels required us to search for a different approach. In Figure 12 we show the new design in terms of the isolating bilayer separating the two apparatus chambers and the NYS molecules in the CIS chamber.

Figure 12. Idealized drawing of the hole separating the CIS and TRANS chambers in our apparatus. The bilayer painted over the hole contains a known mol fraction of CHOL, represented as small red dots. We added NYS to the CIS chamber with a micropipette and stirred the chamber. The NYS molecules were adsorbed on the bilayer and formed CHOL/NYS channels in the bilayer.

DownLoad: Full-Size ImgPowerPoint

The bilayer painted over the hole contained a known mol fraction of CHOL, represented as small red dots. We added NYS to the CIS chamber with a micropipette and stirred the chamber. The NYS molecules were adsorbed on the bilayer and formed CHOL/NYS channels in the bilayer. It was, of course, impossible to incubate the membrane containing CHOL to assure formation of Cholesterol Superlattice Domains (CSLDs). We could only assume that some CSLDs were present.

Initial experiments with very low concentrations of NYS in the CIS chamber produced results that gave us optimism for a program based on the study of the initial rate of formation of CHOL/NYS in the isolating bilayer of the apparatus. Figure 13 illustrates the result of one such experiment.

Figure 13. Experiment in which NYS was stirred in CIS chamber and bilayer contained CHOL with mol fraction 0.26. Data were gathered to maximum current.

DownLoad: Full-Size ImgPowerPoint

We note, however, the peak in the current at $390\ \rm{s}$ and the fact that the activity at time $<$ $300\ \rm{s}$ cannot be attributed to noise. Nevertheless we elected to pursue a proof of principle investigation by conducting the same experiments with ERG in the bilayer rather than CHOL.

3.3. Initial slope experiments: ERG/NYS

In Figure 14 we present our results for the initial slope experiments plotted with Chong's data for DHE fluorescence as a function of ERG mol fraction.

Figure 14. Dehydroergosterol (DHE) in regular (ESLD) regions and the initial slope of the logarithm of the bilayer current plotted as functions of the ERG mol fraction.

DownLoad: Full-Size ImgPowerPoint

Here we have inverted Chong's DHE fluorescence data so that density appears positively. Except for the fact that the initial peak in our data appears at an ERG mol fraction of $0.19$ rather than $0.20$, our data match those of Chong nicely in the mol fraction regime we studied. We then proposed that a series of experiments on bilayers containing CHOL would result in a variation of the slope of the initial log linear portion that would correspond to Cholesterol Superlattice Domains (CSLDs) that may be present in the bilayer.

3.4. Initial slope experiments: CHOL/NYS

Our experiments on CHOL/phospholipid bilayers, however, provided difficulties. As examples we present two experimental plots in Figure 15.

Figure 15. Two separate experiments conducted on bilayers of POPE and mol fraction 0.25 CHOL with 2 microl of 2.5 mg /ml NYS stirred in CIS.

DownLoad: Full-Size ImgPowerPoint

In Figure 15 we have included a limited number of data points to make the difficulty apparent, while not obscuring the plots and fit lines. Even for the low concentration NYS we added and stirred in the CIS chamber, the channel activity fluctuated considerably. Higher densities of NYS only resulted in larger peaks in bilayer current.

3.5. Fluctuation Spikes in CHOL/NYS Bilayers

The fluctuation spikes we had initially identified in the CHOL/NYS channels (see Figure 10) finally attracted our attention. In Figure 16 we present an example of these spikes as recorded at low current shortly after the NYS was added and stirring began in the CIS chamber.

Figure 16. Spikes in bilayer current. Bilayer is 0.33 mol fraction CHOL in POPC.

DownLoad: Full-Size ImgPowerPoint

Realizing that the spike frequency is generally a function of the average current magnitude (see e.g., Figure 10), we made our measurements on the initial portion of the recordings of the bilayer current data, while the current magnitude was very small.

We focused our interest on the lowest frequency spikes, which seemed to indicate single channel activity. The magnitude of these spikes was also considerably above what may ultimately be noise. We could also visually examine the recorded data and identify more fruitful experiments before turning to the FFT and orthogonal polynomial fit. These final steps were not completely automated and human judgement entered in the final evaluation of the frequency from the polynomial fit. To standardize the analysis, we transformed only the first $400$ seconds after the initial current appeared.

In Figure 17 we have plotted our results for the frequency of the lowest frequency spike appearing in the FFT as a function of the mol fraction of CHOL in the bilayer separating the chambers of our apparatus (black) and for comparison we have also plotted the DHE/CHOL fluorescence data from Liu, et al. ^[17] (their Figure 2D) on the same graph (red).

Figure 17. Lowest Spike Frequency from FFT as a function of CHOL mol Fraction compared to ERG Density from DHE Fluorescence as a function of ERG Mol Fraction.

DownLoad: Full-Size ImgPowerPoint

We inverted the DHE/CHOL data so that an increase in CHOL density is upward. The structure in both plots correspond. To clarify the peaks, which appear to us to correspond to one another we have drawn thin vertical black lines. Our data indicate a rather well defined, although broad peak at a mol fraction of $0.333$. This is observed in the DHE fluorescence data of Chong ^[4] and appears in Figure 4 in Liu, et al. for a sample incubation time of $28$ days.

4. Discussion

The results reported here on the mosaics in the ESLDs were to us insightful. They indicate that in general a phospholipid bilayer will contain ERG (and, presumably, also CHOL) arranged into domains of SLs, even though the total mol fraction of the sterol does not match one for which there will be a SL in the entire bilayer. The fact that our MC simulations resulted in ERG densities that matched those from Chong's DHE fluorescence measurements ^[4] is important. It is a cross check of our MC simulations with experimental data from another laboratory. The fact that our MC simulations provided a consistent picture of the AV experiments, including a description of multiple linear regions in the decay of the bilayer currents, provides substantial support for the emerging theoretical pictures of SL domains and channel dynamics.

We must accept the limitations of MC simulations for any studies of nonequilibrium systems. The Markovian master equation on which they are based has no memory. Without the Metropolis algorithm ^[18], which biases the system towards equilibrium and may be problematic for nonequilibrium systems, time duration of any simulation would be unreasonable. And contact with the physics only through the Hamiltonian is valid for equilibrium thermodynamic systems alone, as Gibbs established ^[10]. We may also raise legitimate questions about how realistic the use of a single hexagonal lattice may be. Physiological lattices are not as simple. And there is finally the question of whether the use of a hexagonal lattice alone might cause an unnatural SL.

We must and have accepted these limitations in favor of the simplicity, flexibility, and transparency of the MC simulations. We believe that the detailed agreement of our MC simulations with our own experimental data, as well as those of Chong and his coworkers, indicates that we may place a certain trust in the results of our MC simulations. These limitations, however, in part led us to reject any attempts to try to discover the reasons for the great differences between the results of the AV fusion experiments for ERG and CHOL.

As we pointed out, all hopes that a similar program, in which studies of CHOL/NYS channels and parallel MC simulations would yield to similar understandings, were dashed when we conducted our first AV experiments involving CHOL instead of ERG. It was tempting to suppose that the channel dynamics we observed resulted simply from a different magnitude in certain items in the Hamiltonian. But it seemed that a more realistic scientific approach was to move instead in a purely experimental direction, at least for the time being.

It seemed prudent to begin a set of experiments in which the isolating bilayer between the chambers carried the CHOL and a minimal concentration of NYS was introduced into the CIS chamber. In this way we would possibly be able to avoid the large groups of channels appearing in the AV experiments. The formation of small numbers of channels may also open some new approaches in the study of the channel dynamics. And it seemed logical to attempt this program first with a known system.

The fact that the initial channel formation rate for ERG/NYS channels in POPE and POPC matched the DHE fluorescence measurements of density ^[4] was encouraging, although not surprising. We had anticipated that the initial channel formation rate should be proportional to the perimeter lengths in the ESLDs. The perimeters, in turn, were functionally dependent on the density of ERG in the ESLDs.

The fact that this approach failed so completely with CHOL seemed to be rooted again in the unknown connectedness of the CHOL/NYS channels and the ESLDs. We had avoided speculation on this when we initially rejected any MC modeling of the dynamics. And in these experiments it was soon apparent that rapid openings and closings of groups of channels made any real identification of an initial slope in the current data impossible.

Our final decision to begin an investigation of the spikes we had been observing on our channel data was neither a stroke of brilliance nor an act of desperation. When we realized these spikes were the result of real channel activity we had spent a period of time looking carefully at them to see if there was any way we could coax some physics out of the form of single spikes. Nothing distinctive had emerged from those initial studies. It was only when we looked carefully at the apparent periodicity, particularly when the total bilayer current was very low, that the idea occurred to us to study the frequency of these spikes.

We began our studies as simply as possible to see what might emerge. We even counted spikes on bilayer current records before we decided that any serious study would require bringing in all of the experimental, mathematical, and technological capabilities we had.

The experiments became very difficult because we had to keep the currents low and record from the point at which the NYS was initially added to the CIS chamber through the appearance of the channels. The analysis was also involved. We turned to the technique of Fast Fourier Transform (FFT) as the only logical approach to identify the periodicities. The differences in magnitudes of the spikes made even a reliable identification of types of spikes by any other means impossible. Fortunately FFT routines are readily available. We used an IMSL (IMSL, Inc., Houston, Texas) routine as the basis for the analysis program we wrote. It was evident visually that the spikes of interest to us were the lowest frequency activity on the bilayer. The data were then in the first peak in the FFT output. To improve our estimates of the location of the center line of the peak we fit the first section of the FFT output with a tenth order orthogonal polynomial. For this we used the orthogonal polynomial fit routine available in the graphical software package that has become our standard, Grapher (Golden Software, Inc., Golden, Colorado).

The results we once again compared with the DHE fluorescence data of Chong and his coworkers. We used the results of DHE/CHOL/DMPC experiments reported by Liu, et al. ^[17]. This comparison we provide in Figure 17. There we have included thin vertical black lines to indicate the corresponding regions in each plot. Even considering that there is a dearth of points in some regions of our data, the maxima in the plots coincide. And the broad peak at a mol fraction of $0.333$ is similar to a peak reported by Liu, et al. after $28$ days of incubation.

5. Conclusion

We began this work as an attempt to answer the question whether or not the ERG and CHOL were arranged in a superlattice pattern in phospholipid membranes. The answer based on our coordinated experimental and MC simulation work seems to be a definitive yes for ERG. And our bilayer current MC simulations lend real credence to our claim in a previous publication that the ERG/NYS channels are on the boundaries of the ESLDs ^[12].

We encountered difficulties, however, in the behavior of CHOL/NYS channels. Our final experiments conducted on the fluctuation spikes appearing on top of the bilayer currents from CHOL/NYS channels yielded data with which we found we could make a statement. We found that these data agree well with the DHE fluorescence data of Liu, et al. ^[17]. That fact we claim reflects an underlying Cholesterol Superlattice Domain (CSLD) structure similar to that we found in ERG.

The theoretical basis for the behavior of these spikes may be as simple as the fact that fluctuations increase in ordered and dense regions. There may also be more complex reasons for this behavior rooted in the cause of the rapid changes in groups of channels observed in the CHOL AV experiments. Here we simply cannot and will not speculate. These statistical fluctuations, we expect, are interesting in their own right.

Acknowledgements

We acknowledge the support of the Turner Laboratory for this research as part of the instructional program in physics at Goshen College. The Turner Laboratory provided the facilities and equipment for this work. We also thank the Goshen College Maple Scholars Program for summer support of the majority of the undergraduate students involved in this research.

Conflict of Interest

The author declares no conflict of interest in this paper.