Nonequilibrium patterns and shape fluctuations in reactive membranes

A simple kinetic model of a two-component deformable and reactive bilayer is presented. The two differently shaped components are interconverted by a nonequilibrium reaction, and a phenomenological coupling between local composition and curvature is proposed. When the two components are not miscible, linear stability analysis predicts, and numerical simulations show, the formation of stationary nonequilibrium composition/curvature patterns whose typical size is determined by the reactive process. For miscible components, a linearization of the dynamic equations is performed in order to evaluate the correlation function for shape fluctuations from which the behavior of these systems in micropipet aspiration experiments can be predicted.


I. INTRODUCTION
Polar lipids self-assemble and orient, with the hydrophilic portions facing water. The water may be sandwiched by two lipid layers, as in the black spots in soap bubbles that develop at locations where the two layers are closer than the shortest wavelength of the visible spectrum and that occur just before the bubble bursts, or the water may be outside of the two lipid layers, as in biomembranes. While it is known that the structure of a biological membrane is far more complex than a simple lipid bilayer [1,2] because of embedded proteins, cholesterol molecules, and ions, to name just a few of the components that provide the full functionality to this structure that is crucial to the life of the cell, it is nevertheless acknowledged that the lipid bilayer is the basic structural unit of all cell membranes. An understanding of this simpler system is therefore extremely important in an effort to shed light on the properties and functioning of active transport, signaling, and adhesion in cells. Furthermore, artificial lipid bilayers are widely used in a number of nanotechnological applications ranging from solar energy transduction and biosensors to drug development.
A lipid bilayer is highly flexible and liquid-like (as is a real membrane). It can therefore not be viewed as a static inert boundary but must be recognized as a dynamical structure [3]. In the case of a biological membrane, lipid bilayers serve as quasi-two-dimensional solvents for proteins and all other components, and are intimately involved in many biochemical processes. Moreover, its ability to change its shape (curvature) is an essential property of the lipid bilayer since the formation of vesicles and the permeability properties of the cell depend on it.
From the modeling viewpoint, it became recognized during the 90's that some internal degrees of freedom are necessary to understand the large variety of conformational changes found in cell and synthetic membranes. In particular, the local composition of the bilayer can crucially affect its local curvature, and some models were developed on the basis of this idea [4,5,6,7,8,9]. However, these approaches considered membranes as equilibrium systems. In a biological context, this hypothesis is at best hopeful.
Nonequilibrium conditions are ubiquitous in nature, and for this reason the out-of-equilibrium behavior of membranes, both in cellular systems and in laboratory-prepared vesicles, has attracted much attention over the past few years. The first successful approach to the study of nonequilibrium membranes was introduced by J. Prost et al. [10]. Roughly summarizing their modeling approach, they suggested that some externally activated components (intramembrane proteins) act as "pumps," generating forces on the membrane that locally change its curvature. Variations, improvements and sequels of the initial model [11,12,13,14] as well as related experimental studies [12,15] to a large extent complete the understanding of the nonequilibrium behavior of bilayers with inserted active components.
Our interest in this work lies in a different nonequilibrium that may arise from externally activated chemical processes involving the transformation of the elementary membrane lipid components. As an example, nice experiments by P. G. Petrov et al. [16] show the effect of an ongoing photocontroled chemical reaction on the curvature of synthetic giant vesicles. Other evidence of chemically-induced shape transformations is found in the nervous synaptic process: a reaction that interconverts two differently shaped constituent phospholipids of the membrane plays a crucial role in the fission of vesicles in nervous cells [17,18]. Recent experiments also point to the importance of high-curvature lipids in many cell processes such as membrane fusion [19]. So far, no models for this kind of nonequilibrium situation can be found in the literature. Our aim is to gain physical insight into the role of a generic nonequilibrium reaction acting on a membrane that is, in turn, described through a curvature/composition coupling. We show how the reactive process endows the membrane with a rich pattern formation phenomenology and specific characteristics of its shape fluctuations.
We present a simple kinetic model of a two-component deformable bilayer with a chemical reaction. We model nonactive membranes where the nonequilibrium nature of the system is caused by a chemical reaction that interconverts two differently shaped components of the bilayer rather than by inserted active components. In the model, the membrane is composed of two species: lipid A, which is assumed to be coned-shaped, and lipid B, which has an inverted cone shape (see Fig. 1). We consider the simplest scenario where the outer layer of the membrane is composed of A and B lipids, whereas the inner layer is composed of a single component without any curvature effect. We also prescribe that lipids do not move between the inner and outer layers. Both of these simplifications are made in most membrane models. According to these assumptions, the membrane is simply modeled as a laterally heterogeneous elastic surface with an internal composition order parameter locally coupled to the curvature. Within this framework, we are basically interested in two situations. On the one hand, for immiscible components, we show that phase separation in the membrane leads to the spontaneous development of structures involving heterogeneous distributions of both composition and curvature that finally result in stationary finite-sized nonequilibrium domains. This instability is studied by performing a linear stability analysis of the kinetic equations, and the patternselection role of the reactive process is established. Corresponding numerical simulations will be shown to support these predictions. On the other hand, for miscible components we derive the correlation function for the membrane shape fluctuations. The fluctuation spectrum is mainly determined by the bending energy at small surface tension, and here we find the same expression for nonreactive and reactive cases, but with different bending rigidities. This change in the behavior of the fluctuation spectrum between equilibrium and nonequilibrium states might be easily controlled and observed in micropipet aspiration experiments. This paper is organized as follows. In Sec. II the study of membranes of immiscible components is addressed. We propose a free energy functional, derive the kinetic equations, and perform the linear stability analysis of these equations. Numerical simulations are carried out and some representative results are shown. In Sec. III bilayers of miscible components are analyzed and the change of the height fluctuation spectrum between the reactive and nonreactive situations is established. We conclude with a brief summary in Sec. IV.

A. Model and analytical results
The membrane is defined as a two dimensional surface with a concentration difference order parameter φ and a local extrinsic curvature H. The rigidity of the membrane leads to an elastic energy contribution κ 2 (H − H sp (φ)) 2 dxdy to the total energy, where κ is the bending rigidity modulus, and the spontaneous (equilibrium) curvature H sp (φ) reflects the shape asymmetry between the two lipid components. For simplicity, we adopt a linear dependence on φ, H sp = φH 0 , with H 0 > 0 according to the schematic in Fig. 1. In the Monge parametrization [20] a deformable surface is described by (x, y, h(x, y)), where h(x, y) is the displacement (height) field for the local separation from the flat conformation. This representation is valid for surfaces that are nearly flat with only gradual variations of h, and allows the approximation H ≈ ∇ 2 h. As a function of these variables, the proposed energy functional reads where the first three terms correspond to the typical Ginzburg-Landau expansion responsible for phase separation (α, β, γ > 0), with an equilibrium concentration difference φ eq = ± α/β, and a typical interface length ζ = γ/α. For self-assembled free membranes, the surface tension contribution ( σ 2 |∇h| 2 ) in the free energy can be neglected, and we have not included it in Eq. (1).
The kinetics of φ follows a conserved scheme [21] plus the reaction contributions, where Γ = k + + k − and φ 0 = (k − − k + )/(k + + k − ). k + and k − are the forward and backward reaction rate constants, respectively. Considering a permeable membrane (i.e., ignoring hydrodynamic interactions) we adopt the following relaxational equation for the evolution of the height field, where Λ is a mobility parameter proportional to the inverse of the typical relaxation time τ h . The kinetic equations can readily be adimensionalized: energy is measured in units of k B T , time in units of τ h , and length in units of √ Dτ h . In terms of the new dimensionless parameters, the kinetic equations become At thermal equilibrium Eqs. (4) describe the spinodal decomposition of two immiscible components. Due to the composition/curvature coupling H 0 , both fields will form complementary patterns. As phase segregation progresses, membrane regions with positive (negative) φ deform in such a way that the curvature become positive (negative), as shown in Fig. 2. In the absence of reaction, this coarsening process does not end until there is complete segregation into two large domains. A nonequilibrium reaction such as the one we propose converts one species into the other. This amounts to a large-scale mixing mechanism that counteracts the short-scale ordering effect of phase separation. Therefore, the segregated structures grow only until mixing and ordering effects compensate, resulting in a stationary pattern. These types of nonequilibrium patterns are also found in other systems such as polymer blends [22,23] as well as in monomolecular adsorbtion on metal surfaces [24,25], and have to be distinguished from typical Turing patterns [26]. Even though they emerge from the same kind of instablity (see below), Turing patterns are expected in systems with species with different diffusivities. In the model presented here, the domains result from the competition between a local thermodynamic affinity of equal species and a nonequilibrium reaction mixing effect. As we will see, linear stability analysis and numerical simulations support this idea. The stationary uniform state corresponds to φ = φ 0 and arbitrary h. The linear stability of these uniform solutions is tested by adding small plane-wave perturbations of wave numer q and linearizing Eqs. (4). This procedure determines the 2 × 2 linearization matrix L, with the following coefficients, The eigenvalues ω q of the Jacobian associated with the matrix L correspond to the linear growth rates of the perturbations. Solving the eigenvalue problem we obtain . At the instability boundary, Re(ω q ) vanishes for one finite wave number that is defined as the first unstable mode. If the imaginary part of ω q is not zero at this wave number, we have a wave bifurcation. The condition for this bifurcation is obtained by requiring T r[L] = 0 and ∆[L] < 0. For positive κ, Γ, α, β and γ, these conditions do not apply at any real wave number, and consequently, wave instability is not found in this model. On the other hand, if the imaginary part of the growth rate is zero at the bifurcation point, we have a Turing-like bifurcation. The condition for this bifurcation is Det[L] = 0, whose analytical expression is easily obtained, yielding the following condition for the model parameters, and for the wave vector of the first unstable mode, A phase diagram is shown in Fig. 3. In the unstable region, Turing-like patterns emerging from the competition between phase separation and reaction are predicted. In the stable phase, although the two components are not inherently miscible, the reaction completely mixes the system, and the bilayer becomes stable and essentially flat.
In the unstable phase, Γ < Γ c , there is a range of unstable modes. Comparing with the equilibrium case (Γ = 0) where a continuous range of modes starting at q = 0 is unstable, reaction stabilizes the long wavelength modes so that, as we anticipate, the long time distribution of the system results in a stationary nonequilibrium pattern with a finite size (see further discussion in Sec. II B). Notice in Fig. 4 that increasing Γ reduces the range of unstable modes progressively until Γ reaches the marginal value in Eq. (6), above which there is no longer any instability. The limits of the range of unstable modes are given by which are independent of curvature parameters. Curvature reduces the unstable mode growth rates (see line with circles in Fig. 4) but without changing either q ± or the marginal condition (6). The effect of curvature is exclusively related to the kinetics of the phase separation process (see below). Before numerically solving the kinetic equation, we can seed light on the expected pattern formation process by performing a weakly nonlinear analysis using the amplitude equations technique. This analysis allows us to compute a solution for our problem near the bifurcation threshold, find the universality class of the pattern formation mechanism, and explain some properties such as the spatial arrangement of the patterns found in the numerical simulations.
For simplicity, we restrict the calculation of the amplitude equations to one spatial dimension. In spite of this restriction, we will be able to predict the kind of spatial arrangement found in two dimensions by using the universality properties of the amplitude equations. Details of the derivation of the amplitude equations are presented in the appendix. The analysis reveals that a solution for Eqs. (4) near the bifurcation to pattern formation reads where c.c. stands for the complex conjugate, and the amplitudes A and B satisfy the equations Moreover, since B is purely relaxational, we can adiabatically eliminate it and reduce (9) to a single evolution equation for the amplitude A, Therefore, the amplitudes satisfy the real Ginzburg-Landau equation. If Γ > Γ c then the amplitudes relax to zero and a homogeneous state (φ = φ 0 and arbitrary h = h) is obtained. However, if Γ < Γ c then the amplitudes reach a stationary value and a pattern develops if q c ∈ R. Notice that the amplitudes present a negative-positive aspect, that is, at sites where A reaches a maximum (minimum) B reaches a minimum (maximum). This is simply a consequence of the particular selection of sign for the spontaneous curvature H 0 . With regard to the nonlinear term, A |A| 2 , note first that its coefficient is always negative. As a result, the nonlinearity always plays a stabilizing role, the bifurcation to pattern formation is supercritical for all values of the parameters, and no hysteresis may occur. Secondly, the nonlinear term provides information about the relevant modal interactions. At this point, we can take advantage of the universality properties of the amplitude equations to infer the spatial arrangement of the pattern in two dimensions. Since the Swift-Hohenberg model also shares the same universality class when performing an amplitude equation analysis, namely, it also reduces to the real Ginzburg-Landau equation [28,29], the pattern formation mechanism is the same in both models. It is well-known in that case that the modal interactions in two dimensions are such that if the inversion symmetry is preserved then roll-like patterns develop. On the other hand, if that symmetry is not fulfilled, a hexagonal structure appears. Note that in Eqs. (4) the inversion symmetry, Thus we expect a roll-like pattern in that case, while hexagons will develop otherwise.

B. Numerical results
Numerical integration of Eqs. (4) has been performed in two dimensions using an explicit Euler scheme in a square lattice with periodic boundary condition. Small random perturbations around φ = φ 0 and h = 0 are implemented as initial conditions. The coordinate step ∆x was chosen equal to 1, and the time step was usually ∆t = 10 −4 to assure good numerical accuracy (length and time in dimensionless simulation units). The numerical results presented in this section correspond to highly immiscible components (deep quench, α = 1 and β = 1, leading to an equilibrium value of φ eq = ±1), and an interface thickness of the order of the space discretization (ζ = 1, which leads to γ = 1). The bending rigidity modulus is taken equal to 10 (in units of k B T ) [27], and we consider two constituent lipids of very different shapes by setting H 0 = 0.2. All our numerical results are consistent with the predictions of the linear stability analysis and the amplitude equations. Representative numerical results are presented below.
In Fig. 5 we show the simulation results for three different situations. The first row corresponds to small Γ = 0.05 and a critical quench (φ 0 = 0), showing the development of a laberynthine pattern that is still evolving at t = 2000 (the longest time for the snapshots shown in the figure). The coarsening process, however, is progressively slowed down later on. The φ-field profiles of these domains (horizontal cross-sectional cuts through the pattern) reveal regions where φ eq is +1 or −1 connected by abrupt boundaries that indicate the short spatial range over which the equilibrium order parameter changes from one value to the other. This is a signature of the fact that the phase separation process is dominant for this value of Γ. The second and third sets of panels correspond to large Γ = 0.2 (the marginal condition for the selected parameters is Γ c = 0.25) showing laberynthine (critical quench in the second row) and quasi-hexatic droplet-like (off critical quench, φ 0 = 0.1, in the third set) patterns. A spatial Fourier transform of this stationary pattern shows a clear hexagonal structure. These domains are already stationary at times longer than t = 2000, and can indeed be considered as nonequilibrium stable phases of the system, involving both composition and curvature modulations. Their φ-and h-field profiles have a smooth harmonic shape due to the strength of the reactive process. In both cases, since we are close to the bifurcation boundary, small local deviations from the stationary values are obtained. We especially monitor the variations of the height field and find that |∇h| < 0.05 is satisfied, so that the model has a real physical correspondence.
In order to assess the kinetic ordering process more quantitatively, we monitor the domain size L(t) computed from the composition correlation functions φ(r ′ , t)φ(r ′ + r, t) (the same results are obtained using the height correlation functions). Here, the brackets · · · indicate not only an average over orientations of r and over surface positions r ′ , but also over different realizations of random initial conditions. As has been reported in other studies, curvature considerably slows down the coarsening segregation process. This is specially evident when reaction is absent, and for this situation different kinetic approaches have been invoked. A mean-field kinetic scheme by Taniguchi [7] leads to extremely slow laberynthine stripe growth obeying L(t) ∼ t r with r = 0.1, instead of the usual spinodal decomposition growth exponent of 1/3 [30]. Later, Monte Carlo simulations [8,9] showed how, in the long time evolution, those stripes break up into disconnected buds that subsequently diffuse and coalesce. In our model, when the reactive term is removed, similar results as those presented in Ref. [7] are obtained, and no stripe breakup is observed. The reason for the disagreement between the analytic and Monte Carlo schemes is probably due to the fact that continuum models subjected to a specific parametrization of the surface such as given in Eqs. (4) do not allow for overhangs that are surely crucial in the late stages of membrane phase segregation.
However, the cases in which there is no reaction or in which the reaction is weak are not of interest for us, since in these cases the system evolves in such away that large gradients of the height displacement field occur. As noted above, when this happens the Monge surface parametrization is not valid and the model, although mathematically robust, does not describe the physical behavior of any real system. In the reactive cases, however, the slow kinetics is still observed for the times prior to stationary pattern formation. This is observed in Fig. 6, where a set of curves showing the domain size L(t) is presented for several values of Γ. The situations with and without curvature are compared for each case. Notice how the deformable systems evolve more slowly than the nondeformable ones, although the same final stationary size is achieved. There is no theory to explain such a slowing down effect, but some hand-waving arguments may explain the physical reasons for this behavior. Monitoring of the different energy contributions in Eq. (1) indicates that at early times (when the kinetics with and without curvature are still quite similar) interfacial energy due to the rapid formation of composition domains is much larger than any other energy contribution. In other words, phase separation proceeds and membrane curvature follows the composition change. At this stage the reduction of the interfacial energy governs the coarsening process, leading to the 1/3 Lifshitz-Slyozov exponent. Later, at intermediate times, when composition domains become larger, the height order parameter is no longer able to keep up with the phase separation process and some curvature energy is stored, causing the whole coarsening process to slow down. The final size of the nonequilibrium stationary domains, L f = L(t → ∞), is determined by the reaction parameter. The dependence of the final pattern size on the reaction parameter, L f ∼ Γ −s , has been largely discussed in the literature [22,31,32]. The results of our model also reproduce the two limiting behaviors: s = 1/4 for large Γ (close to its marginal value), and s = 1/3 for small Γ. In Fig. 7, L f is plotted for different values of the reaction parameter. Linear fits for the four first and the last four data points give the slopes 0.323 ± 0.006 and 0.26 ± 0.006, respectively. The derivation of the exponent values for the rigid situation (nondeformable surfaces) is performed by minimizing an effective free energy expression (where the reactive term is included via Green's functions) in a square (small Γ) and a harmonic (large Γ) approximation for the stationary concentration field [31,32]. In the present model, the curvature kinetics relaxationally follows the concentration dynamics. At longer times (when the stationary patterns are already achieved) no significant curvature energy is present in the system, so that we can neglect the curvature contributions in Eq. (1) and the results in Refs. [31,32] are recovered.
In Eq. (3) we ignored the hydrodynamic effects due to the background fluid velocity by considering a permeable membrane through which the fluid drains freely as the system evolves. Solvent hydrodynamic effects [33] can be approached through a renormalized height mobility Λ = (4µq) −1 in the second of Eqs. (4) in Fourier space, where µ is the solvent kinematic viscosity. The inclusion of this q-dependent mobility does not change the main results of the stability analysis and numerical simulations shown so far.
Before ending this section, we check the applicability of our results to real membrane systems. We adimensionalized the kinetic equations to units in which D = Λ = k B T = 1. Redimensionalizing them, taking the typical values D = 10 −7 − 10 −8 cm 2 /s (for lipids in a liquid-phase membrane), µ = 1cp (H 2 O at 20 • C) and the typical unstable modes q ≈ q c = 0.5, we find that the size of the patterns obtained above lies within the range 2 − 20µm, which is accessible in giant vesicle (10 − 100µm in diameter) experiments.

III. MEMBRANES OF MISCIBLE COMPONENTS: LINEARIZATION AND SHAPE FLUCTUATIONS
In this section we consider a parameter region where the membrane is thermodinamically stable. The instability in the previous section was caused by the immiscibility of the two lipid components, while we now consider the case of miscible membrane constituents. Notice that the two situations could correspond to the same physical system but at different temperatures, below (unstable) or above (stable) the critical temperature. For the stable case we adopt the following free energy functional, where α is now negative, and the nonlinear (β) and line tension (γ) terms have been removed since they are irrelevant in the absence of phase segregation. An additional surface tension term has been included in order to study membranes under tension. In the equilibrium situation (no reaction), the concentration field can be integrated out ( ∂F ∂φ = 0) from Eq. (11), leading to φ eq = κH0 κH 2 0 −α ∇ 2 h. Inserting φ eq in Eq. (11) we obtain an effective free energy for h alone, with an effective rigidity modulus [34] Notice that, as a consequence of the composition/curvature coupling, the different preferred curvature of the lipid components acts to reduce the rigidity modulus from κ to κ ef f . Note also that κ ef f > 0 for negative values of α, so that the stability of the membrane is assured above the critical temperature. The dynamics for the immiscibility situation has so far been considered to be strictly deterministic. Here, however, in order to obtain dynamical steady state functions, we add stochastic forces to the kinetic equations and perform an average over an appropriate ensemble. The dimensionless kinetic equation for φ, once the reactive terms are included, reads where the last term corresponds to a conserving Gaussian noise with correlations f φ (r, t)f φ (r ′ , t ′ ) = 2k B T δ(r − r ′ )δ(t − t ′ ). For the curvature order parameter we again adopt the relaxational equation for a permeable surface, with Λ = 1 in dimensionless units, and f h is a thermal equilibrium noise with correlations The kinetic equations become The important quantity to characterize membrane shape fluctuations is the height variance |h q | 2 at wave number q, which is calculated by linearizing the kinetic equations (16) and solving them in Fourier space. Ifφ(q, ω) and h(q, ω) are the Fourier transforms of φ(r, t) and h(r, t), respectively, Eqs. (16) can be written as where Solving these coupled equations, and using the statistical properties of the thermal noises, yields ĥ (q, ω)ĥ * (q, ω) = 2k B T a 2 21 q 2 + ω 2 + a 2 11 ω 2 (a 11 + a 22 ) 2 + (a 12 a 21 − a 11 a 22 + ω 2 ) 2 .
On integrating over ω, the height variance follows the general expression, where ν and η are negative variables that read and In the absence of the reaction, the evaluation of Eq. (20) in the long-wavelength limit for a tensionless membrane leads to |h q | 2 = k B T /κ ef f q 4 . This behavior is truncated and replaced by k B T /σq 2 for a membrane under tension if only the dominant terms at small q are retained. Keeping the dominant and the first subdominant terms, one recovers the well known expression for the height variance in an equilibrium membrane, Notice that this result is obtained in a much simpler way from an equilibrium average using the effective free energy in Eq. (12). One of the usual experiments to study equilibrium membrane shape fluctuations is based on the micropipet aspiration technique [35]. In these experiments, a pressure difference is applied inside a micropipet in contact with a vesicle membrane. This creates a tension in the membrane that pulls the excess area ∆S due to the thermal shape fluctuations inside the micropipet. By means of this technique, the areal strainᾱ ≡ ∆S/S is obtained experimentally for different values of the applied tension σ. According to Eq. (23), the relative areal strain ∆ᾱ ≡ᾱ −ᾱ 0 (ᾱ 0 being the areal strain for a reference value σ 0 ) can be calculated analytically [36,37]: Therefore, the slope of the logarithm of σ versus ∆ᾱ obtained by the micropipet technique yields a measure of the effective bending modulus. Note that in these experiments a certain tension σ is needed, but it has to be small since otherwise the κ ef f q 4 term in Eq. (23) would be insignificant compared to σq 2 . Now the question is, how is the shape fluctuation spectrum affected by the presence of the reaction? In other words, how would the nonequilibrium membranes described here behave in micropipet experiments? To answer this, we evaluate Eq. (20) for Γ = 0 and keep only the dominant terms in the limit q → 0. For tensionless membranes we get |h q | 2 = k B T /κq 4 and for membranes under tension we recover |h q | 2 = k B T /σq 2 . Keeping the dominant terms for both limits, in the weak tension regime the fluctuation spectrum reads Comparing this result with the fluctuation spectrum obtained for membranes with active proteins [12,13,15], no novel nonequilibrium contribution to the height fluctuations is found here. The effect of the nonequilibrium reaction in stable membranes results in a change between a regime governed by κ ef f to another one with the actual rigidity κ. Thus, the reaction makes the membrane more rigid by simply removing the composition/curvature coupling effect that diminished the rigidity in the equilibrium situation. This change is more evident when the two components of the bilayer are rather different in shape and miscible but close to the critical temperature. In this situation (large H 0 andᾱ 0) the difference between κ and κ ef f might be experimentally observable. Some caution, however, must be observed when deriving Eq. (25). In order to evaluate Eq. (20) for small q, we considered the terms Γ 2 to be dominant with respect to the terms 2Γ κH 2 0 − α q 2 (subdominant) and also with respect to the subsubdominant quartic contributions σ + κH 2 0 − α q 4 (which are those that led to Eq. (23)). However, the analysis becomes more difficult if one notices that we are not in the region of asymptotically small wavenumbers since the smallest q accessible to an experiment is of order 1/L, where L is the linear system size. Thus, one must look in detail at the values of the parameters in order to determine exactly which terms initially considered subdominant might indeed be dominant when the reaction is present. Γ 2 is dominant if Γ > 2 κH 2 0 − α q 2 min , with q min ∼ 1/L and thus dependent on the system size. We know that for a sufficiently strong reactive process, Γ 2 dominates and Eq. (25) , and therefore the parameter A ′ in Eq. (22) has to be divided by (4µq) as well. However, with these modifications the evaluation of |h q | 2 again leads to the results in Eqs. (23) and (25).

IV. CONCLUSIONS
Starting with a simple model of a deformable reactive membrane composed of two differently shaped molecules, we show that stationary finite-sized patterns may appear under some parameter conditions for the immiscibility situation as a result of the competition between phase segregation and reaction. These structures involve heterogeneous distributions of composition and curvature whose sizes are determined by the nonequilibrium reactive process. For typical values of the viscosity of water and lipid lateral diffusion constants in bilayers, and at normal room temperatures, such patterns are predicted to have a size of a few microns (see the discussion at the end of Sec. II B). Therefore, this behavior would correspond to a reliable pattern formation mechanism in lipid membranes which we believe to be experimentally accessible in giant synthetic vesicles. The amplitude of these patterns is modulated by the bilayer rigidity and the spontaneous curvature of its components. In our numerical calculations we have used realistic typical values for the rigidity, while the spontaneous curvature depends on the specific geometry of the membrane constituents. We specifically propose that azobenzene compounds, which are known to show amphiphilic behavior in Langmuir monolayers, and whose shapes are strongly modified by means of well-known photoisomerization reactions [38,39], might be suitable to test our predictions. The selection of the applied light wavelength and its intensity may determine both the fraction of the two isomers (φ 0 ) and the strength of the reactive process (Γ), respectively. Thus, experimental work on synthetic vesicle membranes made of these compounds can be specifically designed to confirm the results of this model. In the same context, for miscible components membranes we have found a difference between the equilibrium and nonequilibrium situations. The effects of the composition/curvature coupling and the reactive process on the membrane rigidity are established. Micropipet experiments with the proposed membrane systems might confirm these results.
For convenience, we have introduced in Eqs. (A1) some notation changes and definitions with respect to the ones that appear in Eqs. (4). Thus, in Eqs. (A1) subscripts indicate partial derivatives, the field ϕ = φ − φ 0 has been defined, and we have introduced the control parameter ε = (Γ c − Γ) /Γ c that accounts for the "distance" to the bifurcation between a homogenous state (ε < 0) and pattern formation (ε > 0). The homogeneous state according to this definitions corresponds to ϕ = 0 and arbitrary h = h.
As shown in Sec. (II A), by linearizing Eqs. (A1) one can easily check that if ε > 0 then the homogeneous state becomes unstable and, is a solution in the steady state if q 2 c = α − 3βφ 2 0 / (2γ). It is worth noting that in Eqs. (A2) we have arbitrarily taken h = 0 without any loss of generality. Moreover, by substituting (A2) into (A1) and expanding up to the first harmonic, O (exp (iq c x)), one finds that the amplitudes scale as a function of ε as A, B ∼ √ ε. Thus, we expect that near the bifurcation the following expansion holds, By computing the linear growth rate, ω q , that is, the largest eigenvalue of the linear problem (see Eqs. (5)), we also note that, where C i are constants. Thus, as a function of ε the width of the band of unstable modes scales as ∼ ε 1/2 . Then, since all modes exp (iqx) can be written as exp (i (q − q c ) x) exp (−iq c x), a separation of spatial scales can be performed between the most unstable mode (fast) and the rest of the modes of the unstable band (slow). Let us call the slow modulation spatial scale X, such that X = ε 1/2 x, where x will now stand for the fast spatial scale. We also note that exp (ω q t)|q→qc ε→0 ≃ exp (εt) .
We could continue up to any order with the expansion. In all cases we will obtain a nonlinear equation, such that at order ε n/2 , Lχ n = ψ n ϕ (1) , . . . , ϕ (n−1) ; h (1) , . . . , h (n−1) . (A5) However, at order ε 3/2 we are already able to extract a closed evolution equation for the amplitudes of the pattern and so we will stop at that order. Our task is to solve the hierarchy of equations given by (A5). At order ε 1/2 the problem is homogeneous and with appropriate boundary conditions, is a solution. However, the amplitudes A and B are undetermined at this point. The subsequent orders are no longer homogeneous and therefore their solvability can not be ensured unless one implements the so-called Fredholm alternative theorem [28]. In our case the application of the theorem simply states, as a recipe, that for Eqs. (A5) to have a solution the functions ψ n can not contain the fundamental mode exp (±iq c x). Thus, by substituting the solution (A6) into the next order or the hierarchy and imposing the solvability condition we obtain Once again, the value of A and B can not be determined at this order. However, at order ε 3/2 the application of the Fredholm theorem provides the conditions that determine the values of the amplitudes A, B. These conditions constitute the amplitude equations for our pattern forming system,