Molecular concentration field design using closed-form steady-state solutions

Abstract

Control over spatial concentration fields represents a fundamental challenge in designing synthetic biological systems and programmable soft materials. While nature creates morphogen gradients that orchestrate complex developmental processes, synthetic approaches have largely relied on empirical optimization and computationally intensive simulations. Here, we present an analytical framework for steady-state concentration fields generated by finite-sized localized sources in diffusion–degradation systems and derive closed-form solutions for one-, two-, and three-dimensional geometries. By expressing these solutions in dimensionless form, we show that gradient steepness and spatial structure are organized by the Thiele modulus, which captures the competition between diffusion and degradation length scales. The analysis reveals distinct design regimes: in degradation-dominated systems, gradient shape is governed by exponential decay and becomes dimension-independent, whereas in diffusion-dominated systems, gradient magnitude and extent follow dimension-dependent power-law scaling. Building on these results, we introduce a quantitative design strategy that uses threshold-based criteria to program concentration ranges by tuning physically accessible parameters, most directly the production rate, while holding transport and degradation properties fixed. Comparisons with numerical solutions and reported experimental systems demonstrate consistency with the predicted scaling behavior. Together, this work provides a generalizable and physically transparent framework for designing steady-state concentration fields in synthetic biological and soft matter systems, enabling predictive control of gradient-mediated organization without reliance on extensive numerical optimization.

---

Introduction

Chemical concentration fields in living cells and tissues guide biological processes across scales, from individual cellular behavior to complex multicellular patterning. Cells respond to diverse molecular gradients—for instance, bacteria follow nutrient concentrations,^1,2 immune cells migrate along chemokine gradients,³ neuronal axons grow in response to spatial distributions of guidance molecules,⁴ and morphogen distributions direct embryonic pattern formation.⁵ Similarly, synthetic systems have been developed in which concentration variations tune colloidal phase transitions,^6,7 control filament growth dynamics,^8,9 orchestrate droplet motility,^10,11 direct condensate formation,^12,13 regulate biochemical oscillations,^14,15 and establish spatially controlled patterning in synthetic biological circuits.^16–24 For realizing the full potential of these systems, a theoretical framework capable of predicting and controlling concentration fields is required for the rational design of chemical gradient landscapes tailored to applications ranging from tissue engineering to dissipative materials.²⁵

These chemical gradients arise from the interplay of localized reactions and molecular transport and are commonly described by reaction–diffusion dynamics, in which local production, degradation, and diffusion of molecules collectively shape spatial concentration profiles.^26,27 To capture and predict these behaviors, theoretical frameworks based on reaction-diffusion equations have been developed using a range of mathematical and computational approaches, each with distinct advantages and limitations. Partial differential equation (PDE) formulations with numerical solutions provide detailed spatiotemporal dynamics for complex multi-component systems, but the high computational cost often makes it challenging to explore large parameter spaces or apply them across many system variants.^28–34 To balance accuracy and computational efficiency, semi-analytical models that combine closed-form solutions with numerical evaluation steps have been developed.^35–37 These approaches are particularly useful for exploring parameter spaces and optimizing system designs in specific reaction-diffusion scenarios with the numerical evaluation of integrals or special functions, which can still impose computational overhead and limit their general applicability for different dimensionality, or complex boundary conditions.

Alternatively, analytical approaches have been developed through local accumulation time, which provide morphogen gradient dynamics before reaching steady state,³⁸ while their later work extended their model to multi-dimensional systems.³⁹ Other reaction–diffusion models describing diffusion and first-order morphogen degradation have been studied, in which steady-state solutions under idealized boundary conditions yield exponential concentration profiles characterized by a decay length 1/λ = (D/k)^0.5, where k is the degradation rate constant and D is the diffusion constant.^40–43 In these frameworks, the ratio λ/L—where L is the system size—emerges as a key dimensionless parameter governing pattern scaling and robustness, and closed-form solutions are available primarily for one-dimensional continuous domains with prescribed Dirichlet or Neumann boundary conditions. However, existing closed-form solutions do not directly address systems with finite-sized localized sources, nor do they provide a unified description across dimensionality for the predictive design of steady-state concentration gradients in systems such as dissipative materials that require continuous energy input to maintain far-from-equilibrium organization.

Here, we address this challenge by developing an analytically tractable framework for predicting steady-state concentration fields generated by localized-source reaction–diffusion systems. We derive closed-form steady-state solutions for one-, two-, and three-dimensional geometries in which continuous molecular production at finite-sized sources is balanced by diffusion and first-order degradation. By expressing these solutions in dimensionless form, we show that gradient steepness and spatial structure are organized by the Thiele modulus, which captures the competition between diffusion and degradation length scales and delineates diffusion- and degradation-dominated regimes. Building on these results, we introduce a transferable, threshold-based design strategy that leverages analytically derived scaling relationships to connect application-level requirements—such as target concentration and target active range—to experimentally tunable parameters, particularly the production rate, across dimensionalities. Together, this framework provides a quantitative and transferable foundation for designing biochemical gradients in synthetic materials, tissue engineering platforms, and autonomous soft matter systems that respond to long-term chemical cues.

Results and discussion

Analytical framework

To establish a tractable and physically grounded model of the concentration field, we consider a system composed of three essential components: a localized production source, diffusing molecules, and degradation enzymes or other catalysts (Fig. 1a). The source with diameter 2R is modeled as a highly porous region that is symmetric about its center in one-, two-, or three-dimensional space (Fig. 1b–d). Within the confined region of the source, molecules are homogeneously produced at a constant rate r_p, while simultaneously diffusing and undergoing degradation via a first-order reaction r_c with rate constant k via catalyzed degradation. The interplay of production, diffusion, and degradation generates a concentration gradient, described by the concentration field C(r,t), where r is the spatial coordinate and symmetric from the center of the source and t is time. Concentration fields can be calculated by solving piecewise partial differential equations in the separate regions. For the source region (0 < r < R), the system has production, degradation, and diffusion:

Fig. 1 Conceptual framework for analytical modeling of concentration fields generated through reaction-diffusion systems with a localized source and a uniform sink. (a) Schematic representation of the system, where discrete production sources (green region) continuously generate molecules (blue spheres) that diffuse outward while being consumed by catalysts (orange) that degrade the molecules, establishing a steady-state concentration gradient governed by the balance between production and degradation processes. (b)–(d) System configurations for analytical modeling, showing how the same principles can be applied in different dimensions: (b) a 1D system with a source of radius R, creating concentration profiles C(r) that decay along the channel length L, (c) a 2D system where a cylindrical source with radius R generates radially symmetric concentration fields C(r) where r is the spatial coordinate from the source, and (d) 3D system where a spherical source with a radius R create radially symmetric concentration fields C(r).

Beyond the source (r > R), the system has only diffusion and degradation without the production reaction:

At steady state, where the concentration field no longer varies with time (∂C/∂t = 0), four boundary conditions—symmetry, two continuity conditions, and a dilution condition (see SI for details)—allow us to derive exact solutions across 1D, 2D, and 3D systems (Table S1). The full derivation is provided in the SI.

To generalize the analytical solutions, we used a dimensionless framework defined by three quantities: C′ = C/(r_p/k), r′ = r/R, Φ = R(k/D)^0.5. Here, concentration C is scaled by the ratio of production to degradation rates, the spatial coordinate r is normalized by the source radius R, and the Thiele modulus Φ quantifies the relative strength of degradation versus diffusion.⁴⁴ Specifically, Φ characterizes the relative contributions of diffusion and degradation to concentration decay: when Φ > 1, degradation dominates, whereas when Φ < 1, diffusion is dominant. The resulting solutions with and without dimensionless quantities are summarized in Table 1 and Table S1. In 1D systems, the steady-state concentration profiles are well studied and are governed by hyperbolic functions within the source region and exponential decay outside the source: cosh(Φr′) within the source region r' < 1 and exp(−Φr′) in the outer region for r′ > 1.⁴⁵ However, few solutions in 2D and 3D have been reported. In 2D systems, the solutions involve modified Bessel functions of the first and second kind (I₁ and K₀), capturing the logarithmic behavior. For 3D systems, the solutions are described by sinh(Φr′)/r′ within the source and exp(−Φr′)/r′ outside the source, reflecting the radial decay of concentration in spherical symmetry. The analytical solutions closely match independent numerical simulations across 1D, 2D, and 3D system configurations and a range of parameters (Φ = 0.3, 1, 3), with total deviations below 10⁻⁶ (Fig. 2a–c). While this dimensionless formulation organizes gradient behavior into diffusion- and degradation-dominated regimes, practical gradient design requires additional tuning of physically accessible parameters such as the production rate, as discussed in Design strategy section.

Table 1 Analytical solutions for steady-state concentration fields across dimensional configurations. Closed-form expressions for dimensionless concentration C'(r′) in 1D, 2D, and 3D systems, presented inside (r′ < 1) and outside (r′ > 1) the source. The dimensionless variables are defined as C′ = C/(r_p/k), r′ = r/R, and Φ = R(k/D)^0.5, where r is the spatial coordinate from the center of the source, k is the degradation rate constant, r_p is the production rate, and D is the diffusion coefficient. For 2D systems, I₁ and K₀ denote modified Bessel functions of the first and second kind, respectively

	r′ < 1	r′ > 1
1D	C′(r′) = 1 − e^−Φcosh(Φr′)	C′(r′) = sinh(Φ)e^−Φr′
2D	C′(r′) = 1 − ΦK1(Φ)I0(Φr′)	C′(r′) = ΦI1(Φ)K0(Φr′)
3D

---

Fig. 2 Validation of analytical concentration field solutions in 1D, 2D, and 3D systems through comparison with numerical simulations, and their semi-logarithmic plots. Dimensionless concentration profiles C′ =C/(r_p/k) as a function of dimensionless radial distance r/R demonstrate the precision of our closed-form analytical solutions (solid lines) against independent numerical calculations³¹ (circles) for three distinct system dimensions in Fig. 1b–d: (a) 1D, (b) 2D, and (c) 3D systems, and their semi-logarithmic plots (d)–(f). The systematic variation of the dimensionless parameter Φ = R(k/D)^0.5 from 0.3 to 3 reveals how the characteristic length scale controls gradient steepness with larger Φ values producing sharper concentration drops near the source boundary. The dashed lines in d–f represent analytical extrapolations that extend from the far from the source to near the source. The shaded gray region indicates the source. For numerical calculation, we solved eqn (1) and (2) by formulating them as boundary-value problems and using MATLAB's bvp4c, a fourth-order collocation solver with adaptive mesh refinement.

Spatial decay behavior near and far from the source

We investigate the scaling principles that govern how concentration profiles decay as a function of r/R. To analyze this behavior, we used semi-logarithmic plots for 1D, 2D, and 3D systems across different parameter ranges (Φ = 0.3, 1, 3) (Fig. 2d–f). In the region far from the source (r/R ≫ 1), the profiles for 1D, 2D, and 3D systems each become linear in the semi-log plots. This behavior is consistent with the analytical solutions, which show that the dimensionless concentration derivative (dC'/dr) asymptotically converges to an exponential form proportional to ∼exp(−Φ) (Table 1). This behavior reflects an exponential decay, arising from diffusion and first-order degradation. Importantly, the linear relationship can be used for controlling the gradient steepness by manipulating Φ (= R(k/D)^0.5). Consistent with this interpretation, increasing either the source radius or the degradation rate produces steeper gradient profiles.

For the 1D system, the linear trend appears both near and far from the source boundaries, consistent with classical diffusion theory⁴⁶ (Fig. 2d). In contrast, the 2D and 3D systems exhibit nonlinear transitions near the source boundaries (r/R∼1), where the profiles deviate from linearity (Fig. 2e and f). As dimensionality increases, the available volume or area for diffusion grows, allowing molecules to disperse more rapidly away from the source—e.g., higher diffusion flux in 3D system compared to 1D system where molecules are restricted to a line. These deviations are more pronounced at low Φ, where diffusion dominates over degradation, producing a more gradual initial decay of concentration and corresponding nonlinear behavior near the source. Mathematically, this is reflected in the dimensionless concentration derivatives near the source: in 2D, dC′/dr ∼ exp(−Φ)/r′, while in 3D, dC′/dr ∼ exp(−Φ)/r', up to dimension-dependent prefactors. This trend illustrates how diffusion becomes increasingly sensitive to spatial dimensionality, leading to faster concentration decay near the source in higher-dimensional systems.

Shape of concentration gradients: steepness and magnitude

To characterize the steepness and magnitude of concentration gradients, the simple linear relationship valid in the far-field cannot be applied near the source, where 2D and 3D exhibit nonlinear transitions. To capture these features, we quantified gradient magnitude using the concentration at the source surface and gradient steepness using the half-width at half-maximum (HWHM). The surface concentration at r/R = 1 was chosen, rather than the center concentration, because the gradient outside the source functions as the relevant morphogen signal, and the maximum external concentration occurs at the surface.

The surface concentration was calculated as a function of Φ (Fig. 3a). In the diffusion -dominated regime (Φ < 1), (summarized in Table S2) exhibits distinct power-law scaling behaviors with Φ, where the slope of the log–log plot corresponds to the scaling exponent that characterizes how the gradient peak varies with Φ. Specifically, 1D systems show linear scaling (slope = 1), 2D systems display intermediate behavior (slope ≈ 1.5), and 3D system exhibits quadratic scaling (slope = 2). These different slopes indicate how the radius of the source and degradation rate constant as Φ = R(k/D)^0.5 becomes effective in higher dimensions. By contrast, in the degradation-dominated regime (Φ > 1), remains constant with respect to Φ.

Fig. 3 Dimensional scaling laws for concentration field optimization and design. (a) Surface concentration as a function of dimensionless parameter Φ = R(k/D)^0.5 demonstrates distinct scaling behaviors across dimensional configurations, with 1D system exhibiting linear scaling (slope = 1), 2D system showing intermediate behavior (slope = 1.5), and 3D system displaying quadratic scaling (slope = 2), reflecting how source-to-volume ratios fundamentally determine achievable surface concentrations. (b) Half-width at half-maximum (HWHM) normalized by source radius R reveals the spatial extent characteristics of concentration fields, where 1D system exhibits the strongest spatial confinement effects (slope = −1), 2D system shows moderate spatial expansion (slope = −0.5), and 3D system maintains constant width (slope = 0). HWHM/R values from experimental data were compared. Zadorin et al.⁴³ (△) used morphogen gradients in 1D microfluidic devices to direct particle organization. The others (Kim et al.⁴⁷ (○), Gines et al.⁴⁸ (◊), Karzbrun et al.¹⁴ (□)) built concentration fields in 2D microfluidic devices. Detailed values are summarized in Table S4.

The HWHM, normalized by the source radius R (summarized in Table S3) was defined as the distance from the center of the source at which the concentration reaches half of the surface concentration and was calculated as a function of Φ (Fig. 3b). In the diffusion-dominated regime (Φ < 1), the HWHM/R exhibits power-law scaling with Φ, where the slope indicates how strongly gradient steepness depends on Φ. The 1D system exhibits strongest dependence (slope = −1), indicating that gradients become sharper as Φ increases. The 2D system displays moderate dependence (slope = −0.5), whereas the three-dimensional system maintains a nearly constant HWHM/R ≈ 2. This result means that the concentration falls to half of the surface value at approximately twice the source radius, regardless of Φ (Fig. 3c). Thus, the dependence on Φ weakens with increasing dimensionality. In the degradation-dominated regime, HWHM/R is independent of Φ across 1D, 2D, and 3D systems.

The physical origin of the slope in the low Φ regime likely arises from the enhanced diffusive flux in higher dimensions, as discussed in the previous section, which enables molecules to disperse into larger available volume, leading to stronger dilution near the source and corresponding changes in scaling exponents. For quasi-one-dimensional or quasi-two-dimensional systems, where confinement in width or height is finite but not negligible, the scaling exponents are expected to interpolate between those of the lower- and higher-dimensional systems. For example, a quasi-2D system with finite height but sources confined to the bottom boundary (e.g., sedimented biomolecular condensates¹³) allows diffusion away from the substrate into the bulk, and is therefore better approximated by a hemispherical or effectively three-dimensional geometry, leading to slope values intermediate between the 2D and 3D limits.

For comparison with experimental observations, we compiled experimental values of HWHM/R, from previously reported morphogen and reaction–diffusion systems in 1D and 2D microfluidic configurations. These data points are overlaid with our analytical scaling predictions in Fig. 3b, and the corresponding physical parameters (diffusion coefficients, source sizes, degradation rates, and calculated Thiele moduli) are summarized in Table S4. As shown, the experimental measurements are consistent with the dimension-dependent scaling trends predicted by the analytical solutions, supporting the applicability of the theoretical framework to experimentally realizable reaction–diffusion systems.

Design strategy for concentration fields

Based on these insights, we propose a practical strategy for designing steady-state concentration fields using threshold-based criteria (Fig. 4a). In both synthetic reaction–diffusion systems and biological morphogen frameworks, spatial organization is often governed by whether local concentrations exceed a critical threshold, beyond which specific structures, patterns, or functional responses are activated. This threshold-based interpretation is widely employed in synthetic systems to direct material assembly and particle organization, as well as in biological systems to regulate tissue patterning and morphogenesis.^{5,26,42,43,49} Motivated by this common principle, we introduce a design guideline in which concentration fields are engineered to achieve a target threshold concentration C* over a prescribed spatial extent by tuning physical parameters.

Fig. 4 Workflow for designing steady-state concentration fields. (a) Schematic illustration of the design goal: tuning the production rate r_p to shift the spatial position at which a target threshold concentration C* is reached, thereby controlling the target active range without altering degradation kinetics or source geometry. Lower production rates yield shorter target active ranges, while higher production rates extend the range at which C* is attained. (b) The design process begins by identifying whether the target gradient lies in the degradation-dominated regime (HWHM < 2R) or the diffusion-dominated regime (HWHM > 2R), as exemplified by localized RNA transcription from DNA-immobilized hydrogels. (c) Target application requirements are then specified by defining a target active range and a threshold concentration, illustrated here by a target threshold C* of 1 nM and a target active range L = 100 µm (d) These targets are mapped back onto physical parameters to extend or tune the field range, showing how the effective production rate r_p can be adjusted to achieve the target active range spanning L, 2L, and 3L. (e) Representative steady-state concentration profiles in one-, two-, and three-dimensional geometries. Black, red, and blue curves correspond to profiles with characteristic ranges of L, 2L, and 3L, respectively, highlighting the role of dimensionality in shaping gradient decay.

The first design step is to determine whether the target concentration field is steep or extended by assessing whether the half-width at half-maximum (HWHM) is smaller or larger than twice the source radius (2R). When HWHM < 2R—corresponding to the degradation-dominated regime (Φ > 1)— the concentration profile decays steeply, with gradient shape governed by simple exponential decay (C′ ∼ e^−Φr′), and the surface concentration approaches a constant value (C′_surface ≈ 0.5) independent of system dimensionality. In contrast, when HWHM > 2R —the diffusion-dominated regime (Φ < 1)— the concentration field extends over longer distances, and both gradient steepness and magnitude follow dimension-dependent power-law scaling with Φⁿ, where the exponent reflects the effective dimensionality of molecular transport. Once the appropriate regime is determined, the physical parameters R, k, and D can be selected to achieve the desired gradient profile; for instance, operating in the diffusion-dominated regime is advantageous for generating long-range concentration fields in localized RNA transcription systems using DNA-immobilized objects, where diffusion coefficients, source sizes, and degradation rates are often constrained to narrow experimental ranges (Fig. 4b).^13,47

In the second step, the threshold concentration (C*) and target active range (L) are specified based on upstream biochemical or material constraints (Fig. 4c). For example, a threshold concentration of C* = 1 nM and a target active range of L = 100 µm are representative of systems such as protein-aptamer–based biosensors, where dissociation constants typically span sub-nanomolar to micromolar concentrations,⁵⁰ or linker-mediated particle assembly processes that operate in the nanomolar to micromolar range.^6,43 By tuning the production rate while holding degradation and transport parameters fixed, we determine production rates of 28, 150, and 850 nM min⁻¹ that yield a concentration of 1 nM at a distance of 100 µm for the 1D, 2D, and 3D cases, respectively.

Finally, to extend the target active range to longer distances, the concentration field can in principle be tuned by adjusting the production rate, degradation rate constant, or source radius. Among these parameters, we focus on the production rate as the primary control variable because it modulates the overall magnitude of the concentration field without altering gradient steepness, whereas changes in source size or degradation rate affect both magnitude and spatial decay (Tables S2 and S3). In practical implementations, the production rate can be directly controlled through experimentally accessible parameters; for example, in localized transcription reactions within hydrogels or biomolecular condensates, it can be tuned by varying the concentration of immobilized DNA templates or the activity of the transcription machinery.^13,47

By rearranging the analytical solutions in Table S1, we find that the production rate required to maintain a fixed threshold concentration at distance L scales as exp(λL), 1/K₀(λL), L·exp(λL), for 1, 2, 3D systems, respectively (Fig. 4d). These relationships provide direct guidelines for selecting production rates that extend the target active range, for example to 2L or 3L, without modifying degradation kinetics or source geometry. We note that while analogous design rules based on varying the degradation rate or source size are less analytically tractable, they can be obtained numerically. Using the production rates determined from Fig. 4d, the resulting concentration profiles in one-, two-, and three-dimensional systems exhibit threshold crossing at 1L, 2L or 3L, respectively, as shown in Fig. 4e.

Conclusions

In this work, we developed an analytically tractable framework for predicting and designing steady-state concentration fields generated by localized reaction–diffusion sources in one-, two-, and three-dimensional systems. By deriving closed-form solutions for finite-sized sources and expressing them in dimensionless form, we identified how gradient steepness and spatial structure are organized by the Thiele modulus, which delineates diffusion- and degradation-dominated regimes. These results clarify how dimensionality fundamentally shapes concentration decay and spatial extent, and they establish scaling relationships that are consistent with numerical simulations and reported experimental systems.^14,43,47,48

Beyond analytical description, this framework enables a practical, threshold-based strategy for concentration field design. Rather than relying on empirical tuning or extensive numerical optimization, target concentration ranges can be programmed by mapping application-specific requirements—such as threshold concentration and spatial reach—onto physically accessible parameters. In particular, we show that tuning the production rate provides a direct and experimentally viable means to extend or contract target active ranges while preserving gradient steepness and degradation kinetics, offering a robust design handle across dimensions.

Together, these results provide a quantitative and generalizable foundation for engineering reaction–diffusion concentration fields in synthetic and dissipative materials,^25,51,52 where precise spatial control under nonequilibrium conditions is essential. The framework is broadly applicable to systems ranging from programmable biosensors^31,53 and gradient-guided assembly to fuel-driven soft matter architectures that sustain autonomous spatial organization.^54–56 Future extensions incorporating time-dependent dynamics,^25,52 nonlinear reaction kinetics,⁵⁷ or spatial heterogeneity⁵¹ could further expand the applicability of this approach, positioning analytical gradient design as a key tool for the emerging field of programmable soft matter and synthetic morphogen systems, with direct relevance to tissue engineering and regenerative medicine.^58–62

Author contributions

Conceptualization, R. S. and D. W. K.; funding acquisition, R. S.; investigation, D. W. K.; project administration, R. S.; software, D. W. K.; resources, R. S.; supervision, R. S.; validation, D. W. K.; visualization, D. W. K.; writing – original draft, R. S., D. W. K. and A. G.; writing – review & editing, R. S., D. W. K. and A. G.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data and code associated with the publication are available at both https://github.com/DongWoo94/ConcentrationFieldModels and https://doi.org/10.7281/T1MLLJ33, Johns Hopkins Research Data Repository. Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5sm01121c.

Acknowledgements

We thank Marc Kamionkowski for insightful discussions. R. S. and A. G. acknowledge support from a Vannevar Bush Faculty Fellowship, ONR N00014-23-1-2868. D. W. K. and therefore the analytical and simulation work, were supported by the U.S. Department of Energy through DE-SC0010426.

References

1. J. H. Wheeler, K. R. Foster and W. M. Durham, Nat. Microbiol., 2024, 9, 2308–2322.
2. V. Sourjik and N. S. Wingreen, Curr. Opin. Cell Biol., 2012, 24, 262–268.
3. J. W. Griffith, C. L. Sokol and A. D. Luster, Annu. Rev. Immunol., 2014, 32, 659–702.
4. A. L. Kolodkin and M. Tessier-Lavigne, Cold Spring Harbor Perspect. Biol., 2011, 3, a001727.
5. K. W. Rogers and A. F. Schier, Annu. Rev. Cell Dev. Biol., 2011, 27, 377–407.
6. J. Lowensohn, B. Oyarzún, G. Narváez Paliza, B. M. Mognetti and W. B. Rogers, Phys. Rev. X, 2019, 9, 041054.
7. W. B. Rogers and V. N. Manoharan, Science, 2015, 347, 639–642.
8. R. Schulman, B. Yurke and E. Winfree, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 6405–6410.
9. D. Krishnamurthy and M. Prakash, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2304981120.
10. Y. Zhang, A. McMullen, L.-L. Pontani, X. He, R. Sha, N. C. Seeman, J. Brujic and P. M. Chaikin, Nat. Commun., 2017, 8, 21.
11. N. J. Cira, A. Benusiglio and M. Prakash, Nature, 2015, 519, 446–450.
12. M. Morita, T. Sakamoto, S. I. M. Nomura, S. Murata, M. Yanagisawa and M. Takinoue, Adv. Mater. Interfaces, 2024, 11, 2300898.
13. E. Kengmana, E. Ornelas-Gatdula, K.-L. Chen and R. Schulman, J. Am. Chem. Soc., 2024, 146, 32942–32952.
14. E. Karzbrun, A. M. Tayar, V. Noireaux and R. H. Bar-Ziv, Science, 2014, 345, 829–832.
15. M. Weitz, J. Kim, K. Kapsner, E. Winfree, E. Franco and F. C. Simmel, Nat. Chem., 2014, 6, 295–302.
16. Y. Fang, V. V. Yashin, S. J. Dickerson and A. C. Balazs, J. Appl. Phys., 2018, 123, 215107.
17. A. Manfrin, Y. Tabata, E. R. Paquet, A. R. Vuaridel, F. R. Rivest, F. Naef and M. P. Lutolf, Nat. Methods, 2019, 16, 640–648.
18. S. Kondo and T. Miura, Science, 2010, 329, 1616–1620.
19. B. A. Grzybowski and W. T. Huck, Nat. Nanotechnol., 2016, 11, 585–592.
20. S. Biswas, V. V. Yashin and A. C. Balazs, Soft Matter, 2018, 14, 3361–3371.
21. M. Moradi, O. E. Shklyaev and A. C. Balazs, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2319777121.
22. K. Abe, I. Kawamata, M. N. Shin-ichiro and S. Murata, Mol. Syst. Des. Eng., 2019, 4, 639–643.
23. I. Kawamata, T. Hosoya, F. Takabatake, K. Sugawara, S.-I. Nomura, T. Isokawa, F. Peper, M. Hagiya and S. Murata, 2016.
24. I. Kawamata, S. Yoshizawa, F. Takabatake, K. Sugawara and S. Murata, 2016.
25. E. Del Grosso, E. Franco, L. J. Prins and F. Ricci, Nat. Chem., 2022, 14, 600–613.
26. T. D. Banerjee and A. Monteiro, Commun. Biol., 2025, 9, 87.
27. H. Connahs, S. Tlili, J. van Creij, T. Y. Loo, T. D. Banerjee, T. E. Saunders and A. Monteiro, Development, 2019, 146, dev169367.
28. A. N. Ford Versypt, P. D. Arendt, D. W. Pack and R. D. Braatz, PLoS One, 2015, 10, e0135506.
29. M. S. Johnson, J. N. Mueller, C. Daniels, H. N. Najm and J. Zádor, J. Phys. Chem. A, 2024, 128, 3685–3702.
30. D. Scalise and R. Schulman, Technology, 2014, 2, 55–66.
31. P. J. Dorsey, D. Scalise and R. Schulman, R. Soc. Open Sci., 2022, 9, 220200.
32. D. Scalise and R. Schulman, Nat. Comput., 2016, 15, 197–214.
33. K. Abe, S. Murata and I. Kawamata, Soft Matter, 2021, 17, 6160–6167.
34. S. Chen and G. Seelig, Soft Matter, 2020, 16, 3555–3563.
35. K. Al Noufaey, Results Phys., 2018, 9, 609–614.
36. H. Y. Alfifi, ANZIAM J., 2017, 59, 167–182.
37. H. Alfifi, T. R. Marchant and M. Nelson, J. Math. Chem., 2016, 54, 1632–1657.
38. P. V. Gordon, C. Sample, A. M. Berezhkovskii, C. B. Muratov and S. Y. Shvartsman, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 6157–6162.
39. P. V. Gordon, C. B. Muratov and S. Y. Shvartsman, J. Chem. Phys., 2013, 138, 104121.
40. B. Vasiev, J. Theor. Biol., 2016, 404, 109–119.
41. H.-W. Kang, L. Zheng and H. G. Othmer, Interface Focus, 2012, 2, 465–486.
42. P. C. Bressloff, Phys. Biol., 2022, 19, 066005.
43. A. S. Zadorin, Y. Rondelez, G. Gines, V. Dilhas, G. Urtel, A. Zambrano, J.-C. Galas and A. Estevez-Torres, Nat. Chem., 2017, 9, 990–996.
44. E. W. Thiele, Ind. Eng. Chem., 1939, 31, 916–920.
45. B. Cushman-Roisin, Thayer School of Engineering Dartmouth College, University Lecture, 2012.
46. P. J. Dorsey, M. Rubanov, W. Wang and R. Schulman, ACS Macro Lett., 2019, 8, 1133–1140.
47. D. W. Kim, M. Rubanov, A. Grinthal, P. Moerman and R. Schulman, Matter, 2025, 8(7), 102208.
48. G. Gines, A. Zadorin, J.-C. Galas, T. Fujii, A. Estevez-Torres and Y. Rondelez, Nat. Nanotechnol., 2017, 12, 351–359.
49. H. Inomata, Dev., Growth Differ., 2017, 59, 41–51.
50. H. Lee, T. Xie, B. Kang, X. Yu, S. W. Schaffter and R. Schulman, Nat. Commun., 2024, 15, 7973.
51. B. Dúzs and I. Szalai, Commun. Chem., 2020, 3, 168.
52. J. Boekhoven, W. E. Hendriksen, G. J. Koper, R. Eelkema and J. H. van Esch, Science, 2015, 349, 1075–1079.
53. S. Jia, S. C. Phua, Y. Nihongaki, Y. Li, M. Pacella, Y. Li, A. M. Mohammed, S. Sun, T. Inoue and R. Schulman, Nat. Commun., 2021, 12, 5729.
54. M. Lovrak, W. E. Hendriksen, M. T. Kreutzer, V. van Steijn, R. Eelkema and J. H. van Esch, Soft Matter, 2019, 15, 4276–4283.
55. J. Zenk, D. Scalise, K. Wang, P. Dorsey, J. Fern, A. Cruz and R. Schulman, RSC Adv., 2017, 7, 18032–18040.
56. P. J. Dorsey, D. Scalise and R. Schulman, Angew. Chem., 2021, 133, 342–348.
57. K. L. Narayanan, V. Meena, L. Rajendran, J. Gao and S. P. Subbiah, Appl. Comput. Math., 2017, 6, 143–160.
58. Z. Zeravcic, V. N. Manoharan and M. P. Brenner, Rev. Mod. Phys., 2017, 89, 031001.
59. Z. Zeravcic and M. P. Brenner, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 4342–4347.
60. E. W. Gehrels, W. B. Rogers, Z. Zeravcic and V. N. Manoharan, ACS Nano, 2022, 16, 9195–9202.
61. N. Tompkins, N. Li, C. Girabawe, M. Heymann, G. B. Ermentrout, I. R. Epstein and S. Fraden, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 4397–4402.
62. I. Hunter, M. M. Norton, B. Chen, C. Simonetti, M. E. Moustaka, J. Touboul and S. Fraden, Phys. Rev. E, 2022, 105, 024310.

---