Reaction-driven assembly: controlling changes in membrane topology by reaction cycles

Abstract

Chemical reaction cycles are prototypical examples how to drive systems out of equilibrium and introduce novel, life-like properties into soft-matter systems. We report simulations of amphiphilic molecules in aqueous solution. The molecule's head group is permanently hydrophilic, whereas the reaction cycle switches the molecule's tail from hydrophilic (precursor) to hydrophobic (amphiphile) and vice versa. The reaction cycle leads to an arrest in coalescence and results in uniform vesicle sizes that can be controlled by the reaction rate. Using a continuum description and particle-based simulation, we study the scaling of the vesicle size with the reaction rate. The chemically active vesicles are inflated by precursor, imparting tension onto the membrane and, for specific parameters, stabilize pores.

---

1 Introduction

Amphiphilic molecules are commonly used to artificially produce compartmentalized nanostructures through their ability to form bilayers and vesicles in aqueous solution, may it be lipids,^1,2 diblock copolymers^3–6 or even short peptides.^7,8 Such assemblies have promising applications in therapeutics^9–11 and as models for synthetic cells.^12–18 In living systems, the size of compartments is tightly controlled and compartments undergo topological changes in the course of transport processes, e.g., in membrane vesicle trafficking.

In order to mimic life-like properties, such as e.g., spontaneous adaptation, one can employ beyond-equilibrium dynamics by incorporating mechanisms that allow the system to dissipate energy.¹⁹ One way to do so is to exploit reaction cycles, where a high-energy molecule, a fuel, enables a forward reaction to a product that is (partly) immiscible with the reactant but, after some time, deactivates back into the reactant state.^20,21 Already the simplest of such systems, an immiscible reactant-product system exhibits a variety of effects that differ from immiscible binary mixtures such as an arrest of droplet coalescence^22–24 or a shape instability for emerging droplets, leading to fission.^25–27

It has remained a challenge to reproduce these predicted effects experimentally for which high reaction rates are necessary.

For instance, Heckel and coworkers observed the arrest in coarsening in the course of spinodal decomposition of reacting polymers.²⁸ Tena-Solsona et al., however, observed an accelerated coalescence of a chemically fueled mixture compared to its equilibrium counterpart.²⁹ More complex reaction cycles have shown to produce exciting new features in recent years, such as the transient assembly of fibers,^30,31 colloids and gels^32–38 or polymeric micelles that were used as nano-reactors.^39,40 The transient assembly of vesicles was achieved by the use of an ATP-hydrolizing reaction cycle,⁴¹ peptides that are switched to their cyclic anhydride at the expense of a carbodiimide⁴² or by polymerization induced self-assembly with a pH or light-responsive deactivation.^43,44

In this work, we elucidate which effects may occur for a class of molecular architectures, where the product is amphiphilic. Specifically, we present simulations of macromolecules in aqueous solution that are involved in a reaction cycle, switching between a hydrophilic and an amphiphilic state. In the latter state, the amphiphiles self-assemble into bilayer membranes or vesicles. We explore the effects of this reaction-driven assembly (RDA) by employing two complementary simulation schemes, a continuum model and a particle-based description. We study the stationary states for varying reaction rates, both analytically and by simulations, and describe the kinetic pathways that are taken to form vesicles. Furthermore, we show how compartments formed by RDA become inflated with precursor, which, in turn, stabilizes the formation of pores in their membrane.

2 Methodology

We consider a system that is comprised of an aqueous solvent, S, small fuel molecules, F and amphiphiles, modeled as asymmetric diblock copolymers. The diblock copolymers can, in principle, also represent different molecules, such as short peptides, to which we expect our results also qualitatively apply. The two blocks of the copolymer are denoted A and B. The minority block A is always hydrophilic. In the presence of fuel, the majority block B may react to become hydrophobic, forming an amphiphilic product P that self-assembles to micelles or vesicles. The reaction cycle is completed by a spontaneous reversion of the tail B from hydrophobic to hydrophilic. In this state we denote the effective homopolymer as R (reactant). Hence, there are two types of chemical reactions: (i) a binary forward reaction in which the precursor molecule reacts with the fuel to form the amphiphile, P, and (ii) the spontaneous deactivation of product. Within both models the average fuel concentration is fixed which corresponds to a situation where fuel is refilled continuously and diffusing quickly. A sketch of this reaction cycle is depicted in Fig. 1.

Fig. 1 Sketch of the reaction cycle: a hydrophilic molecule acts as precursor, reacting with fuel to form an amphiphile. In aqueous solution this product self-assembles, forming micelles or vesicles. The product can spontaneously deactivate to become hydrophilic, again.

Note that the reaction simultaneously switches the hydrophilicity of the entire tail block. Such a reaction is feasible experimentally, for instance, with a small molecule where a single reaction at the tail end can change the physical properties of the complete molecule.

2.1 Particle-based simulations

The particle-based simulations employ a soft, coarse-grained model that represents several monomer repeat units by a single particle. Within the soft, coarse-grained model the Hamiltonian is split into strong bonded (b) and weak non-bonded (nb) interactions, . The strong bonded interactions are taken to be harmonic springswhere m runs over all molecules, b over all bonds within the molecule, and r_m,b, r_m,b+1 refer to the positions of the bonded particles. Further, k_B is the Boltzmann constant, T the temperature, N₀ = N_P = N_R = 20 refers to the chain-contour discretization of the polymer, and R_e to its root mean-squared mean end-to-end distance. f_AN₀ = 6 particles of a copolymer are of type A. Solvent and fuel molecules are represented by dimers, N_S = N_F = 2.

The weak non-bonded interactions are expressed in terms of the normalized densities ϕ_c(r) of component c at position r, i.e.,

where κ₀ = 6 characterizes the inverse compressibility of the system, χ_cc′ are the Flory–Huggins parameters describing the binary repulsion between the components c and c′. Solvent and fuel are chemically identical. The minority component of the amphiphile repels solvent (and fuel) less than the majority component, χ_AS = 1 and χ_BS = 4. The components of the amphiphile repel each other, χ_AB = 2. In the reactant state, the polymer is only comprised of A segments. High repulsion between A and S result in unwanted morphologies, e.g., onion-structures⁴⁵ whereas χ_AS = 0, slows down the exchange between aggregates via the diffusion of precursor in the solution and increases the free-energy barrier for merging of aggregates. denotes the invariant degree of polymerization and sets the strength of thermal fluctuations. n denotes the total number of beads in the cubic simulation cell of volume, V, with periodic boundary conditions. The system sizes are in the range 15 ≤ L/R_e ≤ 20. The densities are calculated on a cubic collocation grid with linear spacing Δx = R_e/8.

We employ the single-chain-in-mean-field (SCMF) algorithm^46,47 that temporarily replaces the weak, non-bonded interactions by external fields and thereby exploits the different strengths of strong bonded and weak but computationally costly non-bonded interactions. Particle positions are updated by the smart-Monte-Carlo algorithm, using the strong bonded forces to bias the trial displacement. The time it takes a copolymer to diffuse its own mean end-to-end distance, R_e, in a disordered system is τ₀ = 2616 MCS and serves as time unit. We use the highly parallel and graphics processing unit (GPU)-accelerated software SOft coarse-grained Monte carlo Acceleration (SOMA).⁴⁷

New in this work is the use of polymer-type conversions that model macromolecular reactions. We attempt the forward reaction, R → P, of every hydrophilic precursor molecule to an amphiphilic product not after every MCS but with a period τ^conv. Typically we use 10 ≤ τ^conv/MCS ≤ 50. The attempted conversion is accepted with probability where r_cm denotes the center-of-mass of the polymer. Likewise the deactivation of the product to the reactant is accepted with the concentration-independent probability, . After this sequence of conversion attempts the density fields and external fields are recomputed.

2.2 Continuum model

As a second, complementary tool, we employ a continuum model for micro- and macrophase separation in block copolymer solutions, introduced by Uneyama and Doi.⁴⁵ It operates on the same set of parameters, allowing for a direct comparison. Within the Uneyama–Doi model (UDM) local concentrations, ϕ_pi(r), for each polymer architecture p and its respective blocks i at position r are order-parameter fields, and the free energy is given as a functional of the former. Chain conformations are not treated explicitly at the advantage of significantly shorter simulation times and the possibility of analytical predictions. Conceptually, the UDM is an extension of Ginzburg–Landau theories for multicomponent systems and remains valid even beyond the weak-segregation limit (WSL) by accurately approximating the long-range interactions.

The molecules are comprised of N_p segments and the block fractions are denoted by f_pi. N_S = N_F = 1. The ratios N_P/N_S = N_P/N_F = 10 coincide with the particle-based model. Given the local concentrations, the free-energy functional takes the form^45,48

where A_p,ij and C_p,ij are polymer-architecture-dependent coefficients, is a long-range kernel, which fulfills the modified Poisson equation, , with cutoff length ξ_cut = 2R_e,⁴⁸ and π(r) is a Lagrange field that enforces local incompressibility.⁴⁸ In the free-energy functional the second term describes a long-range interaction due to the covalent bonding of the two blocks of the copolymer. For solvents and homopolymers, we have f_p = 1 and the coefficients are given by A_p = 0, C_p = N₀/N_p. For the diblock copolymer with blocks A and B, these take the formwith f_B = 1 − f_A = 0.7 and the definitionand s(f) as taken in the work of Ohta and Kawasaki.⁴⁹ In accord with prior work,⁵⁰ the parameter C_p,ij differs from the original version of the UDM⁴⁵ to ensure that the order–disorder transition matches the result of the random-phase approximation.⁵¹

The third and fourth terms in eqn (3) describe the entropic contributions. For a homopolymer blend, the second term reduces to the Flory–Huggins entropy of mixing,^52,53 and the third term is chosen to match the order–disorder transition for diblock copolymers. The fifth term represents a square-gradient penalty to match the structure factor at large wavevectors. The sixth term accounts for the local, binary repulsions between two distinct block types. χ_AS = 0.5, χ_BS = 4, and χ_AB = 2.

Given the free-energy functional, we obtain the chemical potentials as the functional derivative of the free energy with respect to the concentrations, , with the joint index c = (p,i). We give their explicit form in section Chemical potentials and time evolution in the UDM of the ESI.† The spatio-temporal evolution of the concentrations follows model-B dynamics⁵⁴ that locally conserves the concentrations. Gradients in the chemical potentials give rise to fluxes, , with Λ_cc′(r) being Onsager coefficients that describe the concentration-dependent mobility of the components. Ignoring nonlocality in space or time,^55,56 assuming that the molecular mobilities do not depend on local concentrations, and enforcing incompressibility via a Lagrange field, we use a diagonal Onsager matrix, . The parameter λ is related to the diffusion coefficient by λR_e² and dictates a time scale, λ⁻¹ = τ₀, which is the diffusion time of the reference polymer in the particle-based model and will be taken as the reference.

The time evolution of the locally conserved concentrations is given by a continuity equation,

∂_tϕ_c(r,t) = −∇·j_c(r,t) + ξ_c(r,t) (7)

with . ξ_c(r,t) is thermal noise, which we take to be Gaussian with its moments dictated by the fluctuation-dissipation relation⁴⁸

〈ξ_c(r,t)〉 = 0 (9)

This concludes the contributions to the time evolution of the system without reactions, which are valid close to equilibrium and with which the system will relax into the equilibrium configuration.

A second contribution to the time evolution of the concentration fields stems from the chemical reactions. Within an infinitesimal time interval, dt, the reaction probability for the precursor is given by . The spontaneous deactivation to the precursor state occurs at probability, . Reactions do not conserve the order parameters and extend the model-B dynamics, eqn (10), by the additional reaction-induced concentration change

Here the index runs over all reactions with component c as a reactant, and runs over all reactions with component c as product and are the corresponding probabilities. This yields the final kinetic equationIn the last term the ratio of reaction rates, r_b and r_f, to inverse diffusion time, λ, appears and the interplay of these two processes dictates a length scale, as demonstrated below.

In both models we apply homogeneous initial conditions, with mean concentrations for the different species c in the system volume V. Unless noted otherwise, we use ρ_R + ρ_P = 0.25, ρ_F = 0.1, and ρ_S + ρ_R + ρ_P + ρ_F = 1. Eqn (12) is integrated with the time step of Δtλ = 4 × 10⁻⁵ in a cubic volume (12.8R_e) covered by a grid with spacing Δx = R_e/10. Further details of the numerical implementation are described in section Numerical implementation of the ESI.†

3 Results and discussion

3.1 Analytical predictions

We start with analytical considerations for the initial structure formation and the stationary state, based on the continuum model, to show that emerging aggregates will arrange on a lattice and take on uniform size. Afterwards, we compare this to the simulation results.

To make progress, we introduce two simplifications: (i) we assume the fuel be homogeneously distributed throughout the system. Hence, the forward reaction of precursor (reactant) to product occurs with rate r_fρ_F and effectively becomes a first-order reaction. In the long-time limit, the mean concentrations become independent of the morphology and are given by and ρ_P + ρ_R = const., as determined by the initial conditions. (ii) Moreover, we assume the hydrophilic reactant, R, be distributed like the solvent, S, and lump the two components into one, hydrophilic species H, with ϕ_H = ϕ_R + ϕ_S and ϕ_R = ρ_Rϕ_H/ρ_H. By virtue of incompressibility, 1 − ϕ_H(r) = ϕ_A(r) + ϕ_B(r), i.e., the system locally has only two degrees of freedom, ϕ_A and ϕ_B.

For the short-time structure formation, we perform a linear stability analysis of eqn (12) around the spatially homogeneous, initial state in terms of the deviations, δϕ_A(r) = ϕ_A(r) − ρ_A and δϕ_B. Details are provided in section Growth rates of phase-separation processes of the ESI.† The linearized time evolution takes the form

with and the linearized evolution operator Σ[∇²]. This is readily solved in Fourier space by the superposition of two exponentially growing modes,

δϕ(q,t) = c₊(q)δϕ₊(q)e^σ₊(q)t + c₋(q)δϕ₋(q)e^σ₋(q)t (14)

where δϕ_±(q) denotes the eigenmodes of the evolution operator at wavevector q and σ_±(q) the corresponding eigenvalues, respectively. The coefficients, c_±(q), determine the initial configuration.

The two local degrees of freedom give rise to two distinct demixing characteristics: (i) phase separation between amphiphile and solvent (termed “PH”, ) and (ii) (micro)phase separation of the two blocks of the amphiphile (denoted “AB”, ). In general, however, these demixing characteristics are not eigenmodes of the evolution operator but linear combinations thereof. The character of an eigenmode can be quantified by the wavevector-dependent angle φ_±(q). Here, the angle corresponds to the process PH, whereas signals the microphase-separation process AB. To distinguish which demixing characteristics dominates, we obtain their initial growth rates, [small sigma, Greek, tilde]_AB/PH, by projecting onto the eigenmodes.

The results of this analysis are presented on the left of Fig. 2 for the typical parameters of this simulation study. It becomes clear that the initial growth of density fluctuations is dominated by the demixing of amphiphilic molecule from hydrophilic solution at length scale ∼R_e because the exponentially growing eigenmode σ₊(q) has an angle close to φ₊ ≈ φ* and a maximum at wave number q*R_e ≈ 2π. In turn, microphase separation between the blocks of the amphiphile will only commence after the amphiphilic product has separated from the solvent, in regions of higher amphiphile concentration. Such kinetics have also been observed in experimental systems of amphiphilic block copolymers in solution.⁵⁷

Fig. 2 Spinodal structure formation of an amphiphilic diblock copolymer in aqueous solution (N_P = 10, χ_ABN_P = 20, χ_AHN_P = 5, χ_BHN_P = 40) with reactions: top left: wavevector-dependent growth rates of eigenmodes σ_±, as well as initial rates for the two demixing characteristics, amphiphile-solvent segregation [small sigma, Greek, tilde]_PH and block–block segregation [small sigma, Greek, tilde]_AB. The inset presents the small-q behavior. Bottom left: Angle, φ_±, between the eigenmodes and demixing characteristics. φ = φ* and correspond to amphiphile–solvent and block–block segregation, respectively. Top right: Effective, inverse structure factors of amphiphile–solvent and block–block concentration fluctuations as function of wavevector q. The families of curves correspond to different reaction rates, r_b, as indicated in the key. Bottom right: Dependence of the characteristic length scales on the reaction rate for weak segregation.

Without reactions, the evolution operator in the linear spinodal regime and the static structure factor, are related by for proper choice of a non-diagonal Onsager matrix as detailed in section Effective free-energy functional of the ESI.† We use this relation to define an effective structure factor of block–block concentration fluctuations, , and likewise for . In the ESI,† we confirm that this identification is compatible with an effective free-energy functional of the reactive system, that can be mapped onto the Ohta–Kawasaki model,⁴⁹ describing a diblock copolymer melt without reactions.

The inverse of these effective static structure factors, and , are presented in the top right panel of Fig. 2 for various reaction rates, r_b. As expected, the inverse block–block structure factor, , exhibits a minimum at a finite wavevector, , that indicates the characteristic inverse length scale of microphase separation between the two blocks of the amphiphile in the WSL.⁵¹ Moreover, is rather insensitive to the reaction rate, r_b.

Intriguingly, the effective static structure factors, also exhibits a minimum at a finite wavevector, ; qualitatively similar to the equilibrium structure factor of a copolymer melt in the WSL.⁵¹ This indicates that the phase separation between product molecules and solvent does not occur macroscopically but, instead, is characterized by a finite, microscopic length scale, . This results from the interplay between diffusion time of precursor molecules through the solution and their lifetime, dictated by the reactions. Thus, the effect of reactions can be conceived as introducing a connectivity between product molecules and solvent that prevents macroscopic phase separation and induces stationary amphiphile-rich domains of a characteristic size. In analogy to the equilibrium behavior of asymmetric copolymers without reaction,⁵¹ the nonequilibrium stationary state of the reactive system at small product concentration will result in amphiphile-rich droplets that densely pack into a body centered cubic (BCC) lattice.

The top right panel of Fig. 2 demonstrates that the location of the minimum, , shifts to large wavevectors and becomes less deep upon increase of the reaction rate, r_b, i.e., larger r_b results in smaller amphiphile-rich droplets and decrease the effective incompatibility between amphiphiles and solvent.

In the ESI,† we explicitly map our reactive system onto the Ohta–Kawasaki model,⁴⁹ describing the equilibrium of a diblock copolymer melt. Using the results for the domain spacing of equilibrium diblock copolymers in the WSL, the characteristic distance between the domains scales like^58,59

which is verified on the bottom right of Fig. 2.

In the strong segregation limit (SSL), i.e., at smaller reaction rates, r_b > 0, the domain spacing is larger and the analogy to the Ohta–Kawasaki model yields^49,58–60

The above analysis yields a prediction for the dense packing of building blocks, arbitrarily shaped product aggregates, assuming that the two phase separation processes evolve separately, as observed for amphiphiles in solution.⁵⁷ In the following, the mean product density, ρ_P, is small such that the aggregates form droplets that ideally arrange on a cubic lattice. At large values of r_b, compact micelles form inside these droplets whose radius, R, scales like the lattice spacing, , of the amphiphile-rich droplet lattice. For smaller reaction rates, however, the lattice spacing becomes larger, and so does the number of amphiphilic molecules in a unit cell, . In this case, the amphiphiles form bilayers that close into unilamellar vesicles. We assume that the bilayer thickness, D, is much smaller than the vesicle radius, D/R ≪ 1, i.e., the number of molecules in a vesicle of size R scales like DR². In this case, the scaling of the radius, R, of densely packed vesicles with reaction rate takes the form

This case is relevant for our work.

When the reaction rate is decreased further, or the average product concentration is increased, the unilamellar vesicle becomes too large to fit into the unit cell, i.e., . In this case, a dense arrangement of multilamellar vesicles will form or in more extreme cases, stacked planar membranes.

3.2 Kinetics of structure formation

To test the above predictions for the emergence of dense packings of vesicles, we start our simulations from a spatially homogeneous distribution of components. We treat the fuel only implicitly, fixing the forward reaction rate of the precursor r_fρ_F to be constant without explicitly considering the fuel. We will later assess the quality of this approximation for a hydrophilic fuel.

For the initial formation of vesicles from a homogeneous solution, two pathways have been proposed: one describes the emergence of vesicles via a disc-like micelle that spontaneously bends to minimize the line tension along its edge and finally closes,^48,61–63 as shown schematically in the inset of Fig. 3(a). We refer to this kinetic mechanism as ‘pathway I’. Alternatively, a vesicle may emerge directly from a spherical micelle via a semi-vesicle, i.e., a micelle in which hydrophilic head groups are enriched in the center, by flip-flopping of amphiphiles to the inside. This way, the vesicle's bilayer emerges gradually,^64,65 see the inset of Fig. 3(b). We refer to this process as ‘pathway II’. Which of these paths is taken, depends on the nature of the amphiphiles. In nonreactive systems, pathway I is favored for a stronger repulsion between head group and tail. Because of this, pathway I appears to be the generic process in experiments.^66–69 Pathway II, in turn, is only rarely evidenced in experiments.^70,71

Fig. 3 Time evolution of the hydrophobic-tail densities in the UDM starting from a homogeneous initial configuration, both, in real space, ϕ_B(r), and in reciprocal q-space, . Reaction rates are chosen (a) r_fρ_F = 4 × 10⁻²λ, r_b = 1 × 10⁻²λ and (b) r_fρ_F = 1.6 × 10⁻¹λ, r_b = 4 × 10⁻²λ, and the time steps are selected to visualize the most relevant events. The insets on the right demonstrate the vesicle formation mechanism, (a) pathway I and (b) pathway II.

To illustrate the time-evolution of the reaction-directed assembly, we take two different reaction rates and show morphological snapshots and the corresponding structure factor of B-density fluctuations in Fig. 3. Additionally, we provide the Videos udm-impl-fuel-rb1e-2.mp4 and udm-impl-fuel-rb4e-2.mp4 in the ESI.†

The time evolution for an exemplary reaction rate is illustrated in Fig. 3(a). Starting from a homogeneous distribution at tλ = 0, small micelles quickly nucleate. Initially, their size is on the order of ∼R_e, as visible in the structure factor at tλ = 2 × 10⁻², and subsequently they quickly coalesce. Rather than forming larger spherical micelles, the structures elongate to form cylindrical micelles, tλ = 2.2 × 10¹, and eventually disk-like micelles, tλ = 5 × 10¹. The latter are unstable. They bend to minimize the line tension along their edge and finally close the remaining pore to form a vesicle. Hence, vesicle formation proceeds via pathway I.

Smaller vesicles continue to coalesce, tλ = 9 × 10¹. In the late stage, coarsening arrests and the vesicles adopt a rather uniform, finite size. Since product molecules are continuously deactivated and precursor molecules are continuously activated within the aqueous solution, the vesicles constantly exchange material, via diffusion of the hydrophilic precursor molecules through the solution. This is in marked contrast to micelles and vesicle in systems without reactions, where exchange of amphiphiles across the solution is extremely slow. This exchange allows for an efficient size matching because all vesicles are equally capable to compete for product as well as to adjust their positioning in the lattice structure. In the stationary state one obtains a structure factor that reflects the lattice periodicity at low wave numbers and is reminiscent of existing predictions for vesicle form factors at high wavenumbers.^72–74

The material exchange and size matching can be quantitatively investigated by using a Hoshen-Koppelman cluster analysis (HKCA)^75,76 to obtain the radii, R, of individual micelles and vesicles. A description is outlined in section Measurement of micelle and vesicle radii in the ESI.† The individual and mean radii as a function of time are presented in Fig. 4. As observed qualitatively in Fig. 3(a), the radius distribution of the vesicles becomes narrow and all vesicles approach a common mean radius by the redistribution of material.

Fig. 4 Temporal evolution of individual and mean aggregate sizes, as well as its variance in the UDM for the system of Fig. 3(a). 0 ≤ tλ < 2 × 10¹: initial coarsening, 2 × 10¹ ≤ tλ < 1.2 × 10²: sporadic fusion of micelles and vesicles. Ultimately, the efficient exchange between aggregates via diffusion of hydrophilic precursor through the solution, results in a well-define and narrowly distributed aggregate size. The different processes are depicted by the pictograms in the insets.

For a second, larger reaction rate, r_b, presented in Fig. 3(b), the behavior of the structure factor, [scr S, script letter S]_BB, is qualitatively similar but we observe two differences: (i) the fast reaction rate, specifically the quicker deactivation rate of product molecules inside the aggregates, leads to an increased densities of precursor material on the inside. For a growing micelle, this facilitates the formation of a semi-vesicle, i.e., a large micelle with a hydrophilic core, made up of precursor molecules and head-groups, such as structures observed at tλ = 2.2 × 10¹–4 × 10¹. By flip-flopping of amphiphiles to the inside, a small vesicle can form immediately, following pathway II. This effect is qualitatively similar to the formation of complex, multiphase condensates formed by RDA.⁷⁷ (ii) Additionally, as predicted by the analytical consideration, coarsening arrests at smaller vesicle sizes.

The corresponding results in the particle-based model are given in section Kinetics of structure formation in particle-based simulations in the ESI.†

3.3 Scaling with reaction rates

To quantify the dependence of the size scale on the reaction rates, we measure the vesicle radii for a range of reaction rates. We start with the implicit fuel.

This size dependence is depicted in Fig. 5 for both the UDM (a) and two different sizes of the simulation box for the particle-based simulations (b). As expected, larger vesicles are observed for smaller reaction rates, while for high reaction rates, the distance, , between aggregates is small and they have smaller size, R. Note that the error bars for the mean radii only show the standard deviations of individual radii, ignoring that different realizations of the structure formation may result in different number of aggregates and thus slightly different mean radii at finite simulation times. Small variations in the numbers of vesicles in the simulation cell are expected because fusion events become rare as the system approaches the stationary state as well as due to finite-size effects. To estimate the spread among realizations, we repeated the simulations for reaction rates r_b = 1 × 10⁻²λ and r_b = 4 × 10⁻²λ in the UDM simulations twice, yielding a deviation of the mean radii within the standard deviation of each distribution, as visible in Fig. 5(a).

Fig. 5 Dependence of stationary micelle and vesicle radii on the reaction rates in (a) for the UDM after time tλ = 10³ and in (b) for the particle-based simulations after t = 3.4 × 10³τ₀. The stationary morphologies for a subset of reaction rates are shown as the insets. We calculated the vesicle and micelle radii by a HKCA,^75,76 assuming a spherical aggregate shape. Radii of individual aggregates are given as semi-transparent dots, whereas the mean and standard deviation are depicted in solid colors. In (a) two additional simulations were performed for r_b = 1 × 10⁻²λ and r_b = 4 × 10⁻²λ, only indicated by the mean radii. In (b) two different sizes of the simulation box are probed, depicted in shades of red and blue, as indicated in the legends. For all three scenarios a power-law fit was performed with mean radii as fit values and standard deviations as weights.

For all three scenarios we fit a power-law R ∼ r^ν and obtain the exponent ν = −0.51 ± 0.06 for the UDM, and ν = −0.44 ± 0.07 for the particle-based simulations with the small simulation box and ν = −0.62 ± 0.08 for the larger one. Thus, within the statistical uncertainty we confirm the scaling of the vesicle size, , predicted by eqn (17) in the SSL.

3.4 Inflation of chemically fueled vesicles

Compartments formed by RDA tend to accumulate precursor on the inside. To analyze the origin of this effect, we set up a single vesicle within the particle-based model and track the evolution of individual molecules. A detailed analysis is given in Sec. Precursor enrichment in RDA vesicles of the ESI.† When amphiphilic product molecules in the inner leaflet of a vesicle deactivate to become hydrophilic precursors, they majorly diffuse to the inside of the vesicle, and vice versa. The diffusive transport of hydrophilic precursor across the vesicle membrane is protracted because the hydrophilic precursor needs to pass the hydrophobic membrane core. Thus, the number of precursors trapped inside the vesicle, V_in, is proportional to the number of molecules in the inner leaflet, and the same holds for the number of precursor molecules in the volume, V_out, between the vesicles. according to eqn (17). For the parameters studied, V_out/V_in > 1.

Hence, vesicles formed by RDA sort precursor molecules, enriching them inside the vesicles. The higher concentration of the precursor inside the vesicle increases the osmotic pressure, Δp > 0, of the inside compared to the outside and thereby inflates the vesicle. The concomitant stretching of the vesicle membrane imparts a tension, σ, onto the membrane.⁷⁸ This membrane tension and the pressure difference between the vesicle's inside and outside are related by the Young-Laplace equation

where R̃ denotes the vesicle radius to the center of the bilayer.

We use the UDM to accurately measure the pressure difference, Δp, and the resulting membrane tension, σ. We doubled the grid resolution for precision but left all other parameters unchanged. The morphology of an isolated, RDA vesicle in a periodic cell of size, is presented in Fig. 6(b). Panel (a) depicts an equilibrium vesicle without reactions but with the same ρ_P as in panel (b). The inflation expresses itself in a visible decrease of membrane thickness, D, compared to the chemically inactive counterpart.

Fig. 6 3D density of hydrophobic tails without (a) and with (b, left) reactions (r_b = 1.5 × 10⁻²λ), in the UDM without thermal fluctuations. (c) Corresponding concentration profiles of the RDA vesicle from (b), where the pressure difference Δp can be read off. (d) Membrane tension calculated from the pressure difference, as well as the direct measurement by equilibration of the bilayer profile, with implicit and explicit fuel.

The corresponding radial profiles of the RDA vesicle are presented in Fig. 6(c). The reactant (precursor) profile, ϕ_R(r), clearly demonstrates the precursor enrichment inside the spherical vesicle. The pressure (or negative grandcanonical free-energy density, g) can be obtained from the Helmholtz free-energy density f of the UDM via Legendre transformation

Note that a more accurate pressure profile involves a proper spatial assignment and the distinction of the pressure in the directions perpendicular and parallel to the membrane.⁷⁹ In the following, however, we are only interested in (i) the difference, Δp, of the isotropic and spatially homogeneous pressure at the center of the vesicle and in the outside and (ii) the excess grandcanonical free energy, ΔG, with respect to the solution.

From the pressure difference and the Young-Laplace eqn (18) we can estimate the membrane tension, σ. Alternatively, we can approximate the membrane tension by that of a planar membrane of the same thickness, equilibrated without reactions. Noting that the membrane tension is the excess grandcanonical free energy, ΔG, per unit area, we estimate σ from 1D pressure profiles

where L_z and A denote the system extent perpendicular to the membrane and the membrane area, respectively.

Fig. 6(d) demonstrates the consistency of the two tension estimates according to eqn (18) and (20) for various reaction rates, both with implicit and explicit fuel. This validates that the precursor enrichment inside the vesicle results in the inflation of the vesicle. In the case of explicit fuel, the fuel concentration is inhomogeneous, decreased inside the vesicle, similar to the solvent concentration in Fig. 6(b). We observe that (i) the inflation increases with increasing reaction rates and (ii) for explicit fuel the tension is higher. (i) is explained by the fact that pressure can be released by diffusion of precursor through the membrane and more of it escapes within its lifetime for smaller reaction rates. (ii) Since less fuel is present inside of the RDA vesicle, less precursor reacts on the inside, raising its concentration and thus the membrane tension.

3.5 Metastable pore anomaly

When exploring RDA of vesicles from homogeneous solution with explicit fuel i.e., we allow for inhomogeneous fuel concentration, we observe an interesting anomaly where vesicles that emerge via pathway I, retain a metastable pore, which causes accelerated coarsening, see Fig. 7. This is related to the fact that inflation and hence increase in membrane tension is accelerated compared to homogeneously treated fuel. A pore arising during the vesicle formation process results from an interplay between two timescales: the reaction-dependent inflation time and the architecture-dictated time it takes for a planar bilayer to bend and close. This will be demonstrated in the following.

Fig. 7 (a) Dependence of the stationary micelle and vesicle radii on the reaction rate for the UDM after time tλ > 10³ and power law fit for data points where meta-stable pores are not observed. (b) Time evolution of vesicle morphology at r_b = 1.68 × 10⁻²λ. The moving vesicle with a metastable pore is highlighted by the yellow contour. At tλ ≈ 4.7 × 10² the vesicle collides with another vesicle and both fuse.

The dependence of the characteristic size on the reaction rates with explicit fuel is depicted in Fig. 7(a). Similar to the previous, implicit-fuel case, the scaling roughly agrees with eqn (17), R ∼ r_b^−1/2 with a measured exponent ν = −0.57 ± 0.06.

The vesicles originate exclusively via pathway II for large forward reaction rates, r_b ≥4 × 10⁻²λ. For smaller rates, when vesicles emerge via pathway I, however, the vesicle radii deviate from the expectation, eqn (17). In this case, disk-like micelles bend and form a vesicle with a pore, see Fig. 7. Remarkably, panel b illustrates that such a pore persists until the vesicle merges with a neighbor.‡

Vesicles formed by RDA are inflated by precursor, and the membrane is under tension. With a pore present, the pressure difference, Δp, can be reduced by exchanging precursor between the vesicle's inside and outside. If the pore is large, this exchange is effective. Thus, the pressure difference, Δp, decreases. This, in turn, reduces the membrane tension, σ, and the balance between the pore's line tension, Λ, and the membrane tension shifts towards shrinking the pore diameter, D_pore = 2σ/Λ.⁸⁰ If, on the contrary, the pore is small, the precursor exchange is limited, the pressure difference is increased and so is the membrane tension. This, in turn, tends to expand the pore size. As small pores expand and large pores shrink, a finite pore size is stabilized.

Additionally, small RDA vesicles with an open pore exhibit a directed movement opposite to the pore. This effect is depicted in Fig. 7(b) and can also be appreciated in the Movie udm-expl-fuel-tumbling-vesicle.mp4 of the ESI.† In the dense fluid of vesicles, the highlighted vesicle initially moves in the direction opposite to its pore. When the vesicle-free space (void) behind it becomes large, the movement slows down and stalls. Since RDA vesicles ideally pack on a dense, regular lattice, vesicle–vesicle interactions give rise to a force that pulls the vesicle back towards the center of the void. The vesicle tumbles and the direction of motion spontaneously changes. At tλ ≈ 4.7 × 10², the highlighted vesicle collides with a neighboring vesicle and both fuse. Such collision-fusion processes greatly enhance the coarsening dynamics and explain the observed increase in vesicle size as a function of reaction rate in Fig. 7(a) at a finite time.

In Fig. 8 we study the properties of a vesicle with a metastable pore, using the UDM and particle-based simulations. The simulation cell contains a single vesicle and we study different reaction rates. r_b. We monitor their average pore diameter and velocity, when the pore is present. Pores remain stable for high reaction rates because of the larger pressure difference and membrane tension. Pores shrink and become unstable as we decrease r_b. For r_b = 10⁻²λ the realization is shown in the Movie udm-metastable-pore-analysis.mp4 of the ESI.†

Fig. 8 Dependence of the pore diameter, D_pore (blue), and the average velocity of a vesicle with a pore (red) on the reaction rate in the UDM. The inset of the main figure shows the same quantities for the particle-based simulations. For three, selected parameter values, identified by (a–c) in the main figure, we depict the concentration fields of the tails, B, and the precursor, R, in the lower panels (a–c).

Additionally, Fig. 8 presents the velocity of RDA vesicles with an open pore, as observed in the UDM and the particle-based simulations. The velocity is opposite to the direction from the center of the vesicle to the pore. In the absence of a pore, i.e., at small r_f, the vesicles diffuse and the velocity is vanishingly small.

The presence of a pore in the vesicle allows for an osmotic-pressure-driven flux of the precursor along the direction from the vesicle center to the pore, reducing the precursor enrichment inside the vesicle. Incompressibility, in return, enforces that this outward flux is compensated by an inward flux of solvent, fuel, and amphiphile in the opposite direction. The latter slightly deforms the vesicle around the pore and induces the movement of the vesicle in the opposite direction of the pore.

Note, that vesicle inflation and pore stabilization are dynamic effects. In the above analysis, vesicles were initialized in the inflated state, allowing for a precise measurement of pore stability. In the dynamic system, a pore only becomes stable if sufficient inflation occurs on the time scales, on which the bilayer bends, forming a vesicle with a pore. This is why the effect is not visible in the simulations with homogeneously distributed fuel, with slower inflation, as well as in the dynamic, particle-based simulations. Alternatively, a pore could potentially form spontaneously in a closed vesicle under tension by thermal fluctuations.

4 Conclusion

Investigating the behavior of compartmentalizing amphiphiles in aqueous solution driven out of equilibrium by a coupling to a reaction cycle, we observe a variety of reaction-driven assembly (RDA) phenomena. The reaction cycle converts a hydrophilic precursor molecule (reactant, R) into a amphiphile (product, P) in the presence of a fuel (F). The amphiphilic product, P, assembles into micelles or vesicle but spontaneously converts back to the hydrophilic reactant, R. We focus on the regime, where a dense arrangement of micelles or vesicles is formed.

Using analytical considerations and simulations of the Uneyama–Doi model (UDM) and a particle-based model, we show that the structure initially forms by microphase separation between the hydrophilic components (solvent S and precursor R) and the amphiphilic product, P. This behavior can be described by an effective free-energy functional, analog to the Ohta–Kawaski model that describes microphase separation in the absence reactions.

Vesicles form in an initial, homogeneous solution by converting spherical micelles to disk-like micelles that bend and close to a vesicle (pathway I) at small reaction rates or via the direct flip-flopping of amphiphiles to the center of the micelle (pathway II) at high reaction rates. Along the former pathway I, we observed the formation of (meta)stable pores in the bilayer that result in the directed motion of vesicles due to the exchange of precursor between the vesicle's inside and its surrounding through the pore.

There are multiple characteristics that distinguish the RDA of vesicles from the formation of vesicles by amphiphilic molecules in equilibrium:

• The size distribution of the vesicles is narrow, and the average vesicle radius scales like with the reaction rate, i.e., it is responsive to thermodynamic parameters that influence the reaction rate.

• There is a fast, efficient transfer of building blocks between micelles that is based on the diffusion of molecules in its hydrophilic precursor state and allows for a rapid adjustment of vesicle sizes. This is in marked contrast to the protracted exchange of amphiphiles between aggregates (micelles or vesicles) in equilibrium self-assembly.

• In contrast to self-assembled vesicles in equilibrium, RDA results in inflated vesicles, characterized by a finite membrane tension. This membrane tension depends on the reaction rate, as well as on the density of aggregates. Vesicle tension, in turn, facilitates topological changes such as pore formation or the fusion of vesicles with membranes that are implicated in various transport processes in cells or subcellular compartments or could be exploited for targeted, tunable release of drugs.

These characteristics allow to control the vesicle assembly by external stimuli without changing the chemistry of the constituents.

Taking a reported diffusion coefficient of ⁸¹ for an amphiphilic diblock copolymer with R_e∼10² nm yields a common self-diffusion time of τ₀ ∼ 10⁻² s and hence the presented simulations cover times of t ∼ 30 s. The corresponding deactivation times in our case would be on the order of for the slowest reaction rates. Commonly reported values of ⁸² are higher, implying the lattice spacing would increase by a factor of . Notice, that lower molecular weights will require higher reaction rates.

Further topics of interest may include the study of fuel that is consumed by the reaction (instead of being a catalyst) and thus needs to be re-supplied. This may give rise to additional spatial inhomogeneities. Additionally, detailed reaction schemes that do not rapidly switch an entire block but allow for segment-wise reactions could be considered. We expect, however, that the qualitative differences between RDA of vesicles will remain unaltered.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was conducted within the Max Planck School Matter to Life supported by the German Federal Ministry of Education and Research (BMBF) in collaboration with the Max Planck Society. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (http://www.gauss-centre.eu) for providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC), as well as the North-German Supercomputing Alliance (HLRN).

References

1. E. Rideau, R. Dimova, P. Schwille, F. R. Wurm and K. Landfester, Chem. Soc. Rev., 2018, 47, 8572–8610.
2. R. Dimova, P. Stano, C. M. Marques and P. Walde, in The Giant Vesicle Book, ed. R. Dimova and C. M. Marques, CRC Press, Taylor & Francis Group, 2020.
3. B. M. Discher, Y.-Y. Won, D. S. Ege, J. C.-M. Lee, F. S. Bates, D. E. Discher and D. A. Hammer, Science, 1999, 284, 1143–1146.
4. K. Kita-Tokarczyk, J. Grumelard, T. Haefele and W. Meier, Polymer, 2005, 46, 3540–3563.
5. D. E. Discher, V. Ortiz, G. Srinivas, M. L. Klein, Y. Kim, D. Christian, S. Cai, P. Photos and F. Ahmed, Prog. Polym. Sci., 2007, 32, 838–857.
6. G. Srinivas, J. C. Shelley, S. O. Nielsen, D. E. Discher and M. L. Klein, J. Phys. Chem. B, 2004, 108, 8153–8160.
7. A. J. van Hell, C. I. C. A. Costa, F. M. Flesch, M. Sutter, W. Jiskoot, D. J. A. Crommelin, W. E. Hennink and E. Mastrobattista, Biomacromolecules, 2007, 8, 2753–2761.
8. V. Haridas, A. Rao Sapala and J. P. Jasinski, Chem. Commun., 2015, 51, 6905–6908.
9. F. Meng, Z. Zhong and J. Feijen, Biomacromolecules, 2009, 10, 197–209.
10. O. Onaca, R. Enea, D. W. Hughes and W. Meier, Macromol. Biosci., 2009, 9, 129–139.
11. F. Meng and Z. Zhong, J. Phys. Chem. Lett., 2011, 2, 1533–1539.
12. K. Göpfrich, I. Platzman and J. P. Spatz, Trends Biotechnol., 2018, 36, 938–951.
13. W. K. Spoelstra, S. Deshpande and C. Dekker, Curr. Opin. Biotechnol, 2018, 51, 47–56.
14. R. Bleul, R. Thiermann and M. Maskos, Macromolecules, 2015, 48, 7396–7409.
15. A. F. Mason and P. Thordarson, J. Polym. Sci., Part A: Polym. Chem., 2017, 55, 3817–3825.
16. P. Schwille, J. Spatz, K. Landfester, E. Bodenschatz, S. Herminghaus, V. Sourjik, T. J. Erb, P. Bastiaens, R. Lipowsky, A. Hyman, P. Dabrock, J.-C. Baret, T. Vidakovic-Koch, P. Bieling, R. Dimova, H. Mutschler, T. Robinson, T.-Y. D. Tang, S. Wegner and K. Sundmacher, Angew. Chem., Int. Ed., 2018, 57, 13382–13392.
17. T. Robinson, Adv. Biosyst., 2019, 3, 1800318.
18. M. Weiss, J. P. Frohnmayer, L. T. Benk, B. Haller, J.-W. Janiesch, T. Heitkamp, M. Börsch, R. B. Lira, R. Dimova, R. Lipowsky, E. Bodenschatz, J.-C. Baret, T. Vidakovic-Koch, K. Sundmacher, I. Platzman and J. P. Spatz, Nat. Mater., 2018, 17, 89–96.
19. A. Walther, Adv. Mater., 2020, 32, 1905111.
20. B. Rieß, R. K. Grötsch and J. Boekhoven, Chem, 2020, 6, 552–578.
21. S. Amano, S. Borsley, D. A. Leigh and Z. Sun, Nat. Nanotechnol., 2021, 16, 1057–1067.
22. D. Zwicker, A. A. Hyman and F. Jülicher, Phys. Rev. E, 2015, 92, 012317.
23. K. A. Rosowski, E. Vidal-Henriquez, D. Zwicker, R. W. Style and E. R. Dufresne, Soft Matter, 2020, 16, 5892–5897.
24. D. Zwicker, Curr. Opin. Colloid Interface Sci., 2022, 61, 101606.
25. D. Zwicker, R. Seyboldt, C. A. Weber, A. A. Hyman and F. Jülicher, Nat. Phys., 2017, 13, 408–413.
26. C. A. Weber, D. Zwicker, F. Jülicher and C. F. Lee, Rep. Prog. Phys., 2019, 82, 064601.
27. J. Bauermann, C. A. Weber and F. Jülicher, Ann. Phys., 2022, 534, 2200132.
28. J. Heckel, F. Batti, R. T. Mathers and A. Walther, Soft Matter, 2021, 17, 5401–5409.
29. M. Tena-Solsona, J. Janssen, C. Wanzke, F. Schnitter, H. Park, B. Rieß, J. M. Gibbs, C. A. Weber and J. Boekhoven, ChemSystemsChem, 2021, 3, e2000034.
30. J. Boekhoven, A. M. Brizard, K. N. K. Kowlgi, G. J. M. Koper, R. Eelkema and J. H. van Esch, Angew. Chem., Int. Ed., 2010, 49, 4825–4828.
31. J. Boekhoven, J. M. Poolman, C. Maity, F. Li, L. van der Mee, C. B. Minkenberg, E. Mendes, J. H. van Esch and R. Eelkema, Nat. Chem., 2013, 5, 433–437.
32. J. Heckel, S. Loescher, R. T. Mathers and A. Walther, Angew. Chem., Int. Ed., 2021, 60, 7117–7125.
33. B. Klemm, R. W. Lewis, I. Piergentili and R. Eelkema, Nat. Commun., 2022, 13, 6242.
34. C. Donau, F. Späth, M. Sosson, B. A. K. Kriebisch, F. Schnitter, M. Tena-Solsona, H.-S. Kang, E. Salibi, M. Sattler, H. Mutschler and J. Boekhoven, Nat. Commun., 2020, 11, 5167.
35. B. G. P. van Ravensteijn, I. K. Voets, W. K. Kegel and R. Eelkema, Langmuir, 2020, 36, 10639–10656.
36. R. W. Lewis, B. Klemm, M. Macchione and R. Eelkema, Chem. Sci., 2022, 13, 4533–4544.
37. K. Dai, M. Tena-Solsona, J. R. Fores, A. M. Bergmann and J. Boekhoven, Nanoscale, 2021, 13, 19864–19869.
38. C. Sharma and A. Walther, Angew. Chem., Int. Ed., 2022, 61, e202201573.
39. X. Lang, U. Thumu, L. Yuan, C. Zheng, H. Zhang, L. He, H. Zhao and C. Zhao, Chem. Commun., 2021, 57, 5786–5789.
40. M. A. Würbser, P. S. Schwarz, J. Heckel, A. M. Bergmann, A. Walther and J. Boekhoven, ChemSystemsChem, 2021, 3, e2100015.
41. S. Maiti, I. Fortunati, C. Ferrante, P. Scrimin and L. J. Prins, Nat. Chem., 2016, 8, 725–731.
42. C. Wanzke, A. Jussupow, F. Kohler, H. Dietz, V. R. I. Kaila and J. Boekhoven, ChemSystemsChem, 2020, 2, e1900044.
43. E. Poros-Tarcali and J. Perez-Mercader, Soft Matter, 2021, 17, 4011–4018.
44. G. Cheng and J. Perez-Mercader, Chem, 2020, 6, 1160–1171.
45. T. Uneyama and M. Doi, Macromolecules, 2005, 38, 196–205.
46. K. C. Daoulas and M. Müller, J. Chem. Phys., 2006, 125, 184904.
47. L. Schneider and M. Müller, Comput. Phys. Commun., 2019, 235, 463–476.
48. T. Uneyama, J. Chem. Phys., 2007, 126, 114902.
49. T. Ohta and K. Kawasaki, Macromolecules, 1986, 19, 2621–2632.
50. O. Dreyer, G. Ibbeken, L. Schneider, N. Blagojevic, M. Radjabian, V. Abetz and M. Müller, Macromolecules, 2022, 55, 7564–7582.
51. L. Leibler, Macromolecules, 1980, 13, 1602–1617.
52. M. L. Huggins, J. Chem. Phys., 1941, 9, 440.
53. P. J. Flory, J. Chem. Phys., 1942, 10, 51–61.
54. P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys., 1977, 49, 435–479.
55. G. Wang, Y. Ren and M. Müller, Macromolecules, 2019, 52, 7704–7720.
56. J. Rottler and M. Müller, ACS Nano, 2020, 14, 13986–13994.
57. A. Ianiro, H. Wu, M. M. J. van Rijt, M. P. Vena, A. D. A. Keizer, A. C. C. Esteves, R. Tuinier, H. Friedrich, N. A. J. M. Sommerdijk and J. P. Patterson, Nat. Chem., 2019, 11, 320–328.
58. J. J. Christensen, K. Elder and H. C. Fogedby, Phys. Rev. E, 1996, 54, R2212–R2215.
59. M. Müller and J. C. O. Rey, Mol. Syst. Des. Eng., 2018, 3, 295–313.
60. S. C. Glotzer, E. A. Di Marzio and M. Muthukumar, Phys. Rev. Lett., 1995, 74, 2034–2037.
61. X. He and F. Schmid, Phys. Rev. Lett., 2008, 100, 137802.
62. J. C. Shillcock, Langmuir, 2012, 28, 541–547.
63. J. Li, H. Zhang, F. Qiu and A.-C. Shi, Phys. Rev. E, 2013, 88, 012719.
64. T. Uneyama and M. Doi, Macromolecules, 2005, 38, 5817–5825.
65. X. He and F. Schmid, Macromolecules, 2006, 39, 2654–2662.
66. J. Leng, S. U. Egelhaaf and M. E. Cates, Europhys. Lett., 2002, 59, 311.
67. J. Leng, S. U. Egelhaaf and M. E. Cates, Biophys. J., 2003, 85, 1624–1646.
68. T. M. Weiss, T. Narayanan, C. Wolf, M. Gradzielski, P. Panine, S. Finet and W. I. Helsby, Phys. Rev. Lett., 2005, 94, 038303.
69. R. Takahashi, T. Narayanan, S.-I. Yusa and T. Sato, Macromolecules, 2022, 55, 684–695.
70. J. Du and Y. Chen, Macromolecules, 2004, 37, 5710–5716.
71. E. V. Korchagina, X.-P. Qiu and F. M. Winnik, Macromolecules, 2013, 46, 2341–2351.
72. J. Gummel, M. Sztucki, T. Narayanan and M. Gradzielski, Soft Matter, 2011, 7, 5731–5738.
73. Z. Ye, Z. Wu and A. Jayaraman, JACS Au, 2021, 1, 1925–1936.
74. V. Chappa, Y. Smirnova, K. Komorowski, M. Müller and T. Salditt, J. Appl. Crystallogr., 2021, 54, 557–568.
75. J. Hoshen and R. Kopelman, Phys. Rev. B, 1976, 14, 3438–3445.
76. J. Hoshen, M. W. Berry and K. S. Minser, Phys. Rev. E, 1997, 56, 1455–1460.
77. C. Donau, F. Späth, M. Stasi, A. M. Bergmann and J. Boekhoven, Angew. Chem., Int. Ed., 2022, 61, e202211905.
78. R. Lipowsky, Adv. Biol., 2022, 6, 2101020.
79. C. L. Ting and M. Müller, J. Chem. Phys., 2017, 146, 104901.
80. T. V. Tolpekina, W. K. den Otter and W. J. Briels, J. Chem. Phys., 2004, 121, 8014–8020.
81. J. C.-M. Lee, M. Santore, F. S. Bates and D. E. Discher, Macromolecules, 2002, 35, 323–326.
82. C. Donau and J. Boekhoven, Trends Chem., 2022, 5, 45–60.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sm00876b

‡ In rare occasions the pore closes by itself. Furthermore, for reaction rates r_b = 2.0 × 10⁻²λ and r_b = 2.8 × 10⁻²λ a cylindrical micelle with cup-shaped ends persists for long times, up to tλ = 2.5 × 10³, where a pure cylindrical micelle remains. Videos of these kinetics are prepared for two exemplary reaction rates in the ESI,† udm-expl-fuel-rb1.67e-2.mp4 and udm-expl-fuel-rb2.00e-2.mp4.

---