Generation of dynamic structures in nonequilibrium reactive bilayers

We present a nonequlibrium approach for the study of a flexible bilayer whose two components induce distinct curvatures. In turn, the two components are interconverted by an externally promoted reaction. Phase separation of the two species in the surface results in the growth of domains characterized by different local composition and curvature modulations. This domain growth is limited by the effective mixing due to the interconversion reaction, leading to a finite characteristic domain size. In addition to these effects, first introduced in our earlier work [Phys. Rev. E {\bf 71}, 051906 (2005)], the important new feature is the assumption that the reactive process actively affects the local curvature of the bilayer. Specifically, we suggest that a force energetically activated by external sources causes a modification of the shape of the membrane at the reaction site. Our results show the appearance of a rich and robust dynamical phenomenology that includes the generation of traveling and/or oscillatory patterns. Linear stability analysis, amplitude equations and numerical simulations of the model kinetic equations confirm the occurrence of these spatiotemporal behaviors in nonequilibrium reactive bilayers.


I. INTRODUCTION
Due to hydrophobic repulsions, amphiphilic molecules such as lipids spontaneously aggregate in water to form bilayer membranes. These bilayers typically exhibit an in-plane fluidlike nature, and are highly flexible surfaces so they can display a large variety of conformations. The experimental study of equilibrium lipid bilayers has attracted a great deal of attention in the past decades [1,2]. Issues such as membrane conformational behavior, shape fluctuations, fusion and fission, and phase segregation in multicomponent bilayers have been extensively investigated [1,2,3]. Among many other manipulation techniques, the fabrication of synthetic giant vesicles and planar lipid membranes [4] as well as micropipet aspiration [5] have been used to perform these studies, leading to a fairly detailed current understanding of equilibrium membranes.
Parallel to the experimental advances, theoretical modeling of equilibrium flexible membranes has also come a long way. Typical theoretical approaches are based on the search for minimum energy conformations according to the Canham-Helfrich bending elasticity description [6,7] of a tensionless membrane, plus (optionally) the corresponding area-difference contributions for closed surfaces [8]. More recently, it has been recognized that internal degrees of freedom such as chemical composition can crucially influence the shape of the bilayer, especially when dealing with multicomponent membranes undergoing phase separation. In particular, models containing a local coupling between composition and curvature have been developed [9,10,11,12,13,14], result-ing in an interesting outcome, namely, the appearance of phase separating domains that grow and adopt distinct membrane curvatures. This is understood as the initiation mechanism for budding phenomena in real membranes. Kinetic mesoscopic schemes, Monte Carlo approaches, and a variety of discrete dynamic algorithms have been proposed on this basis, and have mainly been devoted to the study of the equilibrium behavior and the dynamics toward equilibrium of flexible bilayers.
The ultimate motivation for the study of lipid bilayers lies in their likeness to cell membranes. A biological membrane is a complex mixture of lipids, sterols and proteins that represents the main structural component of the cell architecture. Rather than an inert static boundary, a membrane has to be viewed as a dynamical surface directly and actively involved in many biological cell processes [15]. In addition to their compositional complexity, cell membranes are also continuously subjected to nonequilibrium driving forces, chemical gradients, and energy fluxes, so it should be recognized that nonequilibrium behavior may underlie many aspects of their dynamics [16,17,18]. Therefore, although experimental and theoretical studies on thermally equilibrated membranes have had considerable success and are a good starting point, nonequilibrium approaches are fundamental for a more accurate understanding of the dynamical properties of both natural membranes and synthetic lipid bilayers.
Excellent work on nonequilibrium lipid bilayers has recently been presented by the experimental groups of J. Prost and H.-G. Döbereiner. J. Prost et al. have studied the effects of active proteins inserted in a vesi-cle membrane. Acting as ion pumps externally activated by light, these proteins provide a nonthermal energy source that directly affects membrane shape behavior [19]. As an example of another nonequilibrium source, H.-G. Döbereiner et al. investigated the morphological transformations of a vesicle membrane due to a photoinduced chemical reaction that modifies the local lipid composition of the bilayer [20]. Some progress has also been made on the modeling front. For example, parallel to their experimental findings, a novel nonequilibrium scheme has been proposed by J. Prost et al. that successfully explains the observed experimental phenomenology [21,22]. In addition to this proposal, other nonequilibrium situations have been modeled very recently, such as those of membranes confined between parallel plates [23], membranes with active inclusions [24], and bilayers near repulsive walls [25].
We are mainly interested in the nonequilibrium situation proposed by H.-G. Döbereiner et al., where the bilayer's own lipid constituents are chemically transformed, leading to changes in the membrane curvature. The influence of lipid composition on the curvature of cell membranes has recently been established in cell fusion processes where high-curvature lipids play an important role [26]. Chemically induced lipid modifications are also known to be responsible for membrane shape transformations involved in nervous synaptic processes such as the formation of microvesicles that release the neurotransmitter to the gap between two nerve cells [27]. Furthermore, in addition to biological membranes, an understanding of the coupling between chemical reactions and the interfacial curvature is essential in elucidating dynamical processes in artificial bilayers that may be useful in nanotechnology applications. As an example, selfassembly techniques based on biomineralization can lead to hierarchically organized materials built from reactive amphiphilic interfaces [28].
Based on the nonequilibrium ingredients described above, in a previous paper [29] we presented a model for a two-component reactive membrane that exhibits a particular compositional and morphological organization. In that model the two lipid constituents have a distinctly different shape, so that they are able to produce distinct curvatures in the membrane. Moreover, they are assumed to be thermodynamically immiscible, so they spontaneously induce the development of phase separating domains. Additionally, an externally induced reaction (promoted by a nonequilibrium source such as a chemical flux, an applied light, etc.) interconverts the two lipids. The combination of these ingredients leads to stationary patterns involving heterogeneous modulations of composition and curvature. The novelty of the resulting pattern organization is that although stationary, the patterns are generated in a nonequilibrium context, that is, they are actively maintained by the competition between thermodynamic (equilibrium) phase separation and the environmentally induced (nonequilibrium) reaction. At thermal equilibrium, the low affinity between the two immiscible components generates domains with a characteristic composition and curvature. In the absence of reaction, these domains grow indefinitely and coalesce until complete segregation into two large domains is achieved. However, the reactive process converts one lipid component into the other, resulting in a large-scale mixing mechanism that halts the growth of the segregated structures when the mixing action balances the short-scale ordering effect of phase separation. This balance leads to a stationary state with patterns of a finite size.
In our previous model [29], we assumed that the nonthermal energy contribution that promoted the lipid interconversion reaction simply dissipates to the medium. In this paper we generalize our original model to the case where such an energetic process directly alters the membrane shape by exerting local forces on it, that is, the induced reactive process locally 'kicks' the membrane. Conformational changes of the membrane constituents due to the reaction can reasonably be expected to have mechanical effects on the membrane. We model this effect in a very generic way, and will demonstrate that the inclusion of this new ingredient is essential for the generation of dynamic spatial organization. As a typical feature of soft-matter systems, robust spatiotemporal phenomena such as those displayed in this paper are interesting in their applicability to real bilayer membranes.
Although our ideas are inspired by phenomena observed in real cellular systems, we hasten to disclaim a direct analogy between our model and in vivo systems which are far more complex than our stripped model. We only wish to suggest that the mechanisms proposed herein might play a role in real membranes in which the energy sources for the local membrane-shape-altering processes might include ATP consumption, the asborption of light, or any of a number of other possible energy providing mechanisms. We point out here and again later in the paper that the spatiotemporal structures of interest require parameters that constrain the range of applicability of the model. In particular, we will see that elastic properties (bending rigidity) play a central role, and that the most likely candidates for the observation of our predictions due to their flexibility are in vitro lipid bilayers formed from single-and/or short-chain lipid surfactants, especially those composed of surfactants of rather distinct shapes that are interconverted by a reactive process that can be controlled experimentally. We discuss these and other possibilities for experimental realization later in the paper.
The outline of this paper is as follows. Section II presents our basic mesoscopic free energy description as well as the derivation of the kinetic equations for the composition and curvature variables in the proposed model. The linear stability analysis and the amplitude equation analysis are presented in Sec. III, with the details of the derivation of the latter relegated to the Appendix. Organized spatiotemporal behaviors are predicted in some regions of parameter space provided that the parame-ter describing the activity of the reactive process on the membrane dynamics is sufficiently large. In good agreement with the analytical predictions, numerical simulations exhibiting traveling and/or oscillating domains are shown in Sec. IV. Finally, in Sec. V we summarize the main conclusions of this paper.

II. THE MODEL
In our earlier work we considered the simplest scenario where one layer of the membrane (say the outer one) was composed of two differently shaped lipids A and B, whereas the other layer was composed of a single component without any curvature effect. We also disregarded any lipid flip-flop exchange between the inner and outer layers. Under these assumptions, we modeled the properties of the non-active membrane as a two-dimensional surface characterized by the local concentration difference φ between the A and B components and the local extrinsic curvature H of the surface. At this point we do not specifically identify the components A and B in any further detail since this is to be seen as a schematic and highly streamlined representation of any number of possible realizations. Indeed, A and B need not even be single compounds; for example, A may be a group of lipids or biomolecules some of which induce a particular curvature, while B may be another group that induces a different curvature. A and B could be the isomers of an amphiphillic azobenzene derivative, or groups of lipids as in Fig. 2 in [27] or those suggested in [26], or simply a given lipid with its polar head more or less charged as in the experiments in [20].
The free energy functional in terms of the two order parameters was assumed to be The first three terms correspond to the typical Ginzburg-Landau expansion that leads to phase separation when α, β, and γ are all positive, with an equilibrium concentration difference φ eq = ± α/β and a typical interface length ζ = γ/α. When α is negative the components are completely miscible, leading to a homogeneous state. The last term is the elastic energy contribution due to the rigidity of the membrane, and κ is the bending rigidity modulus. For simplicity we have assumed that the spontaneous (equilibrium) curvature H sp (φ), which reflects the shape asymmetry between the two lipid components, is linear in the order parameter, H sp = φH 0 .
In the Monge parametrization [30] 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 that we have incorporated in the free energy functional. For selfassembled 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 φ and h are assumed to be driven by the free energy functional, and by the chemical reaction that interconverts the two lipid components, A ⇄ B. The effect of the reaction is two-fold: it affects the concentration difference order parameter φ, and it also affects the shape of the membrane. This last effect is our new contribution in this paper. The kinetics is thus a generalized version of that introduced in [29]: The concentration difference φ follows a conserved scheme [31] with diffusion coefficient D, augmented by the reaction contributions. The rate parameter Γ = k + + k − and the stationary concentration difference pa- are determined by the forward and backward reaction rate constants k + and k − . Note that we have assumed local overdamped relaxational dynamics for the membrane height, that is, Rouse-like dynamics, where Λ denotes a mobility parameter proportional to the inverse of the typical relaxatio time τ h . This approximation is only valid when the membrane is immersed in a high viscosity medium and/or when the membrane is highly permeable [32]. When hydrodynamic effects can not be neglected, a more realistic, albeit complex, description of the membrane dynamics (Zimm-like dynamics) is required, where inertial effects and a space dependent relaxation parameter are in order [33]. The last term in Eq. (2b) is based on the following hypothesis. The reaction that interconverts A and B can be understood as an isomerization-like chemical transformation, involving a strong modification in the shape of the membrane constituents. This process implies the displacement of parts of these molecules that could have a mechanical effect on the local membrane shape. If additionally the process is strongly energetically activated by an external energy source, it might exert a local force on the membrane. A simplifying approximation is to assume that the forward and backward reactions exert opposite forces on the membrane. These forces are assumed to act locally for a negligible period of time (the time needed to complete the reaction is much shorter than any other time scale of the system), in the same direction as the preferred curvature of the reaction product component, that is, positive (outwards) for A → B, and negative (inwards) for B → A, see Fig. 1. This is modeled in a generic way [34] by the last term in the membrane height kinetic equation, where ξ is a parameter accounting for the strength of the effect of the reaction on the shape of the membrane. Here we take the parameter ξ to have the same sign as H 0 (in Fig. 1 both are positive). Active pro-teins are known to act as force centers when inserted in lipid bilayers [19], and other experimental studies [17,18] also provide evidence that reactive processes may locally modify the membrane shape in red blood cells. However, we must recognize ours as an experimentally verifiable conjecture of what might happen at the lipid component level since the cited experimental evidence only corresponds to chemical processes involving much larger molecules (proteins).
In the next two sections we show analytically and confirm numerically that the new force term can have profound effects on the instability and pattern formation properties of the model. In particular, this term may lead to dynamical patterns, whereas in its absence only stationary patterns are possible.

III. ANALYTICAL TREATMENTS
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 to the uniform state and linearizing Eqs. (3a) and (3b) in these perturbations. This procedure leads to 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.
The solutions of the eigenvalue problem are given by . 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. If the imaginary part of the growth rate is zero, one has a Turing-like bifurcation. The condition for this bifurcation is Det[L] = 0.
In the absence of the interconversion reaction it is of course well known that immiscibility leads to complete phase separation. This occurs because a continuous range of modes starting from q = 0 is unstable, and phase separation does not stop until there is complete segregation into two large domains. In our previous work [29] we showed that the interconversion reaction provides a mixing mechanism that stabilizes the longest wavelength modes. If the reaction rate parameter Γ is sufficiently large, then mixing is complete, the instability disappears, and the bilayer becomes stable and essentially flat. If Γ < Γ 1 , the reaction no longer stabilizes all modes, but it does stabilize the longest wave vector modes. The first mode to become unstable is The longest wavelength unstable mode now lies at the finite value This means that the phase separation process stops when the separated regimes are of characteristic size ∼ q −1 min . The patterns are Turing-like (stationary) because the imaginary part of the growth rate is zero at the bifurcation point. A typical phase diagram found earlier for this case is shown in Fig. 2. Note that the pattern size is independent of curvature parameters. Curvature reduces the unstable mode growth rates (see the ξ = 0 curves in Fig. 3), but without changing either the characteristic size of the patterns nor the marginal condition (5). The effect of curvature is thus restricted to the kinetics of the phase separation process, as has been extensively analyzed in [29]. The new questions of interest here are the consequences of the reaction-shape coupling in the kinetic equations. To place the analytic results obtained in this section in context, we anticipate the numerical results of the next section and exhibit two different views of a typical phase diagram in Fig. 4. Regime III is the stable regime where the mixing due to the interconversion is sufficiently strong to produce a homogeneous phase. Here "sufficiently strong" depends on the value of the parameter ξ. Beyond this stable regime, there are now two instability regimes. One, called region II in the diagrams, is the Turing instability regime leading to stationary nonequilibrium patterns, that is, the same sorts of patterns observed in the absence of the reaction-shape coupling (upper region of Fig. 2) [29]. The condition for the Turing-like bifurcation, Det[L] = 0, obtained from the linear stability conditions now is Γ < Γ 2 , with This is the curve bounding region II in the figure. The first Turing-like unstable mode occurs at The new regime, region I, is one in which one observes dynamical (i.e., time-dependent) patterns. Provided the reaction parameter is sufficiently small and the reactionshape coupling parameter is sufficiently large, the linear stability analysis in this regime leads to wave bifurcation solutions, T r[L] = 0 and ∆[L] < 0. The explicit conditions for these time dependent solutions obtained from the linear stability conditions are Γ < Γ 0 , with and ξ > ξ 0 , with The first unstable mode in the above expression now is The curves bounding region I in Fig. 4 arise from these expressions. A typical dispersion relation function ω(q) in this regime is shown in Fig. 3. The function now has a nonzero imaginary part as well as a real part, as appropriate for spatiotemporal pattern formation. Typical spatiotemporal field configurations in the unstable regions of the phase diagram can only be found numerically. We do so in the next section. Before doing so, however, we can go a step further in the analytic approach to the problem by introducing a weakly nonlinear analysis based on the amplitude equation formalism. Such a nonlinear expansion is valid near the bifurcation threshold to pattern formation and sheds light on the expected pattern configurations and their dynamical properties. For simplicity, we restrict the calculation of the amplitude equations to one spatial dimension. Nonetheless, we can readily deduce the possible arrangements in two dimensions by an appropriate interpretion of the coupling terms. The details of the derivation of the amplitude equations are presented in the Appendix. The analysis reveals that a solution for equations (3a) and (3b) near the bifurcation to pattern formation reads where c.c. stands for the complex conjugate, ω 0 = |Im(ω q0 )|, and the amplitudes A ± (x, t) and B ± (x, t) satisfy the equations where .

(15b)
Note that no pattern develops if Γ > Γ 0 since the amplitudes then decay with time, A ± , B ± → 0. Note also the coupling term between the rightward, A − , and leftward, A + , traveling waves. This interaction is responsible for generating either standing or traveling waves in the membrane dynamics. Standing waves appear as a result of a kink -antikink dynamics where both waves have the same amplitude, A − = A + . Traveling waves require these amplitudes to be different. Therefore, if the coefficient c 2 of this interaction term is small (strictly speaking, zero) we expect the two amplitudes to be equal, so that that is, the solution is a standing wave. On the other hand, if c 2 is not zero, a traveling wave appears. Observe that if φ 0 = 0, that is, if the composition of the membrane is 50% lipid A and 50% lipid B, then the A + ↔ A − interaction term is always relevant and no standing waves are possible. This scenario is confirmed by the numerical simulations described in the next section (cf. Figs. 5 and 6, where we see a compositioncurvature traveling wave in the membrane. For the parameters in that figure, c 2 ≈ 0.25). On the other hand, this restriction does not apply if the membrane composition is such that φ 0 = 0 (cf. Figs. 7 and 8, where φ 0 = 0.14 and c 2 ≈ 10 −2 ; the simulations show that in those cases there is a standing oscillatory contribution to the pattern). As for the different spatial configurations, one obtains either rolls or hexagons depending on the symmetries of the system. Thus, note that in Eqs. (3a) and (3b) the inversion symmetry is satisfied only if φ 0 = 0 and in that case roll-like arrangements are obtained. On the other hand, if φ 0 = 0 there is no inversion symmetry, and the resulting structures are expected to be hexagonal.

IV. NUMERICAL RESULTS
Representative numerical results are presented in this section, showing good agreement with the predictions of the linear stability and amplitude equation analyses. We numerically solve Eqs. (3a) and (3b) in a twodimensional square lattice of 100 × 100 sites with a mesh size ∆x = 1 and periodic boundary conditions. The spatial derivatives are calculated by means of a simple centered scheme, and a first-order Euler algorithm with time step ∆t = 10 −4 is used for the temporal integration. These coordinate and time steps insure good numerical accuracy. Simulations are started from a slightly randomly perturbed homogeneous distribution φ( r) = φ = φ 0 and h( r) = h = 0. In this section we present numerical results that 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 0.5 (in units of k B T ), and very different shapes are implemented for the two constituents by setting H 0 = 0.2.
When dealing with mesosopic or coarse-grained modeling schemes one has to be concerned about the applicability of the results. The limit on the predictive power of these models is mainly determined by the identification of the correct time, length and energy scales accessible to the experiments. In our case, we have adimensionalized the kinetic equations in order to have D = Λ = k B T = 1, so we can return dimensions to these equations by choosing typical values for these parameters. The diffusion coeficient is known to be in the range 10 −7 −10 −8 cm 2 /s (for lipids in a liquid-disordered phase membrane [15]). In the absence of the nonlinear and mechanical effects considered in this paper, solvent hydrodynamic effects (Zimm dynamics) [3] can be incorporated with a renormalized height mobility Λ = (4µq) −1 in Eq. (2b) in Fourier space [32]. Here µ stands for the solvent kinematic viscosity, which is 1cp for water at 20 • C. For κ = 0.5k B T the linear stability analysis shows typical unstable modes at q ≈ q 0 ≈ 0.5, which corresponds to a pattern size in the range of 2 − 20µm. Such a size may be accessible in experiments on giant vesicles and planar bilayers. However, when mechanical effects are included, a more elaborate analysis such as that provided by dynamical renormalization methods [32] is in order. A simple and plausible hypothesis might be to include such effects by simply renormalizing both Λ and ξ in Eq. (2b), both with a 1/q decay in Fourier space since one might expect a non-local viscous kernel to decay as a power law with distance [33]. In any case, a more detailed consideration of this issue is beyond the scope of this work.
However, in order to accurately delimit the range of applicability of our model, it is useful to add some comments regarding the value of κ. The elastic properties, and thus the bending rigidity, of lipid bilayers are strongly determined by the size, shape and molecular elasticity of their constituents. Membranes composed of phospholipids (which have two hydrophobic chains) are characterized by a bending rigidity of the order of tens of k B T [35]. Cell membranes containing large molar fractions of cholesterol (rather rigid molecules) display even higher rigidities [16]. In our model, the typical sizes of the predicted dynamical patterns emerging from a wave instability for such large values of κ exceed the experimentally accessible sizes. Moreover, for κ 10, wave-like unstable modes are observed to be supressed at long times by the other existing Turing-like unstable modes, so that only stationary structures are asymptotically generated in this parameter region. However, single-and/or shortchain lipid surfactants are known to form much more flexible bilayers, with bending rigidities of the order of k B T or even less [36]. Furthermore, recent theoretical models [37] show how bilayers composed of surfactants of rather distinct shapes (as in our case) may exhibit smaller rigidities than one-component membranes. This, therefore, delimits the range of applicability of the model results in this paper.
Our quantitative assigment of ξ has been made with the energetic feasibility of a number of external sources in mind. To see this, note that the characteristic time for a single reactive event is of order 1/Γ. During this time interval the change in the height of the membrane due to the contribution of the reactive term is of order Moreover, the local energetic "cost" of such a height increase due to the reactive term is of order 2κξ 2 /(∆x) 2 [see Eq. (1)]. In most of our simulations we have used ξ = 5 (in simulation length units), which, for κ = 0.5, corresponds to an energy cost of order 25 (in units of k B T ). The dephosphorylation of an ATP releases about 50kJ/mol or 10 −19 J/ATP, which at room temperature is precisely of order 25k B T . If the energy source is light, the power required to produce 25k B T within a time interval of order Γ −1 for Γ = 0.14 as used in our simulations is only approximately 0.3 × 10 −6 mW/cm 2 , which is much lower than the power produced by typical commercial light sources used, for example, in photosensitive Langmuir monolayer experiments. While some of the energy pro-  Fig. 4 and the following set of parameters: α = β = γ = 1, φ0 = 0, κ = 0.5, H0 = 0.2, Γ = 0.14 and ξ = 5 (dashed curve in Fig. 3). The two upper panels correspond to the φ-and h-field distributions at t = 15000, once the structures are robust. The stripes travel from left to right at constant velocity. In the third panel we plot the temporal evolution (horizontal axis) of a one-dimensional cross-section φ-profile (vertical axis) at y = 50, starting at t = 15000 until t = 19000. In the fourth panel, the corresponding h-profile (horizontal axis) is plotted against time (downward below sheet axis), starting at t = 15000 until t = 15500. Darker (lighter) regions are richer in the A (B) lipid in the concentration snapshots, and some exaggeration along the vertical direction has been applied to the height plots.
vided by these external sources (ATP, light, etc.) would be used for processes other than the purely elastic motion of the membrane, including dissipation into the thermal surroundings and the movement of other masses, this order of magnitude estimate shows that the mechanism is feasible and provides a basis for the value of the parameter ξ chosen for our simulations.
For a highly flexible bilayer, typical phase diagrams are those of Fig. 4. As noted earlier, region I corresponds to a regime where at least one of the unstable modes has nonzero imaginary part, so both types of instablity are present. Region II is the unstable phase also present in Fig. 2. Here, the unstable modes have no imaginary parts, so that stationary finite-sized patterns are predicted. In region III, there are no unstable modes (stable phase). Here we have exhibit spatiotemporal patterns associated with the new region I.
The nature of the spatiotemporal patterns is determined by the value of φ 0 and by the associated relative magnitudes of the amplitudes of waves traveling in different directions. These assertions are supported by the amplitude equation analysis of the previous section. Thus, for a critical quench (φ 0 = 0), the system typically displays some transients with domains that travel in differ-ent directions until they organize into a coherent train of traveling stripes. The post-transient spatiotemporal behavior is described in the different panels in Fig. 5. We plot the concentration and height fields (upper panels), and one-dimensional cross-sections showing the temporal evolution, for both order parameters. The traveling waves are consistent with the predictions of the amplitude equation analysis. Numerical profiles at a given post-transient stage reveal the mechanism that leads to the motion of the generated spatial structures. In Fig. 6 we observe that φ-and h-profiles are slightly displaced, the field φ being ahead in the direction of propagation. That there is such a phase difference between the two fields is consistent with the appearance of the contribution with an imaginary coefficient on the right hand side of Eq. (14b).
Off-critical quenches (φ 0 = −0.14) display different spatial and temporal behaviors, again consistent with the amplitude equation results. As before, we first observe transients, but now involving concentration droplets and bud-like surface deformations. The fields settle into an oscillating pattern of buds rich in the minority species and whose geometry depends on the value of the parameter ξ that measures the mechanical effect of the reaction on the shape of the membrane. This behavior is shown in Fig. 7. For ξ = 5 a spatial Fourier transform of the pattern indicates clear square symmetry. For ξ = 3, the domains are rather hexagonal and they oscillate as well as move in one direction, as shown in Fig. 8.

V. CONCLUSIONS
In a previous model for a flexible reactive bilayer composed of two differently shaped lipids, we demonstrated the occurrence of stationary finite-sized domains of composition and curvature as a result of the competition between thermodynamic phase segregation and a nonequilibrium reaction [29]. A generalization of the model that includes the local effect of the reaction on the membrane shape has been presented here. We have shown how the mechanical influence of the reactive process on the membrane may lead to the formation of spatiotemporal structures. The linear stability analysis of the model equations shows the existence of a wave instability for sufficiently large reaction-shape coupling. A weakly nonlinear analysis based on the amplitude equation formalism provides insight into the spatial and temporal symmetries of the emerging patterns. Correspondingly, numerical simulations display the spontaneous generation of traveling lamelar phases, and oscillating or moving buds, showing an important potential for the formation of spatiotemporal patterns in deformable nonequilibrium reactive lipid membranes.
This extension of the model should be considered as a formal issue in reactive deformable surfaces. However, we believe that experimental work on giant vesicles and on planar membranes can be designed to qualitatively capture the spatiotemporal phenomenology obtained in our model. Azobenzene derivatives, which are also known to display dynamical organization in nonequilibrium Langmuir monolayers [38], may be suitable compounds to test our predictions since their shapes can be strongly modified by means of well-known (externally controlled) photoisomerization reactions [39], and the resulting isomers normally exhibit phase separation.
Finally, we note that it would be interesting to study the crossover between Rouse and Zimm dynamics in our model system so as to elucidate how the different parameters are renormalized due to solvent hydrodynamic effects, and to clarify the resulting consequences for the pattern formation mechanism proposed herein. Such a study has been carried out for tethered membranes, where different coarsening exponents are found [32,33]. For convenience of notation, we rewrite the one dimensional version of the model equations (3a) and (3b) for reactive membranes in terms of the shifted field variable ϕ = φ − φ 0 and the control parameter ε = (Γ 0 − Γ)/Γ 0 that accounts for the "distance" to the threshold of the bifurcation to pattern formation so that a pattern develops if ε > 0: The subscripts indicate partial derivatives. The homogeneous state is ϕ = 0 and arbitrary h =h. With no loss of generality we takeh = 0.
Slightly above the bifurcation to pattern formation, the following expansion holds (see [29,41] and references therein): and a separation of spatial scales can be implemented between the most unstable mode (fast growth), q 0 , and the rest of the modes within the unstable band (slower growth). In terms of the control parameter, we let X = ε 1/2 x denote the spatial modulation scale of the slow modes and T = εt the associated time scale. The separation of scales between the fastest growing mode and the slower modes can be implemented in Eqs. (A1a) and (A1b) by expanding the spatial and temporal derivatives according to the chain rule so that ∂ x → ∂ x +ε 1/2 ∂ X and ∂ t → ∂ t + ε∂ T . Implementation of this scale separation and substitution of Eqs. (A2a) and (A2b) into (A1a) (A1b) leads to a rather cumbersome expansion in terms of ε. The lowest order contribution is of order ε 1/2 , and balancing all terms of this order gives Eqs. (A3a) and (A3b) can trivially be written as The next order of the expansion gives us terms of O(ε) and yields At order ε 3/2 we get Lχ 3 = ψ 3 ϕ (1) , ϕ (2) ; h (1) , h (2) , FIG. 8: Spatiotemporal behavior for the same case in Fig. 5 except for φ0 = −0.14 and ξ = 3 (still in region I). In the upper panels, the droplet-like structures move from the upper-left corner to the lower-right one at constant velocity, while they oscillate. The snapshots correspond to t = 18000 and the temporal profile evolutions go from this time up to t = 22000 for φ, and to t = 18500 for h. The same panel organization and specifications as in Fig. 5 apply.