Exploring thermodynamic stability of the stalk fusion-intermediate with three-dimensional self-consistent fi eld theory calculations

The prospects of compressible Self-Consistent Field (SCF) theory schemes for describing structures in amphiphilic membranes are illustrated by considering the thermodynamic stability of hourglass-shaped, hydrophobic connections (stalks) between apposed bilayers. The membranes are represented by a coarse-grained, solvent-free model. We represent the chain architecture by a Gaussian-thread representation of the chain architecture and capture the non-bonded interactions with a functional, which is of third-order in the densities of the hydrophilic and the hydrophobic segments. Using a three dimensional real-space scheme, we study the thermodynamic stability of the stalk with respect to two planar apposing bilayers as a function of membrane tension and molecular asymmetry. The structure and thermodynamics predicted by SCF theory agree very well with particle-based simulations, which include fluctuations. We discuss how the longer-range perturbations of the membrane induced by the stalk can affect thermodynamic properties.


Introduction
Membrane fusion is a prototypical collective transformation of lipid membranes that is crucial for a variety of basic biological processes, 1-3 e.g., endo-and exocytosis or synaptic release.][6][7][8][9][10][11] Recent experiments 12 suggest that the shape of thermodynamically stable stalks is universal, i.e. it does not considerably vary with lipid species.Indeed, the rhombohedral phase, in which passages between planar lamellae arrange on an hexagonal lattice, has been observed both in aqueous lipids [12][13][14][15][16] as well as in synthetic block copolymer systems. 17,18In marked contrast to the shape, however, the excess free energy of a stalk is a pronounced function of the lipid architecture. 19,20he free energy of stalks has attracted abiding interest and has been studied by phenomenological calculations using curvature Hamiltonian 1,21 and computer simulations 8,10,20 alike.While the former invoke uncontrolled approximations, e.g., small curvature expansion, which may be quite inaccurate for such bent fusion intermediates, and account for the arrangement and molecular shape changes of lipids in the highly curved bilayers only indirectly, the latter are computationally expensive.Thus a systematic exploration of the lipid architecture and complex membrane topologies remains a methodological challenge.
0][31][32] In this context, however, previous applications were limited to two-dimensional calculations for a model mixture of block copolymers and homopolymers.4][35][36] Here we use a three-dimensional real-space approach and apply it to a coarse-grained solvent-free model, 34 which has recently been studied by molecular simulation. 20,37,38We investigate the structure and thermodynamics of stalks and compare the results of SCF theory quantitatively to simulation data. 20We emphasize that the computational scheme, however, can also be applied to more complex threedimensional structures, e.g., the stalk-pore fusion intermediates.

Model and computational approach
The molecular architecture of the amphiphiles is represented through a continuum analog of the bead-spring model employed in the particle-based simulations, i.e., they are considered as Gaussian threads.Within this description the effective Hamiltonian of the bonded interactions reads: where R e denotes the average end-to-end distance of an amphiphile in the unperturbed coil regime and will be used in the following to set the length scale.The non-bonded interactions are captured through a local density functional: where rA (r) and rB (r) are the instantaneous local densities of the hydrophobic and hydrophilic blocks dened as: In the above equation, g a (s) ¼ 1 if the s-th monomer is of type a and g a (s) ¼ 0 if otherwise.Here a ¼ A if s # f and a ¼ B if otherwise; f expresses the molecular asymmetry dened as f ¼ N A /N, where N A stands for the number of hydrophobic monomers in an amphiphile made of N monomers in total.The total number of amphiphiles in the system is denoted by n.
We use the same parameterization as in the particle-based simulation, 20 setting the virial coefficients in eqn (2) to v AA ¼ À5.15, v BB ¼ À0.01, v AB ¼ À1.775, w AAA ¼ w AAB ¼ w ABB ¼ 0.095625, and w BBB ¼ 0. The strategy used for their identication is discussed in ref. 34.Here we briey mention that the coefficients v AA and w AAA are chosen aer considering the equation of state of a homopolymer melt (dened from eqn (2) assuming a one-component system) consisting only of A-blocks and requiring that it has a certain bulk density r A and compressibility; the current choice corresponds to r A ¼ 40 amphiphiles per R e 3 .v AB dictates the repulsion between hydrophilic and hydrophobic segments and can be adjusted to match the oil-water interface tension.Setting v AA ( v BB implicitly accounts for the fact that the solvent is signicantly preferential towards the hydrophilic B heads, and v BB additionally controls the strengths of hydration repulsion in our solvent-free model.The SCF scheme for the solvent-free, compressible model of amphiphilic membranes is a straightforward generalization of standard SCF formulations for multicomponent polymeric systems. 23,25,27,28One decouples the non-bonded interactions, expressing the interactions of the A and the B parts of the amphiphiles with their surroundings through the mean elds, W a (r) (a ¼ A, B) calculated as: In the above equation, r A (r) and r B (r) stand for the average local densities of the two block species and we used the fact that the virial coefficients in our case are symmetric with respect to index permutation to obtain a compact expression.
The self-consistency "loop" is obtained aer linking the average local densities of the A and B blocks to the conformations of the amphiphiles in the mean elds.These can be characterized through the "propagators" q(r, s) and q + (r, s), proportional to the partition functions of two subchains of contour lengths s and 1 À s, connected at point r of space, while their corresponding free ends are located anywhere.The propagators are obtained from the solution of the Edwards diffusion equation, which for q(r, s) reads: A similar expression holds for q + (r, s) with the difference 39 that the right-hand side is multiplied by À1, while the initial condition is: q + (r, 1) ¼ 1.The average local densities are calculated from the propagators as: where Q ¼ Ð q(r,1)dr/R e 3 is the partition function of a single amphiphile in the mean elds.
Apart from describing conformational properties and the spatial distribution of different components, SCF theory allows us to compare the thermodynamic stability of different structures because it yields the thermodynamic potential of the system in a given morphology.In our case, where the computations are performed in the nVT ensemble, this will be the Helmholtz free energy, F(n, V, T), calculated within the mean eld approximation as: Although straightforward conceptually, when exploring the stability of the fusion intermediate, the solution of the set of the SCF equations, i.e., eqn ( 4)-( 6), in three dimensions requires extreme numerical precision.To meet this numerical challenge, we invoke a three-dimensional real-space scheme, 26 comprising an iterative relaxation scheme.Namely, one starts from an initial conguration of the elds, W (1) a (r).These are substituted into the Edwards equation, eqn (5), which is solved numerically using an alternating direction implicit algorithm (ADI) to obtain the propagators q and q + .Subsequently, the densities are calculated from eqn (6).A new value for the eld is established as where W a (r) is obtained by substituting the calculated densities into eqn (4).The magnitude of the relaxational parameter is set here to l ¼ 0.02.The new elds serve as an input for the next iteration step and the whole procedure, described previously, is repeated again.We consider the iterative procedure as converged when the relative accuracy of the scheme 40 is G a $ 10 À10 .During the calculations reecting boundary conditions are employed in all three directions; the dimensions of the box are L x ¼ 3.25R e , L y ¼ 3R e , and L z ¼ 3R e .To reach high accuracy, in the ADI scheme and during the calculation of the densities (cf., eqn (6)), the contour of the amphiphile is discretized with a step ds z 3 Â 10 À3 while the volume is discretized via a cubic lattice with spacing dx z 0.025R e .

Results and discussion
To obtain an optimized stalk morphology, we start the iterations from an initial conguration of the elds, which have a very crude stalk-like pattern.The axis of the stalk is oriented along the longest edge of the box, L x .In the initial eld conguration W A (r) takes two values, one for the hydrophobic and one for the hydrophilic region, while W B (r) is constant everywhere.Here our choice of initial eld conguration sets the distance of the apposed bilayers in the stalk morphology to be around 2.6R e (as measured from their center planes).Due to the reecting boundary conditions the SCF calculations consider only an octant of the stalk morphology.Although two non-connected bilayers present a one-dimensional structure, in order to have the same numerical accuracy as in the stalk case, the SCF calculations for this morphology are also performed in 3D.In this case, the reecting boundary conditions allow the consideration of a single monolayer.
In the main panel of Fig. 1 we show a 3D stalk morphology aer the SCF scheme has converged.It is described by a contour plot of the normalized total density of the amphiphiles, r(r) ¼ (r A (r) + r B (r))/ r A .To clarify the internal stalk structure one quarter has been cut out.The inset of Fig. 1 presents the 3D pattern of the elds used for starting the iterative scheme.The difference is striking, demonstrating that the structure is strongly optimized by the SCF calculation so that the nal morphology corresponds to a local free-energy minimum.Fig. 1 reveals that at the interstitial positions in the hydrophobic bilayer, above and below the stalk, the density drops slightly.The amphiphiles must signicantly stretch to ll the volume at these positions and in a compressible model the drop of the density relieves some conformational frustration.
The molecular stretching and tilting as well as the strength of conformational uctuations is illustrated in Fig. 2 depicting a contour plot of the distribution in space of the hydrophobic end of a molecule, provided that the junction between the A and the B block is located close to the stalk at the position indicated by the black dot on the A-B interface.To highlight the signicant change of the amphiphile conformations (i.e., strong tilt and stretching) near the stalk we also show the contour plot of the distribution of the hydrophobic end of a molecule closer to the "unperturbed" part of the membrane (in this case the location of the junction is marked by a red dot).In contrast to SCF calculations, in continuum models such changes of lipid conformations must be incorporated phenomenologically and are essential for reducing the predicted free-energy cost of a stalk between two bilayers. 21hen comparing the thermodynamic stability of the stalk with that of two apposed, disconnected bilayers, the chemical potential of the amphiphiles, m, in both structures should be the same.To this end, we realize the above morphologies in the nVT ensemble with different numbers of amphiphiles.An example for amphiphiles with asymmetry f ¼ 0.875 is presented in Fig. 3 where the superscripts s and b denote the morphology of the stalk and of the two apposed bilayers, respectively.The upper panel of Fig. 4 demonstrates the dependence of DU of the stalk on the tension of the membrane, for amphiphiles with asymmetry f ¼ 0.875.In this case, the stalk is metastable, being characterized by a positive excess grandcanonical potential, DU z 21k B T at S ¼ 0. The SCF calculations predict that DU monotonously increases with tension (or equivalently, decreases with m).Simultaneously, they conrm the strong dependence of the thermodynamic stability of the stalk on the molecular asymmetry.This is demonstrated by the bottom panel of Fig. 4, showing the dependence of DU on S for amphiphilic asymmetry f ¼ 0.921.For this high asymmetry the stalk becomes actually thermodynamically more preferable compared to two non-connected bilayers, i.e., DU z À12k B T at vanishing S. Therefore, in principle, a spontaneous stalk formation from two non-connected bilayers should be favored.Of course this argument does not account for potential freeenergy barriers along the transition pathway and we do not compare the free energy of the stalk with that of the invertedhexagonal phase, 13 which will compete for stability at such high molecular asymmetries.
The trend in the change of DU with S can be related to the dependence of the excess number of molecules in the stalk Dn(m) on m.In particular, from eqn (8) one obtains vDU/vm ¼ ÀDn(m).Since in all cases considered here DU decreases with m it follows that the stalk morphology has more amphiphiles than the two bilayers at the same tension.This conclusion is less obvious than it might seem at a rst glance, since the difference in the number of molecules between the two morphologies is not determined only by the stalk itself but also by the long-range perturbation of the bilayer it induces.This is demonstrated in the main plot of Fig. 5 showing the area density of the amphiphiles, n(r), as a function of radial distance   It can be seen that indeed the depletion region reduces signicantly the total excess and the effect is more pronounced in the case of the negative tension where the total excess is almost 50% smaller than at the peak of the prole.Interestingly, within the range of the parameters considered here and in agreement with experiments, 12 the shape of the stalk appears to be highly universal when rescaled by the thickness of the bilayer, L b , away from the stalk, i.e., it is independent of tension, S, and molecular architecture, f.This conservation of the stalk morphology is demonstrated by rescaled plots in Fig. 6, showing a side-view of the iso-surfaces r(r) ¼ 0.1 (cf.Fig. 1) calculated for the stalk structure for two amphiphilic asymmetries (f ¼ 0.875 and f ¼ 0.921) and different tensions.The height of the stalk has been measured with respect to the center plane of the bilayer in the "unperturbed" region.

Conclusion
The results above favorably compare with the simulation results of a very similar, particle-based model 20 using thermodynamic integration.The SCF calculations, however, are several orders of magnitude less computationally demanding.This point is particularly important because of the large parameter space that characterizes membrane systems (e.g., molecular architecture, membrane tension, distance between the apposing membranes).Advancements of numerical algorithms for solving the SCF equations 27,[41][42][43] allow the systematic study of large structures in three dimensions with high numerical precision.
For amphiphiles with asymmetry f ¼ 0.875 the simulations predict the excess grand-canonical potential of the stalk in the tensionless state to be DU z 16 AE 4k B T (as apposed to DU z 21k B T in SCF calculations).In agreement with simulations, the SCF theory predicts that the geometric shape is rather universal but that the stability of the stalk is signicantly affected by the molecular architecture, i.e., DU is reduced by making the amphiphiles more asymmetric.This observation is also in agreement with experiments, [12][13][14][15][16] where a lattice of stable stalks (i.e., the rhombohedral phase) with DU < 0 is observed for smaller head groups or strong dehydration (i.e., smaller bilayer distances).In the limit of extreme dehydration, the remaining water between the bilayers is tightly associated with the lipid head groups and, in a coarse-grained solvent-free model, increasing dehydration corresponds to smaller values of f.
The rather marginal differences between the particle-based simulation and the SCF theory regarding the numerical values of the free energies can be explained by the nite interaction range and discretized chain contour used in the former as well as by uctuations of the collective densities, included in the simulations but neglected by the mean-eld approximation.It is actually remarkable that omitting uctuations does not qualitatively change the predictions, which is very inspiring regarding the usefulness of the SCF theory in describing morphologies and distribution of conformations in amphiphilic membranes.
Computer simulation and SCF calculations also agree on how membrane tension affects the stability of the stalk, i.e., the stalk becomes thermodynamically more costly as we increase the membrane tension, S. The results presented here concerning the magnitude of DU and its dependence on the asymmetry of the amphiphiles are in qualitative agreement with earlier two-dimensional SCF calculations 19,29 using a different model combining an incompressible Flory-Huggins density functional with an explicit representation of solvent by hydrophilic homopolymers.This nding corroborates that the magnitude of the excess free energy of the stalk is not strongly affected by the details of the model.However, other aspects of the thermodynamic behavior are more sensitive.For example, DU in the latter model is either almost independent of S or decreases with it. 19By considering the excess number of molecules as a function of distance from the center of the stalk, the present SCF calculations highlight that longer-range perturbations of the membrane induced by "local" structures (i.e., the stalk) can make an important contribution to thermodynamic properties.
, showing the dependence of the free energy per molecule, F/n, on the area density, n h n/A (A ¼ L y L z being the membrane area), for the stalk (open squares) and the two bilayers (open circles).The plot illustrates the high accuracy that is required from SCF calculations: the free-energy difference between the stalk and the two disconnected bilayers is of the order of 10 À3 k B T per molecule.The dependence of the free energy on n is tted with a quadratic polynomial in n, and the chemical potential m ¼ vF/vn and the lateral tension of the membrane S ¼ vF/vA can be calculated.The minimum of the Helmholtz free energy per

Fig. 1 3 .
Fig. 1 Main panel: distribution of the normalized density of the amphiphiles, r(r) ¼ (r A (r) + r B (r))/ r A in a single stalk morphology, where r A ¼ 40 amphiphiles per R e 3 .The complete stalk structure is obtained from the reflecting boundary conditions by mirroring the octant of the stalk in which the optimization was performed (red box).Inset (bottom-right): 3D pattern of the W A (r) field used to start the iterative solution of the SCF formalism.For clarity only the hydrophobic region is shown.

Fig. 2
Fig. 2 An octant of the stalk morphology is shown, where the color contour plot corresponds to the spatial distribution of the hydrophobic end of an amphiphile (red contours correspond to highest probability) provided that the location of the junction of the A and the B blocks is at the position marked by the thick dot and the arrow.Two cases are considered where the junction is located (a) near the stalk (thick black dot and arrow) and (b) away from it (thick red dot and arrow).The iso-surfaces r(r) ¼ 0.5 and r(r) ¼ 1 of the stalk are also shown in gray.

Fig. 3
Fig. 3 Helmholtz free energy per molecule as a function of the total area density of the amphiphiles, n, for the stalk (open squares) and the two bilayers (open circles).

Fig. 4
Fig. 4 Excess grand-canonical potential DU of the stalk as a function of the membrane tension, S, for molecular asymmetries f ¼ 0.875 (top panel) and f ¼ 0.921 (bottom panel).The right axis shows the dependence of the chemical potential, m, on tension (solid line).

Fig. 5
Fig. 5 The main graph shows the area density of the amphiphiles, n(r), as a function of radial distance, r, from the center of the stalk for two different tensions, S, (open and solid squares).n(r) of bilayers at the considered tensions are indicated by open and solid circles.The inset demonstrates the change of the excess number of molecules, Dn, as a function of r.

Fig. 6
Fig. 6 Iso-surfaces r(r) ¼ 0.1 calculated for the stalk for two amphiphilic asymmetries f ¼ 0.875 (circles) and f ¼ 0.921 (squares) at different tensions.For each pair, f and S, the iso-surfaces have been rescaled by the corresponding bilayer thickness, L b , shown at the legend in units of R e .Tensions are expressed in units of k B T/R e 2 .