{'fraction_non_alphanumeric': 0.030885482833149097, 'fraction_numerical': 0.009189751035290624, 'mean_word_length': 4.690199023657529, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Field theoretic study of bilayer membrane fusion: I. Hemifusion mechanism

Self-consistent field theory is used to determine structural and energetic properties of metastable intermediates and unstable transition states involved in the standard stalk mechanism of bilayer membrane fusion. A microscopic model of flexible amphiphilic chains dissolved in hydrophilic solvent is employed to describe these self-assembled structures. We find that the barrier to formation of the initial stalk is much smaller than previously estimated by phenomenological theories. Therefore its creation it is not the rate limiting process. The barrier which is relevant is associated with the rather limited radial expansion of the stalk into a hemifusion diaphragm. It is strongly affected by the architecture of the amphiphile, decreasing as the effective spontaneous curvature of the amphiphile is made more negative. It is also reduced when the tension is increased. At high tension the fusion pore, created when a hole forms in the hemifusion diaphragm, expands without bound. At very low membrane tension, small fusion pores can be trapped in a flickering metastable state. Successful fusion is severely limited by the architecture of the lipids. If the effective spontaneous curvature is not sufficiently negative, fusion does not occur because metastable stalks, whose existence is a seemingly necessary prerequisite, do not form at all. However if the spontaneous curvature is too negative, stalks are so stable that fusion does not occur because the system is unstable either to a phase of stable radial stalks, or to an inverted-hexagonal phase induced by stable linear stalks. Our results on the architecture and tension needed for successful fusion are summarized in a phase diagram.


I. INTRODUCTION
The importance of membrane fusion in biological systems hardly needs to be emphasized. It plays a central role in trafficking within the cell, in the transport of materials out of the cell, as in synaptic vesicles, and in the release of endosome-enclosed external material into the cell, as in viral infection. While proteins carry out many functions leading up to fusion, such as ensuring that a particular vesicle arrives at a particular location, or bringing membranes to be fused into close proximity, there is much evidence that they do not determine the actual fusion mechanism itself. Rather the lipids themselves are responsible for the evolution of the fusion process in which the lipid bilayers undergo topological change (Lee and Lentz, 1998;Lentz et al., 2000;Zimmerberg and Chernomordik, 1999).
The physical description of the fusion process has, until very recently, been carried out using phenomenological theories which describe the membrane in terms of its elastic moduli (Safran, 1994). The application of these theories to fusion has been reviewed recently (Zimmerberg and Chernomordik, 1999). The fusion path which has been considered by these methods is one in which local fluctuations are assumed to cause a rearrangement of lipids in the opposed cis leaflets, resulting in the formation of a stalk (Markin and Kozlov, 1983). To release tension imposed on the membranes by the reduction of the solvent between them, the inner cis layers recede, decreasing their area, and bringing the outer, trans, leaves, into contact. In this way the stalk expands radially to form a hemifusion diaphragm. Creation of a hole in this diaphragm completes formation of the fusion pore.
In recent years, coarse grained models of amphiphiles (Müller et al., 2003a) have been used to provide a microscopic, as opposed to phenomenological, description of membranes. Fusion was studied within two such models. One, in which non-flexible molecules were composed of three segments, was studied by Brownian dynamics simulation (Noguchi and Takasu, 2001), and the other, in which the amphiphiles were modeled as flexible, polymer, chains in solvent, was studied by Monte Carlo simulation (Müller et al., 2002b). Both models showed a markedly different path to fusion than the phenomenological approaches assumed. Along this new path, the creation of the stalk is followed by its non-axially-symmetric growth, i.e. elongation. After the stalk appears, there is a great increase in the rate of creation of holes in either bilayer, and the holes are created near the stalk itself. After a hole forms in one bilayer, the stalk elongates further and surrounds it, forming a hemifusion diaphragm. Formation of a second hole in this diaphragm completes the fusion pore. Alternatively, the second hole in the other membrane can appear before the elongated stalk surrounds the first hole. In this case, the stalk aligns the two holes and surrounds them both forming the fusion pore. This mechanism has also been observed more recently in molecular dynamics simulations (Marrink and Mark, 2003;Stevens et al., 2003). As has been stressed by us recently (Müller et al., 2003b), this alternative mechanism can be distinguished experimentally from the earlier hemifusion mechanism because it predicts at least two phenomena that are not compatible with the earlier hypothesis. The first is lipid mixing between the cis leaf of one bilayer with the trans layer of the other, a phenomena which has been observed (Evans and Lentz, 2002;Lentz et al., 1997). The second is transient leakage through the holes, noted above, which is correlated in space and time with the fusion process. Just such correlated leakage has recently been observed and extensively studied (Frolov et al., 2003).
The Monte Carlo (MC) simulations of the fusion process showed very clearly the nature of the process, and obtained quantitative correlations between leakage and fusion. Unfortunately, simulations are computationally expensive, so that investigation of the fusion process for different molecular architectures and membrane tensions is impractical. Moreover simulations are not well suited for calculating free energy barriers of the fusion intermediates in the mechanism observed over a wide range of tension and amphiphilic architecture. Not only would one like to obtain these barriers, one would also like to compare them to those of the intermediates involved in the original hemifusion mechanism. To do so, we employ a standard model of amphiphilic polymers, which we describe in the next section, and solve it within the framework of Self-Consistent Field Theory (SCFT). In the first paper of this series we examine the original hemifusion mechanism, while in the second we shall consider the new mechanism observed in the simulations.
In Sec. III we present the basic properties of the isolated bilayers and monolayers which result from our calculation. These properties include the spatial distribution of hydrophilic/hydrophobic segments, the area compressibility, the bending rigidity and spontaneous curvature. We compare them to those obtained independently from previous Monte Carlo simulations, and from experiments on liposomes and polymersomes. It should be noted, that these effective macroscopic properties are calculated within our microscopic approach, and are not required as input, in contrast to the common phenomenological descriptions of fusion based on membrane elasticity theory.
In Sec. IV we examine the free energy landscape along the standard hemifusion pathway, one which is shown in Figs. 1 and 2. This is the same path as assumed in phenomenological approaches, but we use a microscopic, molecular, model to calculate along this trajectory the distribution of microscopic components in the system, and the free energy which results. The initial configuration of the system is that of two parallel bilayers ( Fig. 1(a). Hydrophilic portions of the amphiphile are shown dark, hydrophobic portions are shown light, and solvent is shown white. To bring the bilayers into close contact requires energy to reduce the amount of solvent between them. Consequently the free energy per unit area, or tension, of the bilayers increases. Fusion is one possible response of the system to this increased tension. Along the standard pathway, close contact of the cis layers is followed by the creation of an axially-symmetric unstable intermediate, or transition state, shown in cross-section in 1(b), which leads, under conditions detailed below, to a metastable, axially-symmetric stalk, 1(c). If the system is under sufficient tension, this stalk pinches down, passes through another unstable intermediate, 1(d), and expands to form a hemifusion diaphragm 1(e). In order to determine the radius at which a hole forms in this diaphragm to complete the fusion pore, we calculate in Sec. V the free energy of an axially-symmetric fusion pore as a function of its radius. Density profiles of fusion pores are shown in Fig. 2. We assume that when the free energy of the expanding hemifusion diaphragm exceeds that of a fusion pore of the same radius, a hole forms in the diaphragm converting it into a fusion pore, Fig. 2(b). Not surprisingly, the radius of the fusion pore which is formed can not be too small, as in Fig 2(a), for then its free energy would be higher than that of the hemifusion diaphragm. Once the pore has formed, it will expand if the membranes are under tension. Fig. 2(c) shows the profile of an expanded pore, although this one is in a membrane of zero tension. We are able to calculate the free energy of the system at all stages of the pathway for various architectures of the amphiphile and tensions of the membrane.
The most notable results of our calculation are as follows.
i. The free energy barrier to form a stalk, that is, the free energy difference between the initial configuration of Fig. 1(a) and that of the transition state Fig. 1(b), is small, on the order of 10 k B T , much smaller than the estimates of phenomenological theories (Kuzmin et al., 2001).
ii. The more important fusion barrier is encountered on the path between the metastable stalk, Fig. 1(c), through another transition state, Fig. 1 (d), to a small hemifusion diaphragm, Fig. 1(e). The height of this barrier depends strongly on both the effective spontaneous curvature of the amphiphile, as well as on the membrane tension. As expected, the barriers to fusion are reduced as the architecture is changed so as to approach the transition from the lamellar to the inverted-hexagonal phases. The effect of tension on the barriers to fusion is less dramatic, but still very important.
iii. We find that the hemifusion diaphragm does not expand appreciably before converting to the fusion pore.
iv. The small fusion pore formed by rupturing the hemifusion diaphragm can, at low tension, be trapped in a metastable state and not expand further. This result provides an explanation for the flickering fusion events which are observed experimentally. As the tension is increased, however, this metastable state disappears, and the fusion pore, once formed, expands without limit, thus resulting in complete fusion.
v. We observe that the regime of successful fusion is rather severely limited by the architecture of the lipids. If their effective spontaneous curvature is too negative, fusion is pre-empted by the formation of either an invertedhexagonal phase or of a stalk phase. If their curvature is not sufficiently negative, fusion is prevented by the absence of a metastable stalk. We discuss these results further in Sec. VI.

II. THE MODEL
We consider a system consisting of an incompressible mixture of two kinds of polymeric species contained in a volume V . There are n a amphiphilic diblock copolymers, composed of A (hydrophilic) and B (hydrophobic) monomers, while the n s solvent molecules are represented by hydrophilic homopolymers consisting of A segments only. The fraction of hydrophilic monomers in the diblock is denoted f , and the identical polymerization indices of both the copolymer and the homopolymer are denoted N . The hydrophilic and hydrophobic monomers interact with a local repulsion of strength χ, the Flory-Huggins parameter. Provided that these parameters are chosen appropriately, this model is essentially equivalent to that which we simulated earlier (Müller et al., 2003b). We emphasize that the earlier simulation and the calculation presented here are completely independent. Comparison of results obtained from the two different calculations, therefore, are instructive.
The partition function of the system of flexible chains with Gaussian chain statistics can be formulated easily in either canonical or grand canonical ensembles (Matsen, 1995;Schmid, 1998), but is too difficult to be evaluated analytically. Consequently we employ the well established Self-Consistent Field Theory, (SCFT), to obtain a very good approximation to the partition function. This theory has been recently reviewed (Schmid, 1998), so we relegate to the Appendix a brief reminder of the salient features of the approximation. It suffices here to say that the inhomogeneous systems we shall be studying are characterized by a local volume fraction of hydrophilic units, φ A , and of hydrophobic units, φ B . Due to the incompressibility constraint, these two volume fractions add to unity locally. It is convenient to control the relative amounts of amphiphile and of solvent, and therefore the relative amounts of A and B monomers, by an excess chemical potential ∆µ, the difference of the chemical potentials of amphiphile and of solvent. Within this ensemble, the Self-Consistent Field Theory produces a free energy, Ω, which is a function of the temperature, T , the excess chemical potential, ∆µ, the volume V , and, when bilayers are present, the area A which they span. In addition, the free energy is a functional of the local volume fractions φ A and φ B which are obtained as solutions of a set of five non-linear, coupled equations, given in the Appendix, Eqs. 8 to 12. Acceptable solutions are defined by various constraints. For example, we require that all solutions be axially symmetric, an assumption embedded in all previous treatments of the stalk/hemifusion process. We take the axis of symmetry to be the z axis of the standard cylindrical coordinate system, (r, ϕ, z). Further, we require all solutions to be invariant under reflection in the z = 0 plane. Thus φ A (r, ϕ, z) → φ A (r, z) = φ A (r, −z), and similarly for φ B . Other constraints are more interesting. For example, to describe a stalk-like structure, as in Fig. 1, we require that the solution display a connection, along a portion of the z axis, between the hydrophobic portions of the two bilayers. Further, we constrain its radius, R, defined by the condition φ A (R, 0) = φ B (R, 0), to be a value specified by us. In this way we are able to calculate the free energy of stalk-like structures as a function of their radius. To reiterate, once a solution of the set of non-linear coupled equations is found which, for given temperature, excess chemical potential, and volume, satisfies the various constraints imposed, the free energy within the Self-Consistent Field Theory follows directly. We now turn to a description of the results of this procedure for the systems of interest. We begin with those obtained for isolated bilayers.

III. PROPERTIES OF ISOLATED BILAYERS
In this section, we present a range of microscopic and thermodynamic properties of isolated bilayer membranes in excess solvent which follow from our calculation. The most basic structural properties of the bilayer membrane is the distribution of the hydrophilic and hydrophobic segments across it. In Fig. 3 we present the composition profiles obtained within the SCFT approach, shown by solid lines, and compare them to those obtained from our independent simulations (Müller et al., 2003b) which are shown by the symbols. The profiles from simulation are averaged over all configurations. In the Monte Carlo simulation we used amphiphilic polymers consisting of 11 hydrophilic and 21 hydrophobic monomers. This corresponds to a hydrophilic fraction f = 11/32 ≈ 0.34, the value we also used in the SCFT calculations. Fig. 3, corresponds to a system under zero tension. The profiles change quantitatively, but not qualitatively, as the tension is increased.
The overall agreement between the SCFT and averaged MC simulation results is very good. The position and the width of the regions enriched in A (head) and B (tail) segments of the amphiphile and solvent segments are reproduced quantitatively in the SCFT model. The small discrepancies in the A/B interfacial width can be attributed to capillary waves and peristaltic fluctuations present in the MC simulations, but neglected in the SCFT calculations.
Thermodynamic properties of the bilayer can be calculated from the free energy of the system which contains such a bilayer of area A. We denote this free energy Ω m (T, ∆µ, V, A). Similarly, we denote the free energy of the system without the membrane, i.e., a homogeneous amphiphile solution, Ω 0 (T, ∆µ, V ). The difference between these two free energies, in the thermodynamic limit of infinite volume, defines the excess free energy of the membrane: (1) The excess free energy per unit area, in the thermodynamic limit of infinite area, defines the lateral membrane tension In the grand canonical ensemble this tension γ can be related to the temperature and chemical potential by means of the Gibbs-Duhem equation where δs is the excess entropy per unit area, and δσ a is the excess number of amphiphilic molecules per unit area. This relation is quite useful because it shows that one can set the tension to any given value by adjusting the excess chemical potential of amphiphiles at constant temperature. Depending on the applied tension, the thermodynamic behavior of the membrane can be classified into three generic types. For values of the tension which are sufficiently small and positive, the membrane is metastable with respect to rupture, and it is this range of tension which we consider below. For larger positive values of γ, the membrane becomes absolutely unstable to rupture, while for negative γ, the system is unstable to an unlimited increase in the membrane's area, which simply leads to formation of the bulk lamellar phase. In the following, we shall use the dimensionless tension, γ/γ int , where γ int is the interfacial free energy per unit area between coexisting solutions of hydrophobic and hydrophilic homopolymers at the same temperature.
The dependence of the membrane thickness, and membrane tension, γ, on the exchange chemical potential is shown in Fig. 4. The agreement is excellent between the SCFT predictions and the independent simulation results averaged over all configurations. From the membrane tension, the area compressibility modulus, κ A , can be obtained by using any of the following equivalent relations Most of the earlier treatments of membranes relied on elasticity theory, in which a membrane is described solely by its elastic properties, such as the bending modulus, κ M , the saddle-splay modulus, κ G , and the spontaneous curvature, c 0 (Safran, 1994). These moduli are normally taken either from an experimental measurement or from a microscopic theory. The SCFT approach, being based on a microscopic model, allows one to calculate these moduli in a straightforward manner which can be sketched as follows.
One calculates within the SCFT the excess free energies due to small spherical and cylindrical deformations of an interface containing an amphiphilic monolayer. These excess free energies depend upon microscopic parameters, such as the amphiphile hydrophilic fraction f , as well as the curvatures of the deformations. They are then fit to the standard Helfrich Hamiltonian H E of an infinitely thin elastic sheet, which depends upon phenomenological parameters, and the curvatures of the deformations. This Hamiltonian for a saturated, tensionless, monolayer is (Safran, 1994) where M and G are the local mean and Gaussian curvatures of the deformed monolayer. From the fit, one obtains the phenomenological parameters in terms of the microscopic quantities. Details of this procedure can be found in (Matsen, 1999) and (Müller and Gompper, 2002). Because SCFT ignores fluctuations, the moduli obtained are the bare, unrenormalized values. The effect of renormalization is usually small (Peliti and Leibler, 1985) and depends on the lateral length scale of the measurement. Fusion proceeds on a lateral length scale that does not exceed the membrane thickness by a great deal and, on this length scale, we expect the renormalization of the elastic constants to be small. In Fig. 5(a) we show the calculated spontaneous curvature of an amphiphilic monolayer as a function of the hydrophilic fraction f of the amphiphile. The former is a monotonic function of the latter. This is to be expected as the phenomenological spontaneous curvature attempts to capture some of the effects of the differing hydrophilic and hydrophobic volumes in the amphiphile, which are specified by f , on the membrane configurations. The above result provides a direct mapping between the phenomenological property c 0 , which can be measured experimentally (Leikin et al., 1996;Rand et al., 1990) and the microscopic variable f used in our calculations. Our results for the area compressibility, κ A , and elastic moduli κ M and κ G are shown in Fig. 5(b). In agreement with experiments on model lipid systems, they are not very sensitive to changes in amphiphile architecture parameter f or, equivalently, to the spontaneous curvature c 0 . For comparison with our model system, we present in Table I the properties of lipid and amphiphilic diblock membranes determined experimentally.
Lipid membranes are known to exhibit strong mutual repulsion at small separations, usually attributed to so-called hydration forces (Parsegian and Rand, 1994). We calculated the free energy of a system containing two planar bilayers in excess solvent under the condition that the distance between the cis-interfaces be constrained to a given value. This technique has been used before to study monolayer interactions in a similar system (Thompson and Matsen, 2000), and we refer the reader to this work for details.
One finds that there are two generic features of the free energy of the two bilayers as a function of their separation. First at sufficiently small separation, when heads of the amphiphiles come in contact, the membranes experience strong repulsion due to this contact, a repulsion which rises steeply as the separation is further reduced. Second at a somewhat larger separation, there is a very weak attraction between the two membranes. These two features have been exhaustively analyzed (Thompson and Matsen, 2000). It has been determined that the repulsion arises mostly from the direct steric interaction between the head segments of the amphiphiles in the contacting cis monolayers. The weak attraction occurs because the solvent molecules prefer to leave the confining inter-membrane gap in order to increase their conformational entropy. This is simply the well-known depletion effect about which there is a very large literature (Götzelmann et al., 1998). The combination of the short-range repulsion and the depletion effect attraction produces a minimum in the free energy at some separation. The distance at which this minimum occurs, typically on the order of one R g between opposing cis leaves, is the equilibrium separation.
The reasonable description of the properties of self-assembled monolayers and bilayers presented above provides confidence that our model can be used to describe the structural changes which occur in membranes during the fusion process.

IV. ENERGETICS OF THE STALK AND OF ITS RADIAL EXPANSION
The formation and stability of the initial interconnection between the membranes, the stalk itself, has not received much theoretical attention. This is due to the fact that creating the stalk requires the membranes to undergo drastic topological changes, which can not be easily described by continuum elastic models. The common approach has been to assume that two small hydrophobic patches (one in each membrane) are produced by thermal fluctuations in the region of contact. Hydrophobic interactions then drive the connection of these energetically costly regions, and result in the formation of a metastable structure, the stalk (Markin and Kozlov, 1983). The size of the hydrophobic patches and the distance between them can be optimized to minimize the free energy of the unstable transition state. With the use of this strategy, continuum elastic models have estimated the free energy barrier to create a stalk to be about 37k B T (Kuzmin et al., 2001). The SCFT method allows us to obtain solutions for both the unstable transition state to stalk creation as well as for the stalk itself, if it is indeed metastable.
The dependence of the free energy of stalk-like structures on their radius, R, when the bilayers in which they form are under zero tension is shown in Fig. 6(a). It is plotted for various values of the architectural parameter, from f = 0.45, which corresponds to an amphiphile with a very small negative spontaneous curvature, to f = 0.25, an amphiphile with a large negative spontaneous curvature. From Fig. 5(a) it is seen that this f -range includes DOPC and DOPE lipids, which are frequently utilized as components of model lipid membranes in fusion experiments. At f = 0.45 we could not find a stalk-like solution for 2.2 < R/R g < 4. The system would spontaneously rupture in the vicinity of the symmetry axis r = 0, resulting in a fusion pore-like structure.
The extremal points of the free energy of the system with a stalk-like structure with respect to its radius correspond to intermediate structures and transition states. The free energy extrema of the intermediates are saddle points with respect to deformations, and are therefore unstable, while those of the intermediates are local minima, and are therefore metastable. In the inset to Fig. 6(a) we identify three states with extremal free energies. Depending on the values of the architecture parameter f and the membrane tension γ, we find the following solutions that play central roles in the description of the membrane fusion process: the transition state S 0 between the unperturbed bilayers and the metastable stalk, the metastable stalk itself, S 1 , and finally, a transition state which occurs as the initial small stalk is radially expanded into the hemifusion diaphragm. It is denoted S 2 . It is clear from Fig. 6(a), that the metastable stalk solution S 1 , and therefore the transition state S 0 , exists only for sufficiently small values of the hydrophilic fraction f . At zero membrane tension this minimal f is about 0.36, which corresponds to the spontaneous curvature c 0 d = −0.9, roughly that of a 1 : 1 mixture of DOPE and DOPC lipids. To our knowledge, this is the first theoretical prediction of the metastability of the stalk itself. Existence of the stalk intermediate is crucial for the fusion dynamics, since it separates the fusion process into two activated stages: creation of the stalk and its further expansion. We argue, therefore, that for those conditions, specified by f and γ, under which there is no S 1 solution, fusion will be considerably slower or, perhaps, impossible. Fig. 6(a) also shows that, at small stalk radius, R, the free energy barrier to create the stalk, which corresponds to the state S 0 , does not exceed 5k B T . This value is much smaller than the 37k B T predicted by the phenomenological calculation (Kuzmin et al., 2001). Of course, one can ask how applicable our results for block copolymer membranes are to biological ones. If we take as a natural measure of energy the dimensionless quantity γ int d 2 /k B T , where d is the membrane thickness, we obtain a value of about 60 in our system, a factor of 2.5 less than that (≈ 150) for a typical biological bilayer. This suggests that our energies should be scaled up by a factor of about 2.5 in order to compare them with those occurring in biological membranes. However even after doing so, they are still much smaller than the estimates of the phenomenological theories. We infer from this result that the phenomenological approaches based on the continuum elastic description are not accurate in describing the drastic changes in membrane conformations on the length scales comparable to the membrane thickness.
In Fig. 1 (b), we show the calculated density profiles of the different segments in the unstable transition state structure S 0 at f = 0.35 and γ/γ int = 0. As can be seen, the interface between A and B segments is extremely curved and the stalk radius R = 0.6R g is much smaller than the membrane hydrophobic core thickness. For clarity, we show only the majority component at each point. The interfaces between the different components appear to be sharp in such a graphical representation, but actually they are relatively diffuse, as can be seen in Fig. 3.
The behavior of the free energy of the stalk-like structure as it expands into a hemifusion diaphragm under nonzero tension is shown in Fig. 6(b), and one clearly sees a second local maximum at R/R g ≈ 3 which corresponds to the second unstable transition state, S 2 . The maximum results from the competition between the elimination of energetically-costly bilayer area, which reduces the free energy as −γπR 2 , and the creation of diaphragm circumference which increases it as 2πλ tri R. Here λ tri is a line tension. Fig. 6(b) shows that this second maximum is, in general, much greater than that encountered in creating the stalk itself. Therefore we infer that crossing this barrier is the rate limiting step of the hemifusion mechanism. As the tension increases, the height of this barrier decreases. Thus at very high tension, on the order of 0.5γ int for f = 0.3, this local maximum disappears entirely leading to stalk expansion without any barrier at all. The rather strong dependence of this barrier height on tension contrasts with the behavior of the free energy of the metastable stalk itself, which, from the figure, is seen to depend only weakly on the tension. We also note that the minimum in the free energy corresponding to the metastable stalk is exceedingly shallow. This means that these interconnections are easily reversible, and would constantly fluctuate in size.
We recapitulate, in Fig. 7, the dependence of the free energy of the metastable stalk S 1 and of the unstable transition state S 2 on the tension, γ, and on the architectural parameter f . The free energy of the metastable stalk (S 1 ) varies greatly with f , and decreases substantially for smaller f , i.e., as the spontaneous curvature of the amphiphile becomes more negative. Although our calculation applies to membranes composed of a single amphiphile, and not a mixture, it clearly strengthens the argument that one role of such negative curvature lipids as phosphatidylethanolamine, present in the plasma membrane, is to make metastable and thermally accessible the formation of stalks, which are necessary to begin the fusion process (de Kruijff, 1997) . For sufficiently small values of the architectural parameter f , the free energy of the metastable stalk actually becomes lower than that of the unperturbed bilayers. Presumably this leads to the formation of the thermodynamically stable "stalk phase" recently realized experimentally in a lipid system (Yang and Huang, 2002).
In summary, both an increase in the membrane tension γ and a decrease in the hydrophilic fraction f favor stalk expansion into the hemifusion diaphragm. The density profiles in the transition state between the metastable stalk and the hemifusion diaphragm is shown in Fig. 1(d), and that of the unstable, expanding, hemifusion diaphragm in Fig. 1(e). Note how thin is the hydrophobic region on the axis of symmetry in (d) and next to the triple junction in (e) compared to the thickness of the bilayers away from the diaphragm.

V. FORMATION AND EXPANSION OF FUSION PORE
The hemifusion diaphragm is a possible intermediate along the path to fusion, but for complete fusion to occur, a hole must nucleate in the diaphragm leading to the formation of the final fusion pore. Here we consider the energetics of the fusion pore. The manner in which we obtain solutions corresponding to the pore and determine its free energy are similar to those for the stalk described above. A density profile of a fusion pore is shown in Fig. 2(a). In the plane of mirror symmetry between the two bilayers, there are two radii at which the volume fraction of hydrophobic and hydrophilic elements are equal, i.e. for which φ A (r, 0) − φ B (r, 0) = 0. We define the radius, R, of the fusion pore as the greater of these two distances. This choice is consistent with the definition of the radius of a stalk-like structure, the only radius at which φ A (r, 0) − φ B (r, 0) = 0, as can be seen from the following: if one begins with a stalk-like structure of radius R and introduces a hole into it to create a pore, then the radius of the pore, according to the above definition, is also R, as it should be. We are able, therefore, to compare directly the free energies of the pore and stalk-like structures with the same radius.
The free energy of the fusion pore in a system with fixed architecture, f = 0.35, but under various tensions is shown in Fig. 8. As the tension on the membrane is increased from zero, the free energy of the pore decreases as −2γπR 2 for large radius R, similar to the free energy decrease of the hemifusion diaphragm. The factor of two arises because pore expansion eliminates both membranes, whereas hemifusion expansion removes only one.
The free energies of the fusion pores, shown in solid lines, are compared with those of the stalk-like structures under the same conditions, which were shown previously in Fig. 6(b). We assume that a hole forms in the hemifusion diaphragm converting the diaphragm into a fusion pore when the free energy of the pore becomes lower than that of the stalk-like structure. The radius at which the free energies of the expanding stalk and of the fusion pore become equal is R ≈ 2.5 − 3.5R g . Of course, there is not a sharp transition from the stalk to the pore, because the free energies are finite. There is rather a region of radii where the transition from the stalk to the fusion pore occurs, and which is characterized by a free energy difference of the order of thermal fluctuations k B T . Because the stalk and pore structures are so similar at such small values of radius, it is likely that there is only a very small free energy barrier associated with rupture of the hemifusion diaphragm, which converts it to the fusion pore. Therefore, our calculations indicate that the hemifusion diaphragm would hardly expand before the fusion pore would form. This agrees with the conclusion of a recent phenomenological calculation (Kuzmin et al., 2001). Fig. 8 shows that, except at low tension, the hole forms in the diaphragm after the barrier to diaphragm expansion has been crossed. Therefore within the standard hemifusion mechanism, the barrier to hemifusion expansion is, in fact, the major barrier to fusion.
At very low tensions, our results show that the most important barrier to fusion is no longer that governing the expansion of the hemifusion diaphragm, but becomes that associated with the expansion of the fusion pore itself. As a consequence of the large barrier to expansion of the fusion pore, pores of small radius, R ≈ 3.4R g , become metastable for most architectures for tensions γ/γ 0 less than about 0.1. This metastability persists to zero tensions, as seen from Fig. 9. Note that the potential minimum of the metastable pore is quite shallow, therefore one would expect that thermal fluctuations will easily cause it to expand and contract, or flicker, about the minimum free energy configuration with an amplitude of the order of 1R g . Note also that the activation barrier to reseal the fusion pore is quite substantial, on the order of 5k B T , which would translate to about 13k B T in lipid systems. This means that flickering pores can be long-lived metastable states and, depending on the tension, might either reseal or expand. In retrospect, it should not be surprising that a small fusion pore is metastable in a tensionless membrane. Its free energy must increase linearly with its radius due to the line tension cost of its circumference. On the other hand as its radius decreases, the pore eventually becomes sufficiently small so that the inner sides of the pore come into contact, repelling each other and causing the free energy to again increase. Thus there must be a minimum at some intermediate distance. If the radius is decreased below some critical value, the compressed pore structure becomes unstable to the stalk-like structure. This instability is shown with arrows in Fig. 9. The density profile of such a marginally stable pore occurring in a membrane with f = 0.35 and R = 2.4R g is shown in Fig. 2(a). Fig. 9 shows that the small flickering pore would be thermally excited from a metastable stalk for most of the architectures shown, f = 0.40, 0.35, and 0.30. As the architecture changes such that f decreases still further, the stalk free energy becomes negative meaning that it is favorable to create many of them which leads to the formation of a stalk phase. This occurs for non-zero tension as well.
Lastly we note from Fig. 9 that at very low tension and very small values of the architectural parameter, f = 0.25, corresponding to a large and negative spontaneous curvature, both the fusion pore and the stalk-like structure are unstable to another structure, the inverted micellar intermediate (IMI), introduced by Siegel (Siegel, 1986). Interestingly this configuration was suggested previously as a possible player in the fusion process, but has been largely neglected because free energy estimates obtained from elastic approaches were prohibitively high. As will be shown in our second paper, the IMI and related linear stalk intermediates are very important in describing fusion.

VI. DISCUSSION
We have carried out a self-consistent field study of the fusion of membranes consisting of flexible block copolymers in a solvent of homopolymer. The main purpose of this paper was to evaluate the free energy barriers encountered within the standard hemifusion, or stalk, mechanism. We summarize our major findings in Fig. 10 in a form of a "phase diagram" in the parameter space of amphiphile architecture and tension.
We find that as the architecture of the amphiphile f is changed so that its spontaneous curvature decreases from zero and becomes more negative, fusion is enhanced. There are at least two reasons for this. First, the initial stalk becomes metastable only if f is sufficiently small, (i.e., the spontaneous curvature is sufficiently negative). This is the case in the region to the left of the nearly vertical black dashed line in Fig. 10. Presence of the metastable stalk intermediate is crucial for fast fusion, because its formation represents the first of at least two activated steps in the fusion process. Within the hemifusion mechanism, the second activated step is the radial expansion of the stalk into a hemifusion diaphragm. We expect that in the region in which the stalk is not metastable, fusion would be extremely slow, if not impossible.
We also find that for a sufficiently small f the stalk intermediate becomes absolutely stable with respect to the system with unperturbed membranes. There is a small region (shown in the figure between the dotted and solid lines) in which the stalk is not only stable with respect to the unperturbed system, but is more stable than any other intermediate we have considered. Consequently we predict formation of a "stalk phase" in this region. This result is in accord with recent experiments on model lipid systems (Yang and Huang, 2002). As the spontaneous curvature is made even more negative, linear stalks, which are precursors to the formation of an inverted-hexagonal phase, become even more stable than radial stalks. Hence the stalk phase becomes unstable with respect to the inverted hexagonal phase.
We remark that, according to our calculations on bilayers comprised of a single amphiphile, the variation of architectures within which successful fusion can occur is quite small. One way that a system composed of many different amphiphiles can ensure the appropriate effective architecture is to employ a mixture of lipids comprised of those with small spontaneous curvatures, or lamellar formers, and those with larger negative spontaneous curvatures, or non lamellar-formers. Further, to remain within the small region of successful fusion, the ratio of these different kinds of lipids must be regulated, as indeed it is (Morein et al., 1996;Wieslander and Karlsson, 1997).
We also note that the restrictions we have found on the possible variation of curvature, if fusion is to be successful, arise from stability arguments involving the stalk itself. Therefore any fusion mechanism which begins with the formation of a stalk will have this same phase diagram even though the subsequent path to fusion may differ radically from the standard hemifusion mechanism.
Within the range of the parameters where fusion can occur, we find that the activation barriers of the three major stages of the standard hemifusion process, which are stalk formation, stalk expansion, and pore formation, are affected differently by changes in architecture and in tension. The first activation barrier is associated with creation of the initial metastable stalk. Its height is essentially independent of the membrane tension and increases only very weakly as f is increased. Stalk creation is apparently not a rate limiting step, because the corresponding barrier does not exceed 5k B T in our model, or 13k B T if we extrapolate our results to lipid systems. This value is much lower than previous estimates of phenomenological theories (Kuzmin et al., 2001).
In the second stage, during radial stalk expansion, the corresponding barrier is a very sensitive function of the amphiphile architecture, as can be seen in Fig. 7. At small tension it ranges from ≈ 10k B T for amphiphiles with small f up to ≈ 25k B T for more symmetric amphiphiles. Again, this range translates into 25 − 63k B T for biological bilayers. These observations are consistent with well known fusion enhancement effects of lipids with large and negative spontaneous curvature, such as DOPE, when they are added to fusing membranes (Chernomordik, 1996;Chernomordik et al., 1995). The effect of membrane tension on the stalk expansion barrier has also been determined. We find that increased tension lowers the activation barrier to radial stalk expansion and to pore formation. This is in accord with experiment (Monck et al., 1990). We predict that this barrier to fusion can, in principle, always be reduced and eliminated by sufficient tension. However, if the architecture of the amphiphiles is unfavorable, resulting in a very high barrier at zero tension, the tension needed to eliminate this barrier can be prohibitively high, and fusion will be pre-empted by membrane rupture. Results on the thermodynamics of membrane rupture will be presented in the second paper of this series.
The third stage of the process, nucleation of a hole in the diaphragm to convert the hemifusion diaphragm into a fusion pore and the pore's expansion, appears to present no additional barrier to fusion at most non-zero tensions. Thus we find that the largest barrier encountered in this standard hemifusion mechanism is that associated with the expansion of the stalk into a hemifusion diaphragm. The diaphragm does not expand very much, however, before this maximum barrier is reached.
At very low tensions, the controlling barrier is that associated with the expansion of the fusion pore. This leads to the prediction of the transient stability of small fusion pores. We believe that these structures correspond to "flickering pores" observed experimentally in lipid bilayer fusion (Chanturiya et al., 1997;Fernandez et al., 1984;Spruce et al., 1990).
In spite of the agreement between our results on the hemifusion mechanism and the experimental observations mentioned above, there remain experimental observations which this hypothesis does not explain. First, lipid mixing is observed to occur not only between the cis monolayers, but also between cis and trans layers (Lentz et al., 1997). Second, transient leakage is observed, and it is correlated spatially and temporally with fusion (Frolov et al., 2003). In addition to these experimental observations, simulation studies have revealed, both in planar bilayer fusion and in vesicular fusion, a very different pathway subsequent to the formation of the initial metastable stalk. To understand these discrepancies, we have performed calculations similar to the ones presented here for the alternative mechanism proposed recently (Müller et al., 2002b(Müller et al., , 2003bNoguchi and Takasu, 2001). These results will be presented in the second paper of this series. gratefully acknowledged VII. APPENDIX A. Self-consistent field theory for polymer systems The first step in any field-theoretic approach to such systems is to convert the partition sum over all possible molecular configurations into an integration over configurations of corresponding, smooth, collective variables, the density and chemical potential fields. The derivation of the effective field-theoretic Hamiltonian of our polymer model follows a standard prescription, so we give only the final expression here: Here z ≡ exp ∆µ/k B T is the relative activity of the amphiphiles, where ∆µ = µ a − µ s , the difference in bulk chemical potentials of the amphiphile and the solvent. There is only one independent chemical potential because the liquid is assumed to be incompressible. The fixed number density of polymer chains, (n a + n s )/V is denoted Φ. Thus the total number density of chains is fixed, but the relative amounts of amphiphile and solvent chains is controlled by the chemical potential difference ∆µ. The Flory interaction parameter, χ, is inversely proportional to the temperature. The Q s [w A ] and Q a [w A , w B ] are the single chain partition functions of the solvent and amphiphile molecules subjected to the local chemical potential fields w A (r) and w B (r), which act on the A and B segments respectively. The local volume fractions of A and B monomers are given by φ A (r) and φ B (r). The "local pressure" ξ(r) is the Lagrange multiplier field introduced to enforce the incompressibility condition. The mean field approximation, which is at the heart of the SCFT approach, amounts to extremizing the fieldtheoretic Hamiltonian with respect to all the fields upon which it depends. Fluctuations about this extremum are ignored. In systems of relatively long chains, in which composition fluctuations are small, the results of this approximation are very good indeed (Bates and Fredrickson, 1999).
It can be shown that the field configurations that correspond to stationary points of H[T, ∆µ, V, A; w A (r), w B (r), φ A (r), φ B (r), ξ(r)], denoted in the following by an over-bar, satisfy the following set of coupled non-linear equations: where the single chain propagators q c (r, s),q † c (r, s), and q s (r, s) satisfy the usual modified diffusion equations for the flexible polymer chains with the Gaussian statistics; e.g. for the homopolymer solvent it is Here R g is the radius of gyration of the unperturbed Gaussian polymer. The value of the free energy of the stationary configuration is given simply by The system of non-linear equations (8)-(12) together with the equations for the propagators (13) can be solved numerically in real space by a straightforward relaxational iterative algorithm (Drolet and Fredrickson, 1999). The real-space approach has far more flexibility in studying localized structures, such as the fusion intermediates, than does the more numerically optimized spectral approach (Matsen and Schick, 1994). The recently proposed pseudo-spectral techniques (Rasmussen and Kalosakas, 2002) appear to be very efficient, but they rely heavily on the numerical fast Fourier transforms, which are not available for the cylindrical coordinates, r = (r, ϕ, z), we use in our calculations. Because the local volume fractions are required to be axially symmetric, φ A (r, z), φ B (r, z), and to be symmetric under reflection in the z = 0 plane, the problem is two-dimensional and need only be solved in one quadrant. We impose reflecting boundary conditions at z = 0, z = z max , r = 0 and r = r max , where z max and r max determine the size of the computational cell. They were set to 8R g and 15R g respectively. We discretized all the fields on a uniform lattice with resolution ∆r = ∆z = 0.1 − 0.05R g . In solving the diffusion equation (13) we used contour length discretization of ∆s = 0.01 − 0.001. This gave us an accuracy of no less than 0.1k B T in the free energy. On a one GHz, Pentium 3, work station, approximately five minutes were required to obtain convergence for a configuration, φ A (r, z) and φ B (r, z), at a given value of f and γ.
To describe a stalk-like structure, the solutions must smoothly connect the hydrophobic regions of the two bilayers. The radius, R, of the stalk, again defined by φ A (R, 0) − φ B (R, 0) = 0, must take a value that we specify. To facilitate finding such a solution we apply, at the beginning of our calculations, auxiliary external fields which favor hydrophobic segments along the axis of symmetry between the two trans leaves. After the system assembles into the desired structure, these auxiliary external fields are switched off, and the solution is allowed to relax. Far from the axis of symmetry, the membranes reach their equilibrium separation. Again this minimum is due to the short-range repulsion and the attractive depletion interaction between the bilayers. Once a solution satisfying these constraints is obtained, the free energy of the stalk-like structure is then calculated. It is the difference, at constant chemical potential (or, equivalently, constant tension), between the free energy of the system with two bilayers connected by the stalk-like structure, and that of the two unperturbed bilayers without the interconnection.

B. Reaction coordinate constraint
The SCFT strategy and its numerical implementation along the lines presented above are capable, in practice, only of identifying thermodynamic, locally stable, configurations of the system. To clarify this point, consider the following construction. Suppose, we first extremize H with respect to the chemical potential fields w A (r) and w B (r), and the incompressibility field ξ(r). Then, at least in principle, we obtain a free energy functional that depends only on the physical density fields φ A (r) and φ B (r): where we have suppressed the dependence on the thermodynamic variables, T, ∆µ, V and A. An extremum of F [φ A , φ B ] is the system's free energy Ω m . This extremum corresponds to a thermodynamic, locally stable, state if and only if the matrix of the second derivatives δ 2 F /δφ A (r)δφ B (r ′ ) is positive definite in that configuration, i.e., the locally stable configurations correspond to the minima of the free energy density functional F [φ A , φ B ]. To distinguish them from other kinds of solutions, we will refer to these locally stable structures as intermediates.
In studying an activated process, such as membrane fusion, we need to know not only the free energy of the intermediates along some reaction path, but also the properties of the transition states, which correspond to saddle points of the free energy functional F [φ A , φ B ]. Unfortunately, finding a saddle point of a functional poses a serious numerical problem. In particular, commonly used relaxational algorithms prove to be inadequate for this task, because they rely on local stability around a solution, which is obviously lacking at the saddle point. Newton-Raphson type algorithms are capable of finding any extremal points as long as the initial configuration is within the basin of attraction of that point. Unfortunately, we do not know the location of a saddle point in advance, so these methods are also not very practical.
In some cases one can identify the unstable directions of the functional and stabilize it by applying a suitable constraint. As an example, consider a transition state between the stalk and the hemifusion diaphragm. We treat this situation in detail in the main text and use it here simply for illustration. We confine the possible solutions to be axially symmetric and to have a mirror symmetry with respect to z = 0 plane because the unstable mode of either the stalk or diaphragm can be associated with an overall radial contraction or expansion of the structure. We further constrain the solution by requiring that, in the z = 0 plane, the A/B interface (i.e., the locus of points at which φ A (r) − φ B (r) = 0) be located on a circle of some specified radius R. This radius plays the role of a reaction coordinate. To impose the constraint, we employ a Lagrange multiplier, ψ R (Matsen, 1999). The corresponding constrained field-theoretic Hamiltonian is then given by The first two SCFT Eqns. (8) and (9) should be modified accordingly: and the third equation expressing the local density constraint is supplemented by the additional local constraint: where, again, (r, z) are cylindrical coordinates. The Lagrange multiplier ψ R plays a role of the local chemical potential that couples to the density fields. The solution of the SCFT equations optimizes the free energy with respect to ψ R as well as the other fields. The same relaxational iterative approach as that used for the non-constrained case proved to be efficient for solving this modified problem. In general, for an arbitrary position of the constraint R, the value of the Lagrange multiplier ψ R is non-zero, which means that the corresponding field configuration is not a solution of the original non-constrained system. Nevertheless, at the points where ψ R is zero, the constrained and non-constrained sets of equations become identical, i.e., an extremum of the constrained free energy functional is also an extremum of the non-constrained functional.
The free energy of the constrained system, F R , quite generally, satisfies the following "force balance" equation: The condition ψ R = 0 is equivalent to dF R /dR = 0, so that the extremal points of F R as a function of R are also the extremal points of F [φ A , φ B ] as a function of R. Clearly, a minimum of free energy F R of the restricted system corresponds to a minimum of F [φ A , φ B ], whereas a maximum of F R corresponds to a saddle-point configuration of F [φ A , φ B ]. By scanning a range of R values, we can identify all the metastable intermediates and unstable transition states along a particular path. Moreover even those constrained configurations for which ψ R = 0 have a transparent physical interpretation, and provide a wealth of additional information, inaccessible by the non-constrained calculations. For other applications of similar constraints we refer the reader to the literature (Duque, 2003;Matsen, 1999;Müller et al., 2002a).

C. Model parameters
To make a direct comparison with our previous independent Monte Carlo simulations (Müller et al., 2003b), we match corresponding model parameters. The length scale in SCFT calculations is usually set by the polymer radius of gyration R g . For the polymers in the simulations, the radius of gyration was found to be R g = 6.93u, where u is the spacing of the cubic lattice into which the simulation volume is divided. Because of the incompressibility constraint, the volume per polymer, V /(n a + n s ), with n a and n s the number of amphiphilic and solvent polymers respectively, enters the SCFT free energy only as a multiplicative factor and in a dimensionless ratio V /[(n a + n s )R 3 g ]. As this ratio was taken to be 1.54 in the simulations, we take the same value in evaluating the SCFT free energy. The energy scale in SCFT is set by the product of the Flory interaction parameter χ and the polymerization index N . It has been shown previously that the choice χN = 30 corresponds to the simulated system (Müller et al., 2003b).   Discher et al., 1999; and on lipids from (a): Rand and Parsegian, 1989, (b): Rand et al., 1990, (c): Chen and Rand, 1997, and (d): Leikin et al., 1996  The radius of the stalk is R = 0.6Rg. In general the stalk radius R is defined by the condition on the local volume fractions φA(R, 0)−φB(R, 0) = 0. (c) The metastable stalk itself, which we label S1. Its radius is R = 1.2Rg. (d) The unstable transition state between the metastable stalk and the hemifusion diaphragm. This transition state is denoted S2. Its radius is R = 2.8Rg . (e) A small hemifusion diaphragm of radius R = 3.4Rg. The radius is shown explicitly. . This result is compared to those of an independent simulation averaged over all configurations, (solid circles). Also shown is the excess free energy per unit area, or tension, of the bilayer, scale to the right, as a function of the exchange chemical potential. Results obtained in our SCFT calculation are shown by the solid line, and those obtained in an independent simulation are shown by the dotted line. The error bar shows the uncertainty in the exchange chemical potential at which the bilayer is without tension as determined in the simulations. FIG. 6 (a) The free energy, F , of the stalk-like structure connecting bilayers of fixed tension, zero, is shown for several different values of the amphiphile's hydrophilic fraction f . In the inset we identify the metastable stalk, S1, the transition state, S0, between the system with no stalk at all and with this metastable stalk, and the transition state, S2, between the metastable stalk and a hemifusion diaphragm. The architectural parameter is f = 0.30 for this inset. No stable stalk solutions were found for f = 0.45 in the region shown with dashed lines. They were unstable to pore formation. (b) The free energy of the expanding stalk-like structure connecting bilayers of amphiphiles with fixed architectural parameter f = 0.35 is shown for several different bilayer tensions. These tensions, γ/γint, are shown next to each curve. Under the assumption that the stalk-like structure converts to a fusion pore when their free energies cross, one can read off the barrier between metastable stalk and formation of the fusion pore. The membranes are comprised of amphiphiles with fixed architecture, f = 0.35, and are under various tensions. Note that the fusion barrier decreases with increasing tension. The free energies, F , of a fusion pore (solid lines) and of a stalk (dashed lines) of radius R are shown. In contrast to Fig. 8, the membranes here are under zero tension, and are comprised of amphiphiles with various values of f . The instability of the fusion pores at small radius is indicated by arrows. For f = 0.3, 0.35, and 0.4, the stalk-like structure converts into a pore when it expands to a radius R ≈ 2.4Rg at which the free energies of stalk-like structure and pore are equal. For the system composed of amphiphiles of f = 0.25, however, the stalk-like structure converts at R ≈ 2Rg into an inverted micellar intermediate, (IMI), whose free energy is shown by the dotted line. The fusion pore is unstable to this IMI intermediate when its radius decreases to R ≈ 3.4Rg. Thus the IMI is the most stable structure under these conditions. FIG. 10 A "phase diagram" of the hemifusion process in the hydrophilic fraction-tension, (f, γ), plane. Circles show points at which previous, independent, simulations were performed by us. Successful fusion can occur within the unshaded region. As the tension, γ, decreases to zero, the barrier to expansion of the pore increases without limit as does the time for fusion. As the right-hand boundary is approached, the stalk loses its metastability causing fusion to be extremely slow. As the left-hand boundary is approached, the boundaries to fusion are reduced, as is the time for fusion, but the process is eventually pre-empted due to the stability either of radial stalks, forming the stalk phase, or linear stalks, forming the inverted hexagonal phase.