## Abstract

By developing new computational models, we examine how enzymatic reactions on an underlying surface can be harnessed to direct the motion and organization of reagent-laden microcapsules in a fluid-filled microchannel. In the presence of appropriate reagents, surface-bound enzymes can act as pumps, which drive large-scale fluid flows. When the reagents diffuse through the capsules’ porous shells, they can react with enzymatic sites on the bottom surface. The ensuing reaction generates fluid density variations, which result in fluid flows. These flows carry the suspended microcapsules and drive them to aggregate into “colonies” on and near the enzyme-covered sites. This aggregation continues until the reagent has been depleted and the convection stops. We show that the shape of the assembled colonies can be tailored by patterning the distribution of enzymes on the surface. This fundamental physicochemical mechanism could have played a role in the self-organization of early biological cells (protocells) and can be used to regulate the autonomous motion and targeted delivery of microcarriers in microfluidic devices.

- Collective behavior
- microcapsules
- self-organization
- protocell aggregation
- colloidal transport
- enzymatic chemical micropumps
- computational modeling

## INTRODUCTION

Polymeric microcapsules encompass an outer shell and a hollow or fluid-filled core, which can hold a broad range of cargo (*1*–*5*). Given these defining characteristics, microcapsules can serve as models for the earliest biological cells, referred to as protocells, or as microcarriers that can deliver their cargo to specific, targeted sites (for example, in vivo or in microfluidic devices). In both contexts, fundamental studies on the behavior of microcapsules in solution can help address important, unresolved issues. With respect to protocells, it is unclear how the earliest cells, which lacked complex biochemical machinery, could communicate and, thereby, self-organize into larger colonies. This level of self-organization was vital for the cells’ chances of survival and, ultimately, the formation of differentiated structures, such as organelles. Although the diffusion of molecules through the protocells’ porous membranes would have enabled communication among these cells (*6*), the mechanisms that contributed to the collective motion and assembly of the protocells into larger aggregates remain unknown.

In the context of harnessing the microcapsules for targeted delivery, there have been significant advances in designing self-propelled microscopic particles (*7*–*12*). However, there remains a need for effective means to both “fuel” the microcapsules’ self-sustained motion and direct the microcarriers to specific sites and, thus, allow the encapsulated cargo to be delivered in a “programmable” way. For microfluidic applications, it is highly beneficial for the capsules to undergo controllable motion and delivery so that devices could operate in a relatively autonomous manner and, thus, be significantly more portable.

Here, we develop a theoretical model and computer simulations for chemical-releasing microcapsules that undergo self-propelled motion on a reactive substrate to shed light on possible mechanisms that enabled protocells to aggregate into colonies, and provide guidelines for driving capsules to autonomously localize at specific regions in a microchannel. The studies are inspired by recent experimental findings that show that, in the presence of their specific substrates or promoters, surface-bound enzymes act as pumps that drive large-scale fluid flows (*13*). Moreover, essentially all surface-anchored enzymes act as pumps when they turn over their substrates (*13*). Building on these observations, we model microcapsules in solution that gradually release a chemical reagent to a surface that contains a patch of enzymes. The released reagents activate the enzymatic pumps, which generate flows that transport the microcapsules, thereby establishing a useful feedback loop between the pumps and microcapsules.

In the systems considered here, the interactions between the reagents and catalytic patch lead to an exothermic reaction, which results in the decomposition of reagents and the formation of a product. The reaction thus gives rise to spatial variations in the fluid density that drives the convective flow. We focus on the case where the reaction products are less dense than the reagents; in this scenario, an upward flowing fluid forms above the enzyme-coated patch and drags the capsules toward the center of the patch, driving their self-organization into a large cluster. Hence, the model developed here captures the following cascade of interconnected dynamic events: (i) the release of reagents from the mobile capsules, (ii) the decomposition of the reagents through an exothermic reaction at the surface, (iii) the consequent changes in the reagent concentration and temperature within the fluid, and (iv) the resulting fluid-induced movement of the capsules in the solution.

In the prebiotic environment, minerals and inorganic surfaces provided viable catalysts for chemical reactions with protocells (*14*). Given the fundamental, physicochemical nature of the events described here, the resulting reactions could have driven the motion and self-assembly of the protocells. Notably, the time required for diffusive intercellular communications (scaling with the square of the characteristic distance) is greatly reduced if the colony is densely packed. The aggregation of the reagent-laden capsules described below yields a 10-fold increase in the areal concentration of capsules and a corresponding 10-fold reduction in time required for diffusive exchange of signals between neighboring capsules. This reduction in the communication time might have provided protocells with a greater chance to develop collective behavior (*15*–*17*).

With respect to microfluidic applications, the enzymatic pumps can be harnessed to both fuel and direct the motion of the microcarriers and thus enable the efficient operation of fluidic devices. Such pumps become active only when specific analytes are present in the solution and, hence, allow the system to also act as a sensor. In effect, with the aid of these catalytic sites, sensing and controlled delivery can be achieved in one device.

We begin by formulating a model for the capsule’s shell, the equations governing the fluid dynamics, and the motion of the capsules within the fluid. Using this model, we then examine the effects of varying the capsule size, shell permeability, and size of the enzyme-coated patch on convective aggregation. Finally, we demonstrate that the capsules can be assembled into colonies of different shapes by tailoring the design of the enzyme-covered patches. The difficulty of reproducing, or even fully knowing, the prebiotic conditions is worth noting; hence, computational modeling can provide a valuable means of probing how variations in the local environment can promote the self-organization of simple cellular entities.

## RESULTS

### Theoretical modeling

The spherical microcapsules are considered to comprise a multilayered polymeric shell that encloses a core, which contains the dissolved chemical reagent (see Fig. 1). Such microcapsules are commonly synthesized using the layer-by-layer (LbL) approach (*18*–*21*). The permeability of the LbL capsules’ shell can be tailored (*2*) to tune the rate of release of the reagent into the surrounding solution. The diameter of the capsules is taken to be tens of microns, and thus, they are comparable in size to biological cells. We assume that the spherical shell is nondeformable and ignore perturbations of the flow due to the capsules.

To describe the release of the chemical reagent out of the capsule (*22*), we model the capsule as a shell of thickness *d* surrounding a spherical volume of radius *R*. We assume that the reagent’s diffusivity *D*_{S} through the shell is much smaller than its diffusivity *D* in the surrounding solution and that the reagent concentration *C* outside the capsule is negligible relative to the internal concentration *C*^{in}. These assumptions allow us to decouple the mathematical descriptions for the reagent diffusion through the capsule shell and its diffusion in the rest of the fluid domain. Combining the laws for conservation of mass and Fickian diffusion, the rate of change of the amount *M* (in moles) of the encapsulated reagent is modeled as(1)where *P* = *D*_{S}/*d* is the shell permeability and *S* = 4π *R*^{2} is the surface area of the capsule. Writing *M* = *C*^{in}4π *R*^{3}/3, the above equation reduces to , which has the solution *M* = *M*_{0} exp(−3 *P t*/*R* ). The rate of release of reagent into the fluid is, therefore, .

For simplicity, we assume that capsules do not perturb the fluid flow. The motion of each capsule is the superposition of the advection with the local fluid velocity **u** and a contribution determined by local forces acting on capsules. These forces include gravity and excluded volume interactions. The gravity, described by the vector **g**, results in the sedimentation velocity . Here, ρ_{C} is the average mass density of the capsules and *η* is the dynamic viscosity of the host fluid. The excluded volume associated with each capsule is modeled through repulsive interactions between neighbors and between capsules and the walls, as discussed below. In our treatment of the motion of capsules with radius *R* = 10 to 50 μm over distances ~10^{−3} m, we ignore Brownian motion as this would contribute root mean square displacements of the order of ~10^{−6} m over the simulated time scale of a few hours.

Our computational domain is a three-dimensional box {(*x*, *y*, *z*) : 0 ≤ *x*, *z* ≤ *L*, 0 ≤ *y* ≤ *H*} bounded by top and bottom walls at *y* = 0 and *y* = *H*, respectively (see Fig. 1); we impose periodic boundary conditions in the lateral (*x* and *z*) directions. The fluid within this domain has a kinematic viscosity *ν*. The density of this fluid is altered when the reagent diffuses out of the capsules. The density of this mixture can be approximated as ρ = ρ_{0}(1 + β_{C}*C* + β_{T}(*T* − *T*_{0})), where ρ_{0} is the density of the pure solvent at some reference temperature *T*_{0} and variations due to changes in reagent concentration *C* and temperature *T* are characterized by the solutal β_{C} and thermal β_{T} expansion coefficients, respectively. The system is further characterized by the following variables: the fluid velocity **u** = (*u*_{x}, *u*_{y}, *u*_{z}) and pressure *p*, the local reagent concentration *C*, the local temperature *T*, and the position of the *N* capsules, which is specified by the vectors **r**_{i}(*t*) = (*x*_{i}(*t*) , *y*_{i}(*t*), *z*_{i}(*t*) ).

The dynamics of this multicomponent system are governed by the following: the continuity equation, Eq. 2; the Navier-Stokes equation [in the Boussinesq approximation (*23*)], Eq. 3; the expression for the diffusion of reagent, Eq. 4; the heat conduction equation, Eq. 5; and the equation for capsule motion, Eq. 6. These respective equations are written as(2)(3)(4)(5)(6)Here, ∇ is the spatial gradient operator, *g* is gravity, **e** = (0, 1, 0) is a vector specifying the direction of the gravity force along the *y* coordinate, and *V* = |**V**| is the speed of capsule sedimentation. Additionally, μ = (6 π η *R*)^{− 1} is the capsule mobility, specifies the vector from the capsule to the closest point on the wall, and and are forces arising from the steric repulsion between capsules and between the capsules and walls, respectively. In the latter expressions for the force, we use the following Morse potentials(7)(8)The parameters ε and ω characterize the respective potential strength and width.

To maximize the effect of convection, we assume that the enzymes cover certain areas on the bottom wall; the reagent-laden capsules sediment to this bottom wall. The rate of catalytic chemical reaction per unit area is described using the Michaelis-Menten relation(9)where the maximal reaction rate *r*_{max} = *k*_{cat}[*E*] incorporates the maximum reaction rate per enzyme molecule, *k*_{cat}, and the enzyme concentration, [*E*], in moles per unit area. *K*_{M} is the substrate concentration at which the reaction rate is half of *r*_{max}. Using the rate in Eq. 9 and the diffusion flux *D*∇*C*, the reagent decomposition at an enzyme-covered site on the bottom surface *y* = 0 can be described by the boundary condition .

Similarly, the heat flux κ∇*T* was used to obtain the heat generated during the reaction: . Here, we use the thermal conductivity of water κ = 0.58 (W/Km) and reaction enthalpy ΔH = 100 (kJ/mol) for hydrogen peroxide. Elsewhere on the bottom surface, and on the top surface, we assume zero-flux conditions for the chemical concentration and no-slip conditions for the fluid flow. The temperature on the top surface is held constant at the reference temperature *T*(*y* = *H*) = *T*_{0}, which is also the initial temperature throughout the domain. On the bottom surface, we impose the zero-flux condition for temperature. The full set of boundary conditions is given by(10)(11)(12)(13)Here, *A* denotes the enzyme-coated region of the bottom surface at *y* = 0 and *B* denotes the uncoated portion of the surface. In Materials and Methods, we define relevant dimensionless parameters and describe the numerical procedures used to solve the above governing equations. The initial reagent concentration in the fluid was C(*t* = 0) = 0, and the capsules were initially randomly distributed with uniform density throughout the domain. In our approach, we ignored friction and potential adhesion between capsules and the container walls. In practice, such effects may slow down the capsules and hinder aggregation dynamics.

### Choice of chemicals and system geometry

It is beyond the scope of this study to systematically vary each dimensionless parameter. Instead, we choose a particular realizable system to demonstrate the physical mechanisms underlying the aggregation phenomenon. In what follows, we specify a particular chemical reaction and explore the influence of the capsule size, shell permeability, and geometry of the enzyme-coated patch.

As a representative example of an enzymatic chemical reaction, we choose the decomposition of hydrogen peroxide (in an aqueous solution) by catalase to form oxygen and water; the reaction is well studied, and thus, we can obtain experimental values (*13*) for the simulation parameters. In particular, the maximum reaction rate per catalase molecule is *k*_{cat} = 8.48 × 10^{5} s^{−1} (taking into account four active sites per molecule) and the substrate concentration at half-maximal rate is *K*_{M} = 0.093 M. Taking the radius of a catalase molecule to be *R*_{cat} = 5.1 × 10^{−9} m, we estimate the areal concentration of a monolayer of enzyme to be on the order of [*E*] = 2 × 10^{−8} mol/m^{2}. The reaction has an enthalpy (*13*) of ΔH = 100 kJ/mol. Hydrogen peroxide diffusivity is taken to be *D* = 10^{−9} m^{2}/s.

Water has a density of ρ_{0} = 10^{3} kg/m^{3}, a kinematic viscosity of *ν* = 10^{−6} m^{2}/s, and a thermal conductivity of κ = 0.58 W/Km. For a mixture of *N*^{s} solutes at low concentrations *C*^{j}, *j* = 1, 2, …, *N*^{s}, the fluid density can be approximated as , where β_{T} ≈ 2 × 10^{−4} 1/K is the thermal expansion of water and ρ_{0} is the density of pure water at temperature *T*_{0}. For the catalase reaction system, the relevant solutes are hydrogen peroxide and oxygen, which have expansion coefficients and (*24*), respectively. Because the density contribution from oxygen is much less than that from hydrogen peroxide, we neglect oxygen in our model and consider only one solute, hydrogen peroxide, with concentration denoted by *C*.

Choosing the capsule size to be on the order of tens of micrometers, *R* = 10 to 50 μm, which is comparable to the typical size of bacteria and many experimentally fabricated capsules (*1*, *2*), we assume that each capsule initially contains a 30% hydrogen peroxide solution (*C*_{0} = 14 M). The permeability (*P*) of capsules varies from 10^{−5} to 10^{−10} m/s, consistent with experimentally realizable values (*2*). For simplicity, we consider a constant capsule sedimentation speed based on the interior fluid density when half of the reagent has been released, that is, ρ_{C} = ρ_{0}(1 + β_{C}*C*_{0}/2) ≈ 1.1ρ_{0}. This choice does not significantly affect the simulation outcomes because, in most cases, the capsules rapidly sediment to the bottom surface and remain in contact with the surface thereafter. The parameters describing the hard-core repulsion between the capsules and walls were set to ε^{C−C} = ε^{C−W} = 1.2 × 10^{−14} J and ω^{C−C} = ω^{C−W} = 2 × 10^{4} m^{−1}. With these parameters, the center-to-center distance between capsules was always at least 98% of a capsule diameter and the distance between the capsule center and the walls was always at least 98% of the capsule radius. (Note that for numerical stability, the repulsive potentials should not be made too stiff.)

For the above parameters, we find that a few hundred capsules contain enough reagent to induce aggregation over a length scale of a few millimeters. To demonstrate this process, we perform simulations in a periodic rectangular domain with height *H* = 1 mm and equal horizontal sizes *L* = 4 mm, as shown in Fig. 1. A square region on the bottom, indicated by the light blue color, represents a catalase-coated patch with area *A* = 1 × 1 mm^{2}.

## DISCUSSION

The efficiency of capsule aggregation by convective transport depends on the magnitude of the characteristic velocity of the fluid flow, which is generated by density variations in the solution. In addition to the geometry of the microchannel, there are a few other factors that control the flow velocity. First, the magnitude of the density variation is determined by the respective solutal and thermal expansion coefficients β_{C} and β_{T} through the associated Grashof numbers and . Second, decreases in the density of the fluid that arise from the decomposition of the reagent depend on the rate of decomposition (described by the parameter ) at the enzyme-coated surface. Enzymes that exhibit higher rates of catalysis provide flows with higher velocities. Third, the amount *N M*_{0} of fuel released into the domain is set by the number of capsules *N* and the amount of reagent *M*_{0} initially contained in each capsule. Finally, the capsule size and shell permeability control the rate of release of the reagent into the external fluid and should, therefore, also affect the aggregation dynamics. Before discussing the effects of varying the capsule size, permeability of the capsule’s shell, and size of the catalytic patch, we first describe the fundamental physical mechanisms that give rise to the capsule motion and aggregation in this system.

### Physical mechanism of convective aggregation

The main stages of the aggregation process are demonstrated in Fig. 2, where the left and right panels show the vertical (side view) and horizontal (top view) cross sections of the computational domain, respectively. Initially, the domain is assumed to be filled with water, and *N* = 500 capsules with radius *R* = 20 μm are homogeneously distributed throughout the domain, as shown in Fig. 2A. In the absence of convective flow, the capsules sediment with a velocity controlled by the dimensionless number .

As the permeable capsules undergo sedimentation, they release the encapsulated reagent, which increases the solution density. When the released reagent reaches the enzyme-coated region on the bottom wall, it is decomposed by the enzyme. The exothermic decomposition reaction decreases the fluid density by both lowering the local solute concentration and raising the temperature; both these effects result in an upward convective flow. The flow field is indicated by arrows in the left panel of Fig. 2B; the maximal velocity is achieved above the center of the enzyme-covered patch where the buoyant fluid is seen to rise upward.

In the case of many small capsules (which have a low sedimentation speed), the flow may carry the capsules around in loops by the convective vortexes that are formed as a result of the exothermic decomposition reaction. As the reagent is consumed, however, the flow velocity decreases and all capsules eventually settle to the bottom surface. The absence of friction between the capsules and the bottom wall allows the flow to drag the capsules, which follow the streamlines toward the center of the patch.

After all the reagent is consumed, the fluid motion stops. At this stage, the capsules are found to be concentrated on and around the patch as shown in Fig. 2C. The characteristic distance between nearest neighbors can be conveniently defined in terms of the areal concentration of capsules. Counting the number of capsules *m* located on a given square patch of area *A*′ centered on the coordinates (*x*,*z*) on the bottom surface, we can characterize the efficiency of the aggregation process by comparing the areal concentration *n* = *m*/*A*′ relative to the initial uniform distribution value *n*_{0} = *N*/*L*^{2}. The relative capsule areal concentration(14)along the horizontal line passing through the center of the domain is shown with a red line above each image on the left in Fig. 2. The control area was chosen to be *A*′ = 0.4 × 0.4 mm^{2}, which was large enough to smooth out some of the fluctuations associated with counting small numbers of discrete particles. According to the definition in Eq. 14, the relative concentration *n*/*n*_{0} is bounded from above by the value *n*/*n*_{0} = *L*^{2}/*A*′, describing the situation when all the capsules (*m* = *N*) reside within a region of area *A*′. To describe the effect of convection on the spatial distribution of capsules, we track the evolution of the relative capsule concentration over time. As shown in Fig. 2A, a peak in capsule concentration forms at the center of the patch and grows taller with time, reflecting the evolution of the aggregation process. We expect that if friction between the capsule and bottom surface were to be included, then the aggregation dynamics would be slower and the final peak in the capsule concentration distribution would be lower and broader. Details of the aggregation process are demonstrated in movie S1.

### Effect of capsule size on aggregation dynamics

The convective aggregation of the capsules is accompanied by the depletion of the reagent and corresponding changes in fluid velocities and temperature. Because the amount of the reagent that can be encased within a capsule depends on its size, we anticipate that the behavior of the system will depend on *R*, the capsule radius. Hence, in Fig. 3, we examine the effect of varying *R* on the temporal evolution of the following salient features: the relative areal concentration of the capsules *n*/*n*_{0} (above the enzyme-coated patch *A*′ = 1 × 1 mm^{2}), the fluid velocity in the center of the domain, the maximal reagent concentration (within the domain), and the maximal temperature change (relative to the reference temperature).

The slopes of the *n*/*n*_{0} curves in Fig. 3A indicate that larger capsules aggregate faster than the smaller ones. This difference in behavior occurs for two reasons. First, because of the no-slip boundary condition on the bottom wall, the horizontal flow speed increases approximately linearly with height from this surface. The center of a larger capsule lies further from the wall than the smaller one and, thus, experiences this higher fluid velocity. Second, a larger capsule releases a reagent at a higher rate (from its larger surface area as indicated in Eq. 1) than a smaller one and, thus, generates a faster flow field (Fig. 3B). Notably, in all cases, a decrease in the reagent concentration due to the decomposition reaction (Fig. 3C) leads to a decrease in the fluid velocity (Fig. 3B) and a decrease in the maximal temperature change (Fig. 3D).

Aggregation dynamics, quantified by the increase of the capsule concentration *n*/*n*_{0}, can be limited by two mechanisms: insufficient encapsulated fuel, which results in weak fluid flow, and hard-core repulsion between neighboring capsules, which prevents an arbitrarily large number of capsules from occupying a given region. To estimate the upper bound for concentrations due to repulsion, we assume that the capsules on the bottom form a two-dimensional close-packed hexagonal structure. This allows about 91% of the surface to be covered by capsules and gives a limiting number density *n*^{lim} = 0.91/(π *R*^{2}). The aggregation factor *n*/*n*_{0}, therefore, has an upper limit that depends on the capsule radius, the total number of capsules, and the size of the domain.

The dynamics shown in Fig. 3A are for the aggregation of *N* = 500 capsules on the bottom of a periodic computation domain of size *L* = 4 mm in both lateral directions. As indicated in the legend to Fig. 3, we consider the following values for the capsule radius: *R* = 10, 15, 20, 25, 30, 40, and 50 μm. The respective limiting values of the relative areal concentrations *n*^{lim}/*n*_{0} are 92.69, 41.19, 23.17, 14.83, 10.29, 5.79, and 3.71. The final simulated values of *n*/*n*_{0} for the cases with *R* = 40 and 50 μm are close to the respective upper bounds, suggesting that the aggregate is nearly close-packed in these two cases. For the cases with smaller capsule radii, the maximum packing density is not reached, indicating that aggregation can still be enhanced, for example, by supplying more fuel (increasing the total number of capsules *N*). To observe this trend, we plot in Fig. 4 values of maximal relative concentrations *n*/*n*_{0} (with *A*′ = 1 × 1 mm^{2}), fluid velocity in the middle point of the domain, maximal reagent concentration, and temperature change during the reaction as functions of capsule size for four representative values of *N*.

Figure 4A demonstrates that for each fixed *N*, there is an optimal radius *R*^{opt}(*N*) that provides the maximal areal concentration *n*/*n*_{0} of capsules. This optimum is a result of competition between the aggregation strength increasing with capsule size (amount of fuel) and the hard-core repulsion between neighbors limiting the number of capsules per unit of area. When *N* decreases, the optimal radius corresponding to the maximal concentration *n*/*n*_{0} increases, suggesting that maximal values of areal concentration can be achieved by the aggregation of fewer capsules of larger diameters. Figure 4A also shows that differences in the behavior of areal concentration *n*/*n*_{0} divide capsules into two groups separated by *R*^{opt}(*N*). The “small” radius capsules, with *R* < *R*^{opt}(*N*), do not experience a limitation from the hard-core repulsion because *n*/*n*_{0} < *n*^{lim}/*n*_{0} and the concentration *n*/*n*_{0} increases for larger values of *N*. For these capsules, the values of *n*/*n*_{0} presented in the figure are limited by the total amount of the encapsulated reagent. The “large” capsules with *R* > *R*^{opt}(*N*) and concentration *n*/*n*_{0} ~ *n*^{lim}/*n*_{0} are limited by repulsion between the neighbors so that the values of relative concentration *n*/*n*_{0} decrease for larger numbers of capsules *N* (though absolute concentration *n* remains the same).

Figure 4 (B to D) shows that increasing the capsule radius and number of capsules monotonically increases the maximal flow speed, reagent concentration, and temperature, respectively. The typical values of the reagent concentration (Fig. 4C) and temperature (Fig. 4D) indicate that the density contribution from solutal effects β_{C}*C* ~ 10^{−3} (β_{C} ≈ 0.011 liter/mol) is a few orders of magnitude larger than the contribution from thermal effects β_{T}*T* ~ 10^{−6} to 10^{−5} (β_{T} ~ 10^{−4} K^{−1}). Therefore, in the case of hydrogen peroxide decomposition by catalase, convective flow is predominantly driven by solutal density variations.

### Effect of shell permeability on capsule aggregation

To examine the influence of the rate of reagent release from the capsules on convective aggregation, we perform simulations for a set of representative shell permeabilities *P* = 10^{−5}, 10^{−6}, 10^{−7}, 10^{−8}, 10^{−9}, and 10^{−10} m/s, which are within the range experimentally realized by Antipov and Sukhorukov (*2*). The ratio between permeability of the multilayered polymeric shell and the capsule radius determines the rate ∝ *P*/*R* of reagent release from the capsules. The corresponding characteristic release times are τ = *R*/(3*P*) = 0.003, 0.03, 0.3, 3, 30, and 300 min. The resulting capsule aggregation, characterized by the evolution of the areal capsule concentration *n*/*n*_{0} (with *A*′ = 1 × 1 mm^{2}) and maximal reagent concentration in the domain, is presented in Fig. 5 (A and B, respectively).

For low rates of reagent release, the generated convection is too weak to bring capsules together, as demonstrated by the slow increase of capsule concentration *n*/*n*_{0} for *P* = 10^{−10} m/s. Aggregation becomes progressively faster as the shell permeability is increased up to *P* = 10^{−8} m/s. For the highest rates of reagent release, when the shell permeability is *P* ≥ 10^{−8} m/s and release times are *τ* < 3 min, the capsules release most of their encapsulated fuel before they reach the bottom wall. This situation is effectively equivalent to starting from a uniform reagent concentration throughout the domain and is insensitive to the precise value of *P*. Hence, aggregation dynamics are similar in the cases where *P* = 10^{−5}, 10^{−6}, 10^{−7}, and 10^{−8} m/s.

Note that in the case when the capsules are initially dispersed over the top surface (as a result of a process that introduces the capsules into the system), encapsulation of the reagent can accelerate fuel delivery to the enzyme-coated regions on the bottom and, thus, expedite generation of convection. In particular, in the domain with height *H* =1 mm, most of the capsules that sediment (in the absence of convection) with velocity ~10^{−5} to 10^{−4} m/s should reach the bottom wall within a time ~10 to 100 s, whereas the diffusive mechanism, controlled by characteristic diffusion time *H*^{2}/*D* = 1000 s, delivers fuel to the enzyme much more slowly.

### Effect of the size of enzyme-coated patch on capsule aggregation

The enzyme-coated patch acts as an “engine,” where substrate is consumed and the largest fluid density variations are generated. Therefore, in addition to the type of enzyme, the area of the patch is expected to influence the generation of convective flow. Distributions of capsules obtained during simulations with different areas of the square patch are shown in the bottom panels in Fig. 6. The top panels show the corresponding relative capsule concentrations *n*/*n*_{0} (with *A*′ = 0.4 × 0.4 mm^{2} and along the horizontal line *z* = 2 mm) developed after the reagent was consumed and convection was stopped.

For a given set of parameters characterizing fluid domain, reagent, enzyme, and capsules, there is an optimal size of the enzyme-coated region that provides the fastest aggregation and the highest concentration of capsules. Small patches (the first case in Fig. 6 with ) do not provide enough surface area for reagent decomposition and, therefore, limit the intensity of convection. In the other extreme case, large patches (the last case in Fig. 6 with ) consume the reagent quickly, but the reaction is spread out over a larger area; hence, the resulting density gradients, which drive convection, are smaller. The patch size that yields the maximal aggregation rate along with maximal capsule concentrations should be comparable to the height of the microchannel, as the second () and third () cases show in Fig. 6.

### Capsule assembly on different patterns

The above finding suggests that the assembly of the capsules could be regulated by patterning the spatial distribution of the enzyme on the bottom wall, and in this way, the capsules could be delivered in a specific configuration and at specific sites in the domain. We chose three representative patterns (circle, hollow square, and crankshaft) to test the ability of capsules to self-assemble into structured clusters. In Fig. 7, we show the stages of assembly onto a hollow square enzyme pattern for the case of *N* = 500 capsules with radius *R* = 20 μm. The enzyme-coated region is indicated in the top view (right column) by small, blue squares.

Examination of the time-dependent convective flow shown in Fig. 7 reveals that the aggregation of capsules occurs in two stages. At the initial stage, when the reagent is abundant, the flow structure consists of four (unequal) counter-rotating vortexes, which are marked in the vertical cross section through the center of the pattern (Fig. 7A). Consider the two vortexes on the left side of the image (0 < *x* < 2); the flow pattern gathers the sedimenting capsules between these two vortexes, at the region where the fluid rises from the bottom wall. This region is close to the enzyme pattern because the updraft is driven by the low-density products of the enzymatic reaction. The same behavior occurs between the two vortexes on the right (2 < *x* < 4). Hence, this first stage of aggregation drives the capsules to roughly reproduce the shape of the enzyme pattern.

As the reagent becomes depleted from the capsules, the flow structure changes and only two vortexes remain in the *x*-*y* cross section (Fig. 7B). In this second stage of aggregation, the capsules are drawn away from the enzyme pattern toward the center of the domain. After all the reagent has been consumed, the fluid motion stops and the final configuration of the assembled capsules is a square colony with rounded edges; the assembly is shifted toward the center relative to the enzyme-coated pattern (Fig. 7C).

The final distribution of capsules can be tuned by adjusting the balance between the two stages of aggregation. As noted in the previous section, larger capsules tend to generate higher flow speeds and also move faster than small capsules near the bottom wall in the same flow field. Hence, we expect larger capsules to be more drawn to the center during the second aggregation stage than smaller capsules. Therefore, the areal distribution of small-sized capsules should more precisely mimic the shape of the enzyme-coated pattern than the distribution of larger capsules.

To demonstrate the differences in aggregation of small and large capsules, we compare two representative cases. Examples of colonies consisting of *N* = 1000 small capsules with radius *R* = 15 μm and generated by the circular, crankshaft, and square enzyme patterns are shown in the left column in Fig. 8. The correlation between the enzyme-coated patterns and the shapes of the capsule colonies indicates a dominant contribution of the first convective regime, which assembles capsules above the patterns with a minor subsequent shift of the colonies toward the center caused by the second regime. Because the aggregation process is a combination of convective mechanisms with reaction-diffusion processes, the resulting shape of the colony does not coincide precisely with the underlying enzyme pattern. In the present cases, sharp corners become rounded with a radius of curvature comparable to the size of the convective vortexes defined by the height of the microchannel.

Using the same enzyme patterns but with fewer (*N* = 500) capsules of larger radius (*R* = 40 μm), we achieve significantly different final configurations (right column in Fig. 8), indicating a crucial contribution of the second aggregation process, which assembles capsules in the center of the domain. Because the large capsules experience faster flow than the small ones, they form denser clusters that completely fill the hollow circle and square patterns. The contours of the colony induced by the crankshaft pattern are smoothened by the centripetal convective flow. Examples of the convective aggregation of the small- and large-radius capsules into the three different shapes are shown in movies S3 to S5 and S6 to S8.

The above results reveal that besides the design of the enzyme-coated pattern, there are other factors that control the shape of the capsule colony. The capsule size and their number determine the amount of fuel that controls the flow velocities and the relative contributions of the convective regimes aggregating the capsules either onto the pattern or in the center of the domain. The different response of the large and small capsules to the drag imposed by the convective flow might be used to segregate colonies of large and small capsules into ring-like structures of different diameters. Friction between the bottom wall and the capsules (ignored in the present model) could also suppress the effects of the second convective regime, which draws capsules from the enzyme-coated patch toward the center of the domain.

## CONCLUSIONS

We developed a computational approach to simulate an experimentally realizable system where an interlocking dynamic cascade of events leads to the spontaneous self-organization of microcapsules that were initially dispersed in an aqueous solution. The advantage of this mode of self-organization is that it does not rely on externally applied stimuli (such as electric or magnetic fields, or various gradients), but rather harnesses the inherent buoyancy force in fluids. The convective mechanism is triggered by a reagent that is released from the capsules. The reagent increases the local fluid density but is consumed by a reaction that takes place at enzyme-covered regions of the chamber. The release and consumption of the reagent generate variations in the solute concentration, which result in density variations that cause convective flows. These flows draw the capsules along the bottom surface toward the patch, leading to an accumulation of capsules. Higher numbers of capsules allow more fuel to be released into the system and result in stronger aggregation with smaller final intercapsule distance, δ (which scales with the capsule areal concentration *n* as ). The capsule aggregation is limited by the hard-core repulsion between neighbors and the formation of a close-packed monolayer.

The findings revealed a potential mechanism for promoting the spontaneous aggregation of protocells into colonies, where these early cells could communicate and, as a result, exhibited different types of complex collective behavior. Because protocells lacked complicated biochemical machinery, they interacted primarily through the release and uptake of molecules through a porous outer shell. For cells to communicate effectively via the diffusive exchange of molecules, they must come relatively close together. The protocells could have moved and assembled through the mechanism presented here because it relies solely on simple physical principles and does not require the use of a sophisticated internal structure or apparatus.

We found that for the system geometry considered here and the decomposition of hydrogen peroxide by catalase, the aggregation mechanism is capable of increasing the areal concentration of capsules by a factor of 10, potentially even bringing capsules into contact with their neighbors. Scaling considerations reveal that for a domain with spatial scale *L* and containing *N* capsules, the convective aggregation increases the initial uniform areal capsule density from *n*_{0} = *N*/*L*^{2} to *n*_{final} ~ 10 *n*_{0} and, thus, reduces the average intercapsule distance from to . Given that diffusive molecular exchange between neighboring cells is determined by the time scale τ^{diff} = δ^{2}/*D*, the aggregation can reduce communication time by a factor of 10, that is, . Reduction in communication time is a crucial factor in the transformation of a collection of separate individual cell-like objects (capsules) into a colony capable of collaboration and collective behavior.

The catalytic decomposition of hydrogen peroxide by catalase considered here is not the only catalyst/reagent combination that could fuel the assembly of the capsules. For example, the catalytic decomposition of hydrogen peroxide by metallic platinum will cause similar effects. Many other enzyme/reagent combinations, such as glucose oxidase/glucose and lipase/ester, have been shown to pump fluids (*13*). Furthermore, the pumping speed and aggregation behavior could be tuned through the appropriate choice of catalyst and the reagent release rate, which is a function of the size and shell permeability of the capsules. It is worth noting that the abundance of heterogeneous mineral-water interfaces present in the primordial ocean could have provided natural catalytic sites for the decomposition of a range of different reagents (*14*).

We also demonstrated that this convective mechanism can drive the self-organization of capsules into colonies with well-controlled shapes, which are dictated by the spatial patterning of the enzyme layer on the bottom of the microchamber and capsule properties, such as size and shell permeability. In this context, our results yield guidelines for promoting and regulating the self-assembly of microscopic objects in fluids. Moreover, the proposed convective capsule aggregation provides a means of manipulating suspended micrometer-sized cargos (for example, cells, viruses, and large molecules) within microliter fluid samples and directing their autonomous transport through microchannels.

## MATERIALS AND METHODS

Below, we describe the approaches we used to numerically solve the system of equations (Eqs. 2 to 6) with boundary conditions (Eqs. 10 to 13) governing the behavior of the system. First, however, it is convenient to rewrite the relevant variables in dimensionless form. Introducing new variables, we obtain(15)Dropping primes, we obtain the following dimensionless equations(16)(17)(18)(19)(20)where we have introduced the solutal and thermal Grashof dimensionless numbers, as well as the Schmidt number and Prandtl number . We also introduce to characterize capsule sedimentation in the absence of convection. The number and characterize the respective intercapsule and capsule-wall interactions. Finally, and the number describes the rate of reagent release out of capsule. The corresponding dimensionless boundary conditions on the domain walls are(21)(22)(23)(24)where the dimensionless number controls the maximum rate of chemical reaction and *K* determines the reagent saturation concentration.

The system of governing equations (Eqs. 16 to 20) with corresponding boundary conditions (Eqs. 21 to 24) was solved in the following way. First, a lattice Boltzmann algorithm (*25*) was applied to simulate the continuity and Navier-Stokes equations (Eqs. 16 and 17). Then, a second-order finite difference scheme was used to solve the diffusion and heat conduction equations (Eqs. 18 and 19). At each time step, the computed flow field **u** was used for advection terms in Eqs. 18 and 19. The updated concentration and temperature fields were then used for the buoyancy forces in Eq. 17 to advance the lattice Boltzmann scheme to the next time step. Finally, for known **u**, a first-order time integration was applied to solve the equation for capsule motion (Eq. 20).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/2/3/e1501835/DC1

Movie S1. Convective aggregation of *N* = 500 capsules with radius *R* = 15 μm and shell permeability *P* = 10^{−6} m/s on a square-shaped enzyme-coated patch.

Movie S2. Convective aggregation of *N* = 500 capsules with radius *R* = 20 μm and shell permeability *P* = 10^{−6} m/s on a square-shaped enzyme-coated pattern.

Movies S3 to S5. Convective aggregation of *N* = 1000 capsules with radius *R* = 15 μm and shell permeability *P* = 10^{−6} m/s on enzyme-coated patterns into the following shapes: square, circle, and crankshaft.

Movies S6 to S8. Convective aggregation of *N* = 500 capsules with radius *R* = 40 μm and shell permeability *P* = 10^{−6} m/s on enzyme-coated patterns into the following shapes: square, circle, and crankshaft.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Funding:**This work was supported by the Charles E. Kaufman Foundation (for the development of the theoretical model) and the Center for Bio-Inspired Energy Science, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under award no. DE-SC0000989 (for the computational studies).

**Author contributions:**A.C.B. and A.S. conceived the idea and supervised the work. O.E.S. and H.S. performed theoretical analysis and numerical simulations.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2016, The Authors