Article Information

• PDF (730 kb)
• Full Text with Large Figures
• Cited by in Scopus (1)

PubMed

• Articles by G. Richardson
• Articles by P. O’Shea

Related Articles

• Parameter Inference for Biochemical Systems that Undergo a H...

  Parameter Inference for Biochemical Systems that Undergo a Hopf Bifurcation Biophysical Journal, Volume 95, Issue 2, 15 July 2008, Pages 540-549 Paul D.W. Kirk, Tina Toni and Michael P.H. Stumpf Abstract The increasingly widespread use of parametric mathematical models to describe biological systems means that the ability to infer model parameters is of great importance. In this study, we consider parameter inferability in nonlinear ordinary differential equation models that undergo a bifurcation, focusing on a simple but generic biochemical reaction model. We systematically investigate the shape of the likelihood function for the model's parameters, analyzing the changes that occur as the model undergoes a Hopf bifurcation. We demonstrate that there exists an intrinsic link between inference and the parameters’ impact on the modeled system's dynamical stability, which we hope will motivate further research in this area. Abstract | Full Text | PDF (594 kb)

• Origin of Bistability in the lac Operon

  Origin of Bistability in the lac Operon Biophysical Journal, Volume 92, Issue 11, 1 June 2007, Pages 3830-3842 M. Santillán, M.C. Mackey and E.S. Zeron Abstract Multistability is an emergent dynamic property that has been invoked to explain multiple coexisting biological states. In this work, we investigate the origin of bistability in the operon. To do this, we develop a mathematical model for the regulatory pathway in this system and compare the model predictions with other experimental results in which a nonmetabolizable inducer was employed. We investigate the effect of lactose metabolism using this model, and show that it greatly modifies the bistable region in the external lactose () versus external glucose () parameter space. The model also predicts that lactose metabolism can cause bistability to disappear for very low . We have also carried out stochastic numerical simulations of the model for several values of and . Our results indicate that bistability can help guarantee that consumes glucose and lactose in the most efficient possible way. Namely, the operon is induced only when there is almost no glucose in the growing medium, but if is high, the operon induction level increases abruptly when the levels of glucose in the environment decrease to very low values. We demonstrate that this behavior could not be obtained without bistability if the stability of the induced and uninduced states is to be preserved. Finally, we point out that the present methods and results may be useful to study the emergence of multistability in biological systems other than the operon. Abstract | Full Text | PDF (441 kb)

• Spontaneous Creation of Macroscopic Flow and Metachronal Wav...

  Spontaneous Creation of Macroscopic Flow and Metachronal Waves in an Array of Cilia Biophysical Journal, Volume 92, Issue 6, 15 March 2007, Pages 1900-1917 Boris Guirao and Jean-François Joanny Abstract Cells carrying cilia on their surface show many striking features: alignment of cilia in an array, two-phase asymmetric beating for each cilium, and existence of metachronal coordination with a constant phase difference between two adjacent cilia. We give simple theoretical arguments based on hydrodynamic coupling and an internal mechanism of the cilium derived from the behavior of a collection of molecular motors to account qualitatively for these cooperative features. Hydrodynamic interactions can lead to the alignment of an array of cilia. We study the effect of a transverse external flow and obtain a two-phase asymmetrical beating, faster along the flow and slower against the flow, proceeding around an average curved position. We show that an aligned array of cilia is able to spontaneously break the left-right symmetry and to create a global average flow. Metachronal coordination arises as a consequence of the internal mechanism of the cilia and their hydrodynamic couplings, with a wavelength comparable to that found in experiments. It allows the cilia to start beating at a lower adenosine-triphosphate threshold and at a higher frequency than for a single cilium. It also leads to a rather stationary flow, which might be its major advantage. Abstract | Full Text | PDF (5944 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 12, 4145-4156, 15 June 2007

doi:10.1529/biophysj.106.090233

Biophysical Theory and Modeling

Toward a Mathematical Model of the Assembly and Disassembly of Membrane Microdomains: Comparison with Experimental Models

G. Richardson^*, ,  , L.J. Cummings^*, H.J. Harris^†, 1 and P. O’Shea^†

^* School of Mathematical Sciences and University of Nottingham, Nottingham, United Kingdom
^† Cell Biophysics Group, School of Biology, University of Nottingham, Nottingham, United Kingdom

Address reprint requests to G. Richardson, School of Mathematical Sciences, University of Nottingham, Nottingham, NG7 2RD UK.

¹ H. J. Harris’ present address is Division of Infection and Immunity, Institute of Biomedical Research, University of Birmingham, Birmingham, B15 2TT UK.

The term “raft” was first coined by Simons and Ikonen 1 in 1997 to describe lipid/protein microdomain structures that are observed within eukaryotic plasma membranes. Since then research interest in this area of cell biology has grown exponentially. The nature of these structures (and in fact even their experimental detection) has been debated vigorously (see, e.g., Lagerholm et al. 2), as has their role in controlling processes such as membrane trafficking, signal transduction within the cell, endo/exocytosis (including virus entry into cells); and in many biochemical reactions occurring within the cell membrane.

Rafts may be characterized in various ways, but perhaps they are mainly defined by their property of detergent-insolubility in, for example, Triton X-100 or Brij98 2. They may also be thought of as “viscous patches” on the cell membrane since, when the plasma membrane is viewed as a 2D fluid sheet, the raft phase is more viscous than the non-raft phase. Another characterization of rafts is, using the terminology of phase diagrams for multi-component membrane systems, as so-called “liquid-ordered” cholesterol- or sphingomyelin-enriched domains within cell membranes.

The phenomenology of the various kinds of biological processes in which rafts are involved has been discussed at length, and some consensus has emerged on their possible biological roles. As (essentially) phase-separated regions within the 2D membrane, it is thought that rafts can recruit certain reactants and prevent their interaction with other reactants in the “fluid-mosaic” membrane 3, or, conversely, bring desired reactants (particularly proteins) into close proximity, thus promoting certain reactions 4,5,6. In each case proteinaceous receptors and smaller ligands or protein-protein interactions may underlie the biological response, but both cases are conceptually analogous.

Thus rafts may play many very important roles in cell biology, although it must be conceded that some questions remain 2,7. Some of this uncertainty may well reside in the enormous complexity of real cell membranes, and the very small size of rafts in vivo (100–200nm in diameter). Similarly, the basic principles that control the formation and function of rafts within cells remain poorly understood 8,9. To shed light on these fundamental issues we have embarked on a series of studies of simple in vitro model systems, generated in the laboratory in which larger “rafts” can be observed, and the key processes governing their dynamics identified. The hypotheses made as a result of the experimental observations can then be tested by formulating and solving a mathematical model of the experimental system.

The model system we study is the simplest possible that gives rise to interesting and informative “raft” dynamics. (Although we utilize the term “raft” relatively indiscriminately in this article, we are aware that this terminology has certain meanings outside of the simple definition we take here, namely, that it represents a distinct microdomain within the membranes.) We work with planar lipid bilayers composed of precise ratios of phosphatidylcholine (PC; egg lecithin in our experiments) and cholesterol. This model of a cell membrane can then be successively elaborated until genuine “raft” behavior is identified. PC-cholesterol mixed composition bilayers have been studied in detail by many authors. Fig. 1 is adapted from Fig. 11 A of Silvius et al. 10 and shows the now well-known phase diagram for ternary PC-cholesterol mixtures, in which the PC contains variable proportions of dipalmitoyl phosphatidylcholine (DPPC) and a bromo-substituted derivative, 12BrPC (indicated along the bottom horizontal axis of the triangle). Although not strictly representative of the mixture of PC lipids in egg lecithin, the phase diagram provides valuable information on what might be happening in terms of phase in the PC-cholesterol membranes used in this study. The ratio of saturated and unsaturated PC lipids in egg lecithin is ∼45% saturated, 55% unsaturated. If we allow DPPC to represent the saturated PC lipid in the ternary phase diagram, then an indication of what the membrane phases are at certain cholesterol concentrations is given by the dashed line in Fig. 1.

Display large version of this figure

Figure 1

Ternary phase diagram for the 12BrPC/DPPC/cholesterol system at 25°C. The various phases found in this ternary system are liquid disordered, l_d, liquid ordered, l_o, and gel phase, P′_β. Further details may be found in Silvius et al. 10.

We observe cholesterol-rich microdomains (the rafts) forming spontaneously within the bilayer. In terms of the phase diagram of Fig. 1, we identify the raft phase with the liquid-ordered phase. The bilayer was left to equilibriate for 24h, until a final (dynamic) equilibrium distribution of raft sizes (surface areas) was obtained. This raft size distribution was recorded, using fluorescence microscopy (in which fluorescent labeling molecules associate preferentially with the raft phase). We also carried out atomic force microscopy (AFM) studies, which detect the rafts by virtue of the different molecular chain lengths associated with the different lipid molecules. These results are described in separate publications (H. J. Harris, S. M. Rigby-Singleton, M. C. Davies, S. Allen, and P. O’Shea, unpublished, and 30), but we have included a representative image in Fig. 2. We note here that related imaging studies have been carried out by other authors 12,13.

Display large version of this figure

Figure 2

(Upper) Fluorescent images of FPE-labeled PC_100% (left), and PC_67%/cholesterol_33% (right). The images were obtained at room temperature with the bilayers submerged in Tris buffer (10mM Tris, pH 7.4). (Lower) AFM image of the bilayer, illustrating cholesterol-enriched microdomains. The grape corresponds to the membrane topography at that location. Note the topology change (highlighted), where two microdomains are sintering together (or possibly splitting apart).

A mathematical model, based on the Smoluchowski theory of coagulation and fragmentation 14, was formulated to describe the interactions (binding and unbinding) of the cholesterol molecules. The expressions for the binding and dissociation rates of cholesterol molecules were derived using thermodynamic principles. The mathematical model was solved to predict the dynamics of raft formation and disassociation in the simple model system, and the results of numerical simulations compared with the experimental data. A number of key features that the model yields appear to be similar to what is observed experimentally, and thus may also feature in cellular membranes.

We also formulated the problem in terms of the Gibbs free energy, by making the usual assumptions about the entropy and enthalpy of a mixture of reacting molecules. However, the results of this calculation (the Gibbs free energy minimization) gave raft size distributions that were very heavily weighted in favor of very small cluster sizes, and did not agree at all with the experimental observations. We believe this is because the standard theory does not take account of the different dynamics of “molecules” (microdomains, in our system) of very disparate sizes.

One of the virtues of our approach (both experimental and theoretical) is that it lends itself to systematic further elaborations, such as the incorporation of additional lipid types and the inclusion of membrane proteins. These could be modeled in a similar manner to the cholesterol modeling, with averaged properties, to embody more realistic models of biological membranes. The obvious advantage of such a systematic approach is that it allows the effect of each new additional “complication” to be elucidated, as we build up the complexity from the very clean, simple system considered here.

The following sections detail the experimental methods and results, and also the assumptions underlying the mathematical model we use to describe the experiments. The results of the simulations versus experiments are then discussed, as well as the implications of our results for more complicated systems.

PC was purchased from Lipid Products (Kent, UK). Polycarbonate filters (100-nm pore size) were purchased from Nucleopore Filtration Products (Pleasanton, CA). A pressure extruder bomb for model membrane preparation was obtained from Lipex BM Inc. (Vancouver, Canada). Fluorescein phosphatidylethanolamine (FPE) was synthesized from 1,2-dipalmyitoyl phosphatidylethanolamine and fluorescein purchased from Sigma Chemicals (Pool, Dorset, UK), as described by Wall et al. 15. Calcium chloride, cholesterol, dimethyl sulphoxide, DPPC, dioleoyl phosphatidylcholine, magnesium chloride, and Tris (hydroxymethyl) methylamine, were purchased from Sigma Chemicals.

Phospholipids (13mM), dissolved in chloroform and methanol, were mixed in a round-bottomed flask and dried under a stream of nitrogen until a thin lipid film was formed. The dried lipid film was rehydrated with Tris buffer (10mM Tris, pH 7.4), which was quickly frozen in liquid nitrogen and thawed five times. Finally, the suspension was extruded 10 times through a 25-mm-diameter polycarbonate filter (100-nm-diameter pores). For all lipids to be in a fluid phase, the lipid suspensions were heated to 45°C before and during extrusion, resulting in a monodisperse, unilamellar suspension of 100-nm-diameter phospholipid vesicles (PLVs) 16.

PLVs were labeled exclusively in the outer lamella with FPE, as described in Cladera and O’Shea 17. The PLVs (13mM lipid) were incubated in the dark with FPE (30μM in ethanol) at 37°C for more than 1h. The unincorporated FPE was removed by size exclusion chromatography using a Sephadex PD10 column, equilibrated with Tris buffer (10mM Tris, pH 7.4). This procedure leads to the incorporation of 30–50% of the externally added FPE to the PLVs.

Clean 22-mm Ø glass cover slips were treated with magnesium chloride (10mM) at room temperature for 1h. Fluorescently labeled liposomes (1.3mM), heated to 45°C, were then added to the coverslip and left in the dark at room temperature for more than 5h. The fusion of the liposomes, resulting in the unilamellar bilayer covering the glass slide, was pioneered by Watts et al. 18 and results in a bilayer forming on a layer of water on the solid support (the glass coverslip).

The single and two-photon fluorescence microscopy was carried out with a Leica (Wetzlar, Germany) SP2 MP, which utilizes a laser scanning confocal microscope (LSCM). The pinhole required for single-photon microscopy was set at one Airy unit. Alternatively, a Leica DMIRBE inverted fluorescence microscope equipped with a laser source and a mercury lamp/monochomator assembly with a LaVision (Goettingen, Germany) cooled charge-coupled device (CCD) camera was used. Both imaging systems were equipped with a thermostatted stage for temperature control. The optical slicing capability of the LSCM system was not necessary; we simply suggest that our experimental protocol can also be implemented with such a system if available instead of high-sensitivity (e.g., CCD-based) fluorescence microscopy. To prevent undersampling, and to increase the number of intensity levels of images, 512×512-pixel (12-bit) images were taken. Both single- and two-photon microscopy imaging were carried out with objectives with a magnification of 63 and a numerical aperture of 1.32, in illuminating wavelength of 488 or 490nm. However, no differences were found between the image collection regimes. The laser power was kept at a minimum, and the offset and the photomultiplier voltage or CCD output were optimized at the beginning of each experiment and kept constant throughout to determine any contribution from photobleaching and for comparisons between experiments.

The images obtained from LSCM experiments (see Fig. 2 for an example) had a total of 1024² pixels and 0–4095 intensity levels (12-bit images). The images underwent particle analysis using the Scion Image software (Scioncorp, National Institutes of Health, Bethesda, MD) yielding the areas of the more fluorescent patches present in the membranes. The background “noise” was removed to avoid recording spurious microdomains; this process is illustrated schematically in Fig. 3. Patch areas were recorded for any connected set of at least five fluorescing pixels. The raw data was recorded as a column of figures: raft number versus its area.

Display large version of this figure

Figure 3

Method of image analysis. The intensity of the individual pixels in the fluorescent images (a) were plotted (c), and the background noise was removed, resulting in image b. A histogram of the various patch sizes was then plotted.

Mixed composition bilayers were left to equilibrate for 24h at 20°C before fluorescence scanning images were obtained. Three different controlled compositions were used: 33% cholesterol/67% PC; 20% cholesterol/80% PC; and 5% cholesterol/95% PC. A typical fluorescence scanning picture is shown in Fig. 2, for a 33% cholesterol bilayer (this figure also shows one of our AFM images for the PC-cholesterol system, which gives a clearer idea of the raft structure).

The fluorescence images were analyzed as described above, and the data for microdomain areas thus obtained was represented in histogram form, with the bin size chosen as desired for comparison with the theoretical model. Histograms for each of the three data sets are presented later. The data is represented in the dimensionless form used for the theoretical model, which is explained in the following section.

We model the experimental system mathematically, applying the Smoluchowski theory of coagulation and fragmentation 14 to an idealized system in which a large number of cholesterol molecules in 2D clusters of differing sizes (the rafts, or patches) are diffusing around in an otherwise inert 2D fluid (the PC bilayer). Although this approach has been widely used to describe coagulation and fragmentation in a variety of systems 19,20,21,22,23,24,26,27, previous applications to the problem of lipid-raft formation, even in a very simple model system of the kind considered here, are almost nonexistent. The only previous such model we are aware of is by Turner et al. 28, who use such an approach to study the effect of membrane recycling on microdomain size distributions within living cells.

Consider first coagulation (or binding) events between groups or clusters of cholesterol molecules. We assume that these are rate-limited (slow) rather than diffusion-limited (fast), i.e., that relatively few of the collisions that occur between clusters of sizes i and j (clusters containing i and j cholesterol molecules, respectively) result in a coagulation event (a reaction). The concentration of clusters of any given size is thus fairly uniform throughout the bilayer. The reaction rate is proportional to the number of collisions, and hence to σ(v_i+v_j)[c_i][c_j], where v_i and v_j are the mean cluster velocities, [c_i] and [c_j] are the concentrations of i and j clusters, and the collision cross-section σ is taken as the sum of the cluster radii (proportional to (i^1/2+j^1/2)). A schematic is shown in Fig. 4. Postulating that kinetic energies are proportional to thermal energies for the clusters gives v_i ∝ (2kT/(mi))^1/2, v_j∝(2kT/(mj))^1/2, where k is Boltzmann’s constant, T is the system temperature, and m is the mass of an individual cholesterol molecule. Including an Arrhenius activation energy term to account for the fact that only a certain proportion of collisions result in coagulation gives the coagulation rate coefficients

Display large version of this figure

Figure 4

Schematic illustrating collision of clusters of size i and j before coalescing. Cluster radii scale as i^1/2 and j^1/2.

For fragmentation of an (i+j) cluster into an i and a j cluster, we postulate that the rate should be proportional to the circumference of the fragmenting cluster, again with an Arrhenius term. Thus, we take fragmentation rate coefficients

In terms of dimensionless variables c_i*, t* (time) and the parameter M, defined by

(1)

(2)

(3)

The experimental observations have millions of molecules per raft, which corresponds to the parameter M being very large in the above model. Solving this discrete system numerically is prohibitively expensive in this regime, but if we introduce the scalings

(4)

(5)

(6)

(7)

Note that the parameter ϵ characterizes the width of the dimensionless fragmentation kernel b(x,y), in the sense that for there is almost no fragmentation. Thus, in our model, the size of fragments that can break off is limited. Translating back into the dimensional variables, we find that almost no fragments of size

The experimental data gives the areas α of all rafts (above the microscopic resolution size α_min), in a given total area a_T of the bilayer, when the system is judged to be in a steady state. A histogram is made of the numbers of rafts n_m with areas in the range α_m to α_m+1=α_m+Δ, where Δ is some small area increment (the bin size). The mathematical model predicts that the number of rafts n_m in the mth bin is

Thus, we should compare experimental plots of versus δα_m/A (the histogrammed data, scaled appropriately), with theoretical plots from the model (Eq. (4)) of c(x) versus x, taken at large enough times that the model solution has reached its steady state.

The dimensional constants Δ, a_T and are known. The average area occupied by a cholesterol molecule in a raft, A, can then (in principle) be evaluated by considering the total area of all rafts measured experimentally, Ω(α_min) (where Ω(α) denotes the combined total area of all rafts with areas greater than α). We should have

However, in addition to the experimental limitation of the fluorescence microscopy resolution size α_min, there is a numerical restriction due to the fact that the solution of Eq. (4) can only be evaluated on a finite domain 0≤x≤x_max. When comparing the numerical solution to the experimental data it is crucial that the total raft area between x=δα_min/A and x=x_max (or equivalently α=α_min and α=x_maxA/δ) is the same in both numerical solution and experiment. This amounts to the condition

(8)

(9)

Given δ, we can determine K* in terms of the experimental parameters α_min, Ω(α_min), a_T, and With x_max, α_min, and K* known, we can solve Eq. (8) to find x₀, whence A can be determined from Eq. (9). We are then left with just two parameters we can adjust to fit the data: δ=M^−2/3 and ϵ, which characterizes the width of the dimensionless fragmentation kernel b(x,y); once this fitting is done, the resulting value for A provides a further check on the reasonableness of the model.

Best fits of the theoretical model to three sets of experimental data are shown in Figure 5 and Figure 6 and Figure 7, corresponding to PC-cholesterol bilayers containing, respectively, 33%, 20%, and 5% cholesterol. In each case plots are shown of the experimental histogrammed data, scaled as explained above, compared directly with best-fit plots of c(x), both for rafts of small size (where most of the cholesterol mass in the system resides) and for larger raft sizes. We also plot the cumulative cholesterol mass for each case, that is, the proportion of the cholesterol in the system that resides in rafts of (dimensionless) size ≤x, versus x. From these cumulative mass plots it is clear that the model predicts that most cholesterol resides in rafts of dimensions that are less than diffraction-limited spots (i.e., rafts that are too small to measure). An alternative experimental approach that would have enabled us to measure the size of all rafts present would have been AFM; and indeed we did carry out a small number of AFM studies on model membranes 11. However, the disadvantage of this procedure is that the size of membrane that can be scanned is very limited. Thus, it is difficult (with current AFM technology) to obtain a statistically significant set of raft areas with which to compare the mathematical model.

Display large version of this figure

Figure 5

Comparison of theoretical predictions with experimental measurements for the 33% cholesterol membranes. Experimental measurements were made over a total membrane area of 40μm×40μm. Best-fit parameter values for the theoretical model are δ=1/8000, ϵ=0.2375. (a and b) Experimental data shown as histograms of number of rafts versus number of cholesterol molecules per raft (both scaled as discussed in text), against the theoretical curve. Different bin sizes are used in the two histograms. (a) The focus is on smaller rafts, where most of the cholesterol mass in the system resides. (b) The whole range of raft sizes is represented. (c) Cumulative cholesterol mass in the system as a function of raft size, with the solid line representing the theoretical prediction and the dots the experimental data points.

Display large version of this figure

Figure 6

Same description as for Fig. 5, but for the 20% cholesterol membranes. Experimental measurements were made over a total membrane area of 240 μm × 240 μmm. Best-fit parameter values for the theoretical model are δ = 1/60000, ε = 0.2375.

Display large version of this figure

Figure 7

Same description as for Fig. 5, but for the 5% cholesterol membranes. Experimental measurements were made over a total membrane area of 40 μm × 40 μmm. Best-fit parameter values for the theoretical model are δ = 1/300, ε = 0.23375.

For the 33% cholesterol membrane (Fig. 5), the experimental measurements were made on a region of bilayer of size 40μm×40μm, and the best fit to the experimental data was found to be δ=1/2000, ϵ=0.2375, giving a value A_33%=38Ų for the area occupied by a cholesterol molecule within the raft. In this case the theory predicts that the experimental measurements captured ∼25% of the total raft area present. This allows us to estimate the total proportion of membrane area occupied by rafts, P_33%, from the experimental data as

For the 20% cholesterol membrane, measurements were made over a larger region of membrane, 240μm×240μm, and the best-fit parameter values were δ=1/15000, ϵ=0.2375, giving A_20%=80Ų for the area/cholesterol molecule. Here we estimate that the experimental measurements captured ∼15% of the total raft area present, from which we find the total proportion of membrane occupied by rafts, P_20%, as

hence, the actual cholesterol molecule area is estimated as

We note that the area/cholesterol molecule within the raft, A, increases as the cholesterol concentration in the bilayer mixture decreases. This is in line with the above results that show the rafts becoming much more dilute in cholesterol as the bilayer as a whole becomes more dilute. The actual molecular areas show no definite trend, however. Our obtained values for may be compared with data from various sources summarized by Edholm and Nagle 29. This article reports values of cholesterol molecule areas in mixed DPPC/cholesterol bilayers with varying cholesterol concentrations. Using the above estimates for the cholesterol concentration within the rafts (93.5%, 31.5%, and 6.5%, respectively, for the three data sets), we take the figures that Edholm and Nagle report for bilayers that are very dense in cholesterol (27Ų, to be compared with our value ), for bilayers that are 33% cholesterol (40.5Ų, to be compared with our value ), and for bilayers that are 6% cholesterol (51.8Ų, to be compared with our value ), these being the results in Edholm and Nagle closest to those required.

Thus, the molecular area for the 33% cholesterol experiment is remarkably close to the desired result, whereas the other values are certainly qualitatively acceptable—remarkably so, given the simplicity of our model and the minimal fitting.

We observe also that the model is able to capture very closely the cumulative mass as a function of raft size in each case. This is particularly true of the results for the 33% cholesterol bilayers, where the experimental data is sufficiently good to demonstrate the sharp decay of the cumulative mass toward zero as the raft size decreases to zero (Figure 5c). The sharp decay of the data toward zero is almost perfectly captured by the numerical results.

An interesting peculiarity of our results is that, though the fitted value for the parameter δ varies widely between the three data sets, the best-fit value of ϵ is almost the same in all cases, ϵ=0.2375, ϵ=0.2375, and ϵ=0.23375, for the 33%, 20%, and 5% cholesterol experiments, respectively. Further numerical calculations of the solution to Eq. (4) show that, for ϵ <∼ 0.23, the total mass within the computational domain decreases with time as mass is lost to large values of x. It is not possible to tell from the computations whether this is as a result of a genuine bifurcation as ϵ passes some threshold value, in which mass is lost to infinity, or whether it is a consequence of the finite size of the computational domain. However, we show explicitly in the Appendix that a closely related (though simpler) toy model has a bifurcation in its steady-state solution. By comparing the small ϵ limit of the full model (Eq. (4)) (in the steady case) with the toy model, we are able to argue that we believe such a bifurcation is highly likely in our system as ϵ decreases, and, moreover, that all the experiments lie rather close to this bifurcation point.

Although we have good reasons for formulating a dynamic mathematical model (to compare theoretical results to future dynamic experimental measurements), it is of interest to consider whether a steady-state model based on thermodynamic considerations can be formulated to describe the experimental results presented here. To do this, we need to write down a Gibbs free-energy density for the system.

We postulate that the dimensionless internal energy of a cluster of size j is given by

(10)

(11)

We have to minimize the Gibbs free energy, subject to the constraint that the total amount of cholesterol in the system is conserved. Nondimensionalizing as earlier, scaling concentrations with γ₀, and h with kT, and G with and working in the dimensionless variables henceforth, we must minimize

(12)

(13)

(14)

(15)

It is immediately obvious that the above approximate solution (Eq. (12)) for c_i can never be made to fit the experimental data, even approximately, since the decay in c_i with i is always rapid unless the concentration of DPPC is absolutely tiny relative to that of cholesterol (corresponding to ). In this context, we note that there are various other scenarios in which it is inappropriate to apply equilibrium thermodynamics to determine equilibrium constants (e.g., diffusion-limited reactions).

Finally, we note that exactly the same solution for the c_i can be derived by the alternative approach of using reaction coordinates ξ_ij to parameterize the equilibrium in the reaction

We have presented three experimental data sets detailing microdomain (“raft”) formation in PC-cholesterol bilayers, and have derived a simple mathematical model based on plausible thermodynamic arguments that describes the experimental findings well for appropriate choices of two fitting parameters. A further check on the reasonableness of the theoretical results is provided by the predictions of the areas of cholesterol molecules within the rafts, which are qualitatively correct. Moreover, our results exhibit the right trend whereby rafts in more dilute cholesterol bilayers are themselves more dilute in cholesterol.

The mathematical model makes several simplifying assumptions, which may not be fully justified and which could be addressed by a more sophisticated model. One obvious point is that the phosphatidylcholine is assumed not to play a dominant role as far as the model is concerned. Although the PC may not play an important role in the simple experiments described here, we will in future consider more complex membrane systems that are more representative of real cell membranes. This will require development of theoretical models that can describe the simultaneous interaction of several bilayer components.

Our results may be discussed in light of the classical theory of phase transitions, referred to briefly in the Introduction. The phase diagram for PC-cholesterol mixed bilayers given in Silvius et al. 10 is reproduced in Fig. 1. The dashed line indicates where in the phase diagram we believe our experiments to lie, according to the fraction of DPPC in the PC component of the mixture. From this diagram we infer that, for the experiments with 33% cholesterol in the bilayer, our experimental results lie in the shaded region where liquid disordered and liquid ordered regions should coexist,. The 20% cholesterol experiments lie in the triangular region where, in addition to the liquid-ordered and disordered regions, we expect some gel phase. Finally, for the 5% cholesterol experiment, according to the phase diagram, there should only be liquid-disordered regions, together with some gel phase. The phase diagram suggests, then, that, at least for the 5% cholesterol experiment, the identification of the rafts with liquid-ordered regions may not be a perfect analogy, but that the experiments are also detecting and measuring gel phase, which may be enriched in cholesterol. This latter point has been examined, as it is simple to prepare FPE with various kinds of fatty acyl chains. In other words, we are able to prepare the FPE with unsaturated or saturated fatty acids, and this predisposes the FPE to “prefer” the raft regions of the fluid-mosaic membrane. In any event, however, this may not be a problem, as the mechanisms by which the gel phase becomes cholesterol-enriched are probably largely the same as those we postulate in our model. However, the difference in the physical structure of the cholesterol-enriched regions in the two cases may mean that trends in the data are obscured (see below).

It is interesting to note that although the value of the fitting parameter δ varies significantly from one set of experimental results to another, the best-fit value of ϵ is almost the same for all data sets. Recalling the definition of ϵ, this means that the quantity

The trend for the 33% and 20% cholesterol bilayers is as we would expect. The clusters in the 33% bilayer are denser in cholesterol (as is reflected in the A-value); thus, the cholesterol molecules within the clusters are more strongly held together, leading to a higher effective surface tension for the cluster. The anomalous result for the 5% case (higher surface tension than the 20% case) we attribute to the fact, discussed above, that we are probably in a different region of the phase diagram, and so we are not really comparing like with like. The different nature of the molecular packing in the 5% cholesterol bilayer makes it difficult to draw direct comparisons between this experiment and the other two.

The last part of this article was concerned with attempts to formulate a steady-state model to describe the experiments, since all experimental measurements made were at steady state. We saw that the predictions of such thermodynamics-based models could not possibly explain the experimentally observed distributions, showing that the simple thermodynamic approach based on minimization of the Gibbs free energy for the system of clusters is not appropriate for our system, essentially because it does not take into account the different mobilities of the clusters of molecules, which have hugely varying sizes in our system.

In future, we intend to extend the experimental study 1), to measurements taken at different time points (to capture the dynamics of microdomain formation); and 2), to thermodynamically open systems (such as a real cell). In both respects, a formulation such as ours has considerable advantages over a purely equilibrium theory. To make it more biologically realistic, additional elaborations to the mathematical model that we plan include, in the first instance, modeling several lipid types (variations in headgroup identity/chemistry and fatty acyl saturation). We would also include the role of sphingomyelin. This would then put us in a position where we would be able to investigate the role of membrane proteins in modulating these interesting aspects of macroscopic membrane structure. We will also investigate the model for the diffusion-controlled limit (in which clusters react on collision, as opposed to the reaction-limited case considered here, in which the probability of coagulation on collision is relatively small) as an alternative hypothesis. The diffusion-limited problem is much more challenging, since one has to solve a local problem (Laplace’s equation with appropriate boundary conditions) around each cluster to describe the “depletion zone” that surrounds it due to other clusters reacting with it. However, such local solutions in two space dimensions have logarithmic divergence in the far field, making construction of the full solution very delicate.

Appendix A:

A Bifurcation in the steady solution of Eq. (4)

Numerical calculations of the solution to Eq. (4) show that, for ϵ <∼ 0.23, the total mass within the computational domain decreases with time as mass is lost to large values of x. It is impossible to tell from the computations whether this is as a result of a genuine bifurcation, in which mass is lost to infinity, or whether it is a consequence of the finite size of the computational domain. However, by considering the limit of the steady state equations for Eq. (4) as we can show that such a bifurcation is likely as ϵ decreases. We start by considering the toy model for which Here, there is an exact (detailed balance) steady solution to Eq. (4), obtained by setting both integrands in Eq. (4) to zero by writing

In turn, this gives the solution

(16)

where k is a constant that lies in the range 0<k≤1 (note that if k>1, the integral diverges). The value of this constant is determined from the mass conservation condition

Furthermore, it is straightforward to show that It follows that, for it is no longer possible for a solution of the form demonstrated in Eq. (16) to satisfy the mass conservation condition; in terms of the unsteady model, some of the mass is lost to infinity below the bifurcation point at .

As in the full model (Eqs. (4)) the exponential term in the expression for b(x, y) becomes dominant. This suggests looking for a steady solution of the form

In our case, it proves convenient to introduce a further rescaling,

which leaves the steady equations for Eqs. (4) in the form

(17)

(18)

It is obvious from Eq. (18) that g≥O(1/ϵ³) for mass conservation to be satisfied. If this is the case, the dominant terms in Eq. (17) are those that are quadratic in g, and it is fairly clear that, for sufficiently small E and ξ, no balance in which g≥O(1/ϵ³) is possible in Eq. (1)7. This suggests that the real model, like the toy one, exhibits a bifurcation as ϵ decreases through a critical value (our computations suggest ∼0.23) in which mass is lost to infinity.