,disordered microstructure ,generative model ,machine learning ,multiscale optimization

1 Introduction↩︎

Mechanical metamaterials derive their properties mainly from the geometry of their microstructures [1]–[3]. They can achieve mechanical behaviors inaccessible to conventional materials [4]–[6], such as negative Poisson’s ratio [7], highly tunable energy absorption [8], [9], and superior impact resistance [10]. With their excellent performance and tunability, mechanical metamaterials have transformed applications ranging from lightweight structures to soft robotics [11], acoustic control [12], [13], and biomedical devices [14]. Most existing mechanical metamaterials assume ordered microstructures, featuring deterministic and regular arrangements of predefined building blocks [1], [2], [15], [16]. This order in assembly makes the building-block microstructures and their multiscale systems easy to parameterize and design [17], [18]. However, it also constrains the design space, limiting access to heterogeneous and more morphologically diverse microstructures [19]. Natural materials, such as trabecular bone [20], sea sponge [21], coral [22], and nacre [23], are instead built on disordered microstructures [24], featuring randomness, irregularity, and interconnectivity. Such structural disorder helps suppress deformation localization [25], deflect crack propagation [8], [26], [27], redistribute stress [28], and dissipate impact energy [29], leading to greater robustness and performance unavailable in ordered structures. These unique advantages reveal only the tip of a much broader design space and performance envelope hidden within disorder [19], [30], [31].

Despite their superior performance, disordered microstructures are more difficult to design than ordered ones, due to their inherent randomness and irregular geometries. A growing number of computational methods have emerged to address this challenge [19], [32]–[35]. One class of approaches generates disordered microstructures through prescribed rules or basis functions, including Voronoi tessellations that partition the domain around randomly placed seed points [36], Gaussian random fields or trigonometric basis functions to create spinodal-like bicontinuous microstructures [37]–[42], and wave function collapse methods that tile local pattern fragments to assemble globally consistent textures [10], [28], [43]. These rule-based methods address the parameterization challenge of disordered structures, yet each is confined to a narrow morphological class constrained by its construction rule, with no straightforward way to go beyond, transition between, or mix different microstructure types. A substantial portion of the vast design space of disordered materials thus remains untapped. Moreover, they require human-specified rules, which are nontrivial to obtain for the diverse, highly complex microstructures found in nature. However, these microstructures are precisely the ones that have great potential for exceptional properties and niche applications.

Deep generative models offer an alternative. Variational autoencoders [44], [45], generative adversarial networks [46], [47], and diffusion models [48], [49] learn the statistical distribution directly from disordered samples to reproduce microstructures without explicitly specified rules, which is more flexible than rule-based methods. However, this capability introduces several new limitations. First, training these models typically requires a large amount of microstructure samples, which are difficult to acquire [50]. Second, their ability to generalize beyond the training distribution is limited, as they struggle to generate morphologies that are underrepresented in the training data, or truly novel ones unseen in their training data [46]. Third, the trained models are also strictly tied to the domains and discretizations on which they were trained (oftentimes a square or cube with a uniform grid) and transfer poorly to other settings, typically requiring complete retraining.

In summary, rule-based methods can transfer across domains and discretizations but are restricted to manually prescribed rules, while deep generative models directly learn from data but require large datasets and remain tied to a fixed domain and discretization. Therefore, unlocking the full potential of disordered metamaterials requires a generative framework that can learn from data while remains data-efficient and domain-adaptable [1], [19].

Nature, the very thing that inspired the exploration of disordered metamaterials, provides a clue to address this gap. Many disordered microstructures we see in nature are formed by self-organization, where individual cells or chemicals respond only to their neighbors, yet global patterns emerge from these local interactions without any central blueprint. [51], [52]. This self-organizing principle creates microstructures whose diversity and complexity far exceed those of our engineered systems, adapting freely across domains and scales in contrast to our brittle models that cannot generalize. For example, leaf veins form through local transport and reinforcement between neighboring cells, yet these local interactions collectively produce a connected vascular network across the entire leaf with varying shapes [53]. This principle inspires us to propose a fundamentally different design paradigm built on Neural Cellular Automata  [54]–[56], which embeds such self-organizing dynamics to generate disordered microstructures (Fig. 1). Instead of learning the static microstructures themselves, as existing AI-driven methods do, we aim to find an underlying local interaction rule that dynamically grows them. And instead of constructing such a rule by hand, as in existing rule-based methods, we use neural networks to learn it from samples. In this way, the framework combines the locality of rule-based methods with the expressive power of neural networks. The former makes it data-efficient and generalizable to unseen domains, while the latter allows it to accommodate complex microstructures and different types. With it, we can grow disordered microstructures in a controllable way across scales, domains, and discretizations, just as nature grows its microstructures, but now in silico.

Specifically, we apply this self-organizing principle at three progressively increasing levels of complexity. We first focus on generating homogeneous microstructures (Fig. 1 (a)). An NCA model, trained from a single template, learns a dynamic growth process in which each cell iteratively extracts neighboring information through local perception and maps it to a state update through a shared network, growing microstructures that match those of the template. We further introduce post-training control mechanisms into the learned dynamics to manipulate orientation, anisotropy, characteristic scale, and directional thickness of generated microstructure. By combining these controls, we expand the mechanical property space far beyond that of the training template [57].

We then extend the framework to heterogeneous microstructures across different domains and discretizations (Fig. 1 (b)). It produces spatially varying microstructures with smooth morphological transitions between different growth controls and microstructure types. By adapting the local perception in the self-organizing rule, we transfer the learned rule to irregular domains and arbitrary discretizations different from the training setting.

Finally, we couple the growth rules with topology optimization to enable multiscale design with optimized system performance(Fig. 1 (b)). The optimizer assigns spatially varying microstructure properties, which are transformed into spatially varying growth control parameters, akin to growth cues in natural materials. This allows us to dynamically grow the entire multiscale structure using the same learned NCA while following spatially varying growth cues. We use mechanical cloaking design as an extreme case to stress-test the method: it produces smoothly varying, disordered multiscale structures that accommodate local property requirements and deliver excellent cloaking behavior.

Together, this study makes the following contributions:

• We establish a lightweight one-shot generative framework that learns self-organizing rules from a single sample and grows disordered microstructures across different domains and scales.

• By treating the NCA as a generalized PDE system, we achieve flexible control of both morphologies and mechanical properties beyond the training template without new data collection or retraining.

• The learned self-organizing rules transfer across irregular domains and discretizations different from the training setting. Their asynchronous growth further enables scalable generation in large domains, as well as self-repair.

• We integrate the NCA into a multiscale optimization pipeline, generating heterogeneous multiscale structures that meet location-specific property requirements with smooth transitions between them.

The remainder of this paper is organized as follows. Section 2 introduces the NCA model, its one-shot training strategy, and post-training steering mechanisms for the generation of homogeneous disordered microstructures. Section 3 extends the framework to generate heterogeneous microstructures and adapt to a 3D surface mesh. Section 4 presents the multiscale optimization pipeline and demonstrates its application to the mechanical cloaking problem. Section 5 concludes the paper with a discussion of limitations and future directions.

2 One-shot learning for control growth of homogeneous microstructures↩︎

In this section, we establish the NCA-based generative model for disordered microstructure design. We first establish the NCA model and demonstrate its flexibility across distinct disordered microstructure types. We then implement four steering mechanisms to control both morphologies and mechanical properties over the generated microstructure without retraining the model. Finally, we explore how these controls broaden the mechanical property space of the generated microstructures, laying the foundation for the multiscale design in Sec. 4.

2.1 One-shot learning with Neural Cellular Automata↩︎

2.1.1 NCA models↩︎

Neural Cellular Automata (NCA) provide a natural framework for learning self-organizing dynamics to grow microstructures from local interactions. We consider each discretized component in a microstructure as a cell, which can be a pixel, voxel, element, or node. In an NCA, each cell iteratively updates its state based on information perceived from the states of its neighborhood [54]. Instead of manually prescribing the update rule, the model represents the update rule with a neural network shared across all cells, whose parameters are learned directly from data. By repeatedly applying the learned update rule across space and time, complex global morphologies emerge from purely local interactions, making NCA well suited for modeling disordered microstructures with diverse and irregular morphologies.

Specifically, each cell in the growth domain carries a state vector \(\mathbf{u} \in \mathbb{R}^{C}\), where the first channel stores the solid-void phase and the remaining \(C-1\) channels are hidden channels that allow cells to encode and exchange latent information, providing the memory and context needed for coherent global patterns to emerge, as illustrated in the local region and state channels in Fig. 2. In each update step, the state of every cell \(\mathbf{u}^{(t)}\) is updated according to \[\mathbf{u}^{(t+1)} = \mathbf{u}^{(t)} + \delta^{(t)} \cdot F_{\theta}\!\left(\mathit{percept}(\mathbf{u}^{(t)})\right), \label{eq:stochastic95update}\tag{1}\] where \(\mathit{percept}\) is a fixed local perception operator to collect neighboring information, \(F_{\theta}\) is a trainable neural network shared across cells that maps this neighboring information to an incremental state update, and \(\delta^{(t)}\) is a stochastic update mask. The stochastic update mask \(\delta^{(t)}\) in Eq. 1 is sampled independently for each cell from a Bernoulli distribution with \(p = 0.5\) at every update step. This stochastic masking breaks spatial symmetry to excite emergent dynamics. It also functions as a form of spatial dropout during training, improving the robustness and generalization of the learned rule to large, irregular, and partially initialized domains, which will be exploited in the asynchronous growth capability introduced in Section 3.3.

While the model is general, we first focus on microstructures represented on pixelated regular grids for ease of illustration and extend it to more general representations in later sections. Given a pixelated grid, the local perception operator extracts neighborhood information in each cell through a set of fixed \(3\times3\) depthwise convolution kernels: Sobel filters \(\mathbf{K}_{x}\) and \(\mathbf{K}_{y}\) approximating the spatial gradients \(\nabla_x\mathbf{u}\) and \(\nabla_y\mathbf{u}\), and a Laplacian filter \(\mathbf{K}_{\mathrm{lap}}\) approximating \(\nabla^2\mathbf{u}\), as shown in Fig. 2. Their outputs are concatenated with the cell’s own state to form the perception vector \[\mathit{percept} = \bigl[\,\mathbf{u},\; \mathbf{K}_{x} \ast \mathbf{u},\; \mathbf{K}_{y} \ast \mathbf{u},\; \mathbf{K}_{\mathrm{lap}} \ast \mathbf{u}\,\bigr] \in \mathbb{R}^{4C}. \label{eq:perception}\tag{2}\]

Note that other sets of operators can be chosen for more general settings. We then pass the perception vectors to the shared neural network \(F_{\theta}\), which outputs a residual state update \(\Delta\mathbf{u} = F_{\theta}(\mathit{percept})\). The same network is shared across all cells, so that each cell executes the same local update rule on different neighbors, keeping the number of trainable parameters small and, crucially, independent of the grid size. A key insight of this study is that the resulting model can be interpreted as a generalized partial differential equation (PDE), allowing us to leverage rigorous operator theory and a rich set of analytical tools to improve controllability and generalization. Specifically, the update rule can be viewed as a forward Euler discretization of the following PDE: \[\frac{\partial \mathbf{u}}{\partial t} = \mathcal{F}_{\theta}\!\left( \mathbf{u},\; \nabla \mathbf{u},\; \nabla^{2} \mathbf{u} \right), \label{eq:pde}\tag{3}\] where \(\mathcal{F}_{\theta}\) is a nonlinear evolution function acting on the local PDE operators, approximated by \(F_{\theta}\) on locally perceived information in Eq. 2 . This PDE perspective unifies the NCA with classical self-organizing systems such as reaction-diffusion models [51], [52]. In both cases, the state \(\mathbf{u}\) evolves through purely local interactions represented by the local perception operator. Both also share the same emergent character: coherent global patterns arise from the cumulative effects of repeated local interactions across the domain. The key distinction is that classical PDE systems rely on manually prescribed rules for \(\mathcal{F}\), which confines them to a narrow class of morphologies, whereas the NCA learns \(\mathcal{F}_{\theta}\) directly from data, replacing hand-crafted rules with a trainable approximation that can capture the complex local dynamics to produce any target morphology. We can also view classical PDE systems as special cases of NCA with a prescribed (e.g., polynomial) form of the evolution operator \(\mathcal{F}_{\theta}\). Because the learned rule relies only on local interaction and does not depend on the resolution, representation, or geometry of the training domain, it can be applied to any target resolution or domain with the same network parameters. This reflects a fundamental distinction from conventional deep generative models: rather than learning a static mapping from the latent distribution to an entire structure, the NCA learns the evolution dynamics that grow the microstructure progressively, with the global pattern emerging from cumulative local interactions over time.

2.1.2 Model training↩︎

Instead of requiring a large dataset of microstructure samples, the NCA update network is trained in a one-shot manner from a single template. The goal of the training is not to replicate each individual pixel in the template, but to learn the local update rule that can grow statistically equivalent morphology [58], [59], which is more meaningful for disordered microstructures as we are learning the universal growth rules rather than a single realization of the template. Therefore, we use the style loss in the training, which measures similarity through feature distributions across multiple spatial scales, allowing positional variation while enforcing statistical similarity. We adopt a frozen ImageNet-pretrained VGG-16 as a fixed feature extractor and apply a sliced optimal transport loss to compute \(\mathcal{L}_{\mathrm{style}}\) by matching the feature distributions of the generated microstructure to those of the template across multiple intermediate layers without enforcing any spatial correspondence [60]. The total training loss combines the style loss \(\mathcal{L}_{\mathrm{style}}\) with a state overflow regularization term: \[\mathcal{L} = \mathcal{L}_{\mathrm{style}} + \sum_{x,y,c} \bigl| u_{x,y,c} - \mathrm{clamp}(u_{x,y,c},\,-1,\,1) \bigr|, \label{eq:total95loss}\tag{4}\] where \({clamp}(\cdot,\,-1,\,1)\) limits each state value to the range \([-1,1]\). The second regularization term therefore penalizes any state value that exceeds \([-1,1]\) to keep the state values bounded and stabilize training over long update sequences. The pretrained VGG-16 parameters are frozen during the training, which is only used for loss calculation. Parameters \(\theta\) in the NCA update network \(F_{\theta}\) are the only trainable parameters. They are optimized with the Adam optimizer with a stepwise decaying learning rate, with the key hyperparameters summarized in Table 2.

With these setups, we validate the framework by training it on a leaf vein microstructure, with the growth process illustrated in Fig. 2 (see Vid. [vid:nca growth] for the whole process). Starting from an empty seed state, scattered local fragments first emerge through the learned update rule and gradually connect into a coherent network, continually refining until the morphology stabilizes. The resulting microstructures closely mimic the template. To provide further quantitative validation, we apply two-point correlation and lineal path functions [61] to quantify the morphological features between the template and the generated microstructures in Fig. 3. The two-point correlation function measures the probability that two points separated by a given distance both lie in the solid phase, capturing the volume fraction and characteristic length scale, while the lineal path function measures the probability that an entire line segment of a given length lies within the solid phase, characterizing the connectivity of the solid network. The close agreement between the generated microstructures and the template across both metrics confirms that the learned update rule has successfully reproduced the key morphological features and underlying distribution of the template, not merely a single pixel-level realization.

2.1.3 Validation on multiple types of disordered microstructure↩︎

To demonstrate the flexibility of our framework, we selected multiple different types of disordered microstructures with distinct morphological features (Fig. 4): homogeneous and hierarchical network topologies, as in the leaf vein and vascular tissue; porous architectures with rounded strut junctions, as in the trabecular bone; dendritic branching systems, as in the biofilm; sparse fiber networks characterized by long-range connections and local conjunctions, as in the spider silk; and granular microstructures with sharp interfacial boundaries, as in the metal fracture surface. We trained a separate model for each template using the same architecture and hyperparameters.

The comparison of generated microstructures and their corresponding templates is shown in Fig. 4 (a). Each of the generated microstructures maintains the key morphological features in those of the corresponding template. The successful generation of these distinct disordered microstructures demonstrates the flexibility and robustness of our proposed generative framework. The homogenized elastic surfaces of these microstructures are shown in Fig. 4 (b). The six generated microstructures exhibit a wide range of mechanical responses from nearly isotropic to strongly anisotropic, demonstrating that the morphological diversity of the microstructures directly translates into a broad mechanical property space. This is enabled by the flexibility of neural networks over predefined rule-based methods.

2.2 Post-training microstructure control↩︎

The model structure and training strategy introduced in the previous section enable the NCA to learn a local update rule that reproduces the morphology of a given microstructure. In this section, we aim to extend this generative capability to go beyond reproduction to create new microstructures with controllable variation in the morphology and corresponding mechanical properties without retraining the model. As noted earlier, the learned self-organizing rule in NCA uses local operators in the same way that PDEs use them to describe pattern evolution. In PDE and computational geometry theory, there is a well-established analysis of how changes in these local operators alter the behavior of the PDE evolution. We can therefore transfer this analysis to the learned NCA model, steering the growth process to control the generated morphology (see Vid. [vid:individual control]), by modifying these local perception operators while keeping the learned update network fixed. Based on this principle, we introduce a set of post-training steering mechanisms to manipulate orientation, anisotropy, characteristic scale, and directional thickness, expanding the morphological diversity and mechanical property space.

2.2.1 Orientation control through gradient rotation↩︎

To control the orientation of growth, we apply a rotation matrix \(\mathbf{R}(\theta)\) to the spatial gradients in the local perception operation, rotating them by a target angle \(\theta\): \[\begin{pmatrix} \widetilde{\nabla}_x\,\mathbf{u}^{(t)} \\ \widetilde{\nabla}_y\,\mathbf{u}^{(t)} \end{pmatrix} = \mathbf{R}(\theta) \begin{pmatrix} \nabla_x\,\mathbf{u}^{(t)} \\ \nabla_y\,\mathbf{u}^{(t)} \end{pmatrix} =\begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix}\begin{pmatrix} \nabla_x\,\mathbf{u}^{(t)} \\ \nabla_y\,\mathbf{u}^{(t)} \end{pmatrix}, \label{eq:rotated95gradients}\tag{5}\] which can be written compactly as \(\widetilde{\nabla} = \mathbf{R}(\theta)\nabla\). As illustrated in Fig. 5, these steered gradients replace the original gradients as inputs for the steered update rule \[\mathbf{u}^{(t+1)} = \mathbf{u}^{(t)} + \delta^{(t)} \cdot F_{\theta}\!\left(\mathbf{u}^{(t)},\; \widetilde{\nabla}_x\,\mathbf{u}^{(t)},\; \widetilde{\nabla}_y\,\mathbf{u}^{(t)},\; \nabla^2\mathbf{u}^{(t)}\right), \label{eq:orientation95update}\tag{6}\]

which provides an orientation-steering mechanism for the self-organizing dynamics, producing microstructures with orientation closely following the prescribed angle (Fig. 6). The resulting effective stiffness also rotates accordingly. This gives direct control over both morphology and mechanical anisotropy. Notably, the steering reorients the principal stiffness axes without changing the degree of anisotropy: the elastic surfaces retain nearly the same shape and size across all steering angles.

2.2.2 Anisotropy and scale control↩︎

Although orientation steering rotates the generated microstructure and its elastic surface, it cannot control the degree of anisotropy. To further control the anisotropy and feature scale, we reshape the local coordinate system in which the NCA extracts neighborhood information using a unified Riemannian metric framework [62]. The steered update rule becomes \[\mathbf{u}^{(t+1)} = \mathbf{u}^{(t)} + \delta^{(t)} \cdot F_{\theta}\!\left( \mathbf{u}^{(t)},\; \mu^{-1/2}\mathbf{g}^{-1}\nabla\mathbf{u}^{(t)},\; \widetilde{\nabla}^2\mathbf{u}^{(t)}\right), \label{eq:anisotropy95scale95update}\tag{7}\] where the Riemannian metric tensor \(\mathbf{g}\) and the anisotropic Laplacian \(\widetilde{\nabla}^2\) are jointly parameterized by \(s_1\), \(s_2\), and \(\mu\): \[\mathbf{g}(s_1, s_2) = \exp(\mathbf{S}), \qquad \mathbf{S} = \begin{pmatrix} s_1 & s_2 \\ s_2 & -s_1 \end{pmatrix}, \qquad \widetilde{\nabla}^2\mathbf{u}^{(t)} = \sqrt{|\mu\mathbf{I}|}^{-1}\,\nabla \cdot \bigl(\mathbf{g}^{-1}\sqrt{|\mu\mathbf{I}|}\, \mu^{-1/2}\mathbf{g}^{-1}\nabla\mathbf{u}^{(t)}\bigr). \label{eq:riemannian95metric}\tag{8}\] Intuitively, the axial stretching parameter \(s_1\) controls the elongation ratio between the two local coordinate axes in calculating the local perception operation, steering the NCA to grow features that are directionally dominant along the horizontal or vertical direction. Similarly, the shearing parameter \(s_2\) introduces diagonal distortion of the local coordinate system, controlling the dominant direction between the off-diagonal and diagonal orientations. The scale parameter \(\mu\) uniformly changes the length scale of both local coordinate axes in the perception operation, controlling the size of morphological features. When \(\mu = 1\), this Riemannian control reduces to pure anisotropy steering, controlled by \(s_1\) and \(s_2\) alone. When \(s_1 = s_2 = 0\), the metric \(\mathbf{g}\) becomes the identity, and the Riemannian control reduces to pure scale steering. The scale and anisotropy controls are therefore naturally unified in this Riemannian metric and can be applied together or separately.

We use the leaf microstructures to illustrate this control mechanism. As shown in Fig. 8, \(s_1\) and \(s_2\) provide complementary control over axial and diagonal anisotropy of the generated microstructures, as intended. With increasing \(s_1\), the generated microstructure continuously transitions from vertically dominated to horizontally dominated, and the principal direction of stiffness changes correspondingly from the vertical to the horizontal axis. As \(s_2\) increases, the microstructure transitions from the off-diagonal to the diagonal direction, with the principal stiffness direction rotating accordingly. Together, \(s_1\) and \(s_2\) span the full directional anisotropy space through two control axes, producing continuous changes in both morphology and elastic response. Meanwhile, as shown in Fig. 9, increasing \(\mu\) progressively refines the feature scale, while the degree of stiffness anisotropy and its magnitude remain largely unchanged.

2.2.3 Directional thickness control↩︎

The three steering mechanisms introduced above operate on the perception stage and control the morphological features of the generated microstructure, but do not directly regulate the solid volume fraction or feature thickness that governs the overall stiffness level. We therefore implement directional thickness control by making a few additional growth steps in the binarized solid-void field \(\mathbf{u}\), as shown in Fig. 10. At each evolution step, the solid boundary is detected via the same Sobel filters as in the perception stage and updated as \[\tilde{\mathbf{u}} = \mathbf{u} + \bigl(s_x \left| \nabla_y \mathbf{u} \right| + s_y \left| \nabla_x \mathbf{u} \right|\bigr). \label{eq:grassfire95update}\tag{9}\] The absolute gradient magnitudes confine the update strictly to the solid-void interface, leaving the interior of each phase unaffected. And \(s_x\) and \(s_y\) are signed speed parameters controlling the boundary propagation in the horizontal and vertical directions, respectively. Positive \(s_x\) or \(s_y\) drives outward propagation of the solid boundary, dilating the solid phase and increasing the thickness, while negative values drive inward erosion. The two parameters act independently, enabling anisotropic control over the directional thickness of the solid phase. The evolution is iterated for a fixed number of steps and followed by adaptive binarization to recover a clean binary field.

As shown in Fig. 11, both \(s_x\) and \(s_y\) produce two coupled effects on the mechanical properties: increasing the solid volume fraction raises the overall stiffness level, while the directional boundary evolution preferentially strengthens the stiffness along the controlled axis. The two parameters act independently along the horizontal and vertical directions, providing simultaneous control over the global stiffness level and the degree of mechanical anisotropy.

2.3 Property space exploration with combined control↩︎

In the previous sections, we have introduced four different post-training steering mechanisms, including orientation, anisotropy, scale, and directional thickness steering. Each of them can control the generated microstructures without retraining the NCA model. Here, we combined these mechanisms by randomly sampling control parameters \(s_1\), \(s_2\), \(s_x\), and \(s_y\) to explore the property space covered by their joint control, which lays the foundation for the multiscale design presented later.

Note that orientation and scale steering are intentionally excluded from the exploration. Orientation steering is excluded because it would introduce a rotational degree of freedom that couples with the directional effects of \(s_1\), \(s_2\), \(s_x\), and \(s_y\), making it difficult to isolate each parameter’s contribution to the stiffness distribution. Moreover, as shown earlier, such orientational changes essentially amount to rotating the stiffness tensor. Scale steering is excluded because it has negligible effect on the mechanical properties. The generated dataset therefore spans a four-dimensional control space \([s_1, s_2, s_x, s_y]\), with each sample consisting of the control parameters, the generated microstructure, and its stiffness matrix \(\mathbf{C}\).

Fig. 12 shows the pairwise and marginal distributions of the six unique stiffness components, \(C_{11}\), \(C_{12}\), \(C_{16}\), \(C_{22}\), \(C_{26}\), and \(C_{66}\), across the property space that spans a broad range with different levels of anisotropy. We use the ratio \(C_{11}/C_{22}\) to quantify anisotropy, where values above one indicate a horizontally dominant response and values below one a vertically dominant one. This ratio ranges from \(0.027\) to \(19.88\), a variation of nearly three orders of magnitude, further quantitatively validating the broad range of anisotropy achieved by the proposed control.

The contribution of each individual parameter to the mechanical property space is illustrated in Fig. 13. The stretching parameter \(s_1\) manipulates the anisotropy determined by \(C_{11}\) and \(C_{22}\), controlling whether the stiffness is horizontally or vertically dominant, while the shear parameter \(s_2\) independently governs the shear stiffness \(C_{66}\), extending the reachable property space toward shear-dominated stiffness that \(s_1\) alone cannot access. The thickness parameters \(s_x\) and \(s_y\) produce complementary directional effects on the overall stiffness level, with \(s_x\) preferentially raising \(C_{11}\) and \(s_y\) preferentially raising \(C_{22}\). Together, the four parameters provide complementary and continuous access to a broad property space, as confirmed by the smooth color gradients observed across all parameter sweeps in Fig. 13. This organized property space forms the foundation for the multiscale design problem presented in Section 4.

3 Spatially varying growth and discretization transfer↩︎

In this section, we extend the generative framework beyond homogeneous microstructures and regular grids. We first show that, by using spatially varying control parameter fields, trained NCAs can grow heterogeneous microstructures with smooth transitions between morphologies and types. We then transfer the trained model from the regular grid on which it was trained to the complex 3D surfaces with irregular mesh by adapting the local perception, without retraining the model. Finally, we demonstrate the programmable growth process by user-defined sources and sequence, enabling growth in irregular geometry and in-situ repair of damaged microstructures.

3.1 Spatially varying control and microstructure transition↩︎

As demonstrated in Sec. 2.2, uniform control parameters enable controlled growth of homogeneous microstructures with statistically equivalent morphology throughout the domain. Many applications, however, require heterogeneous microstructures with spatially varying morphological features to accommodate local requirements. To enable such designs, we extend the steering strategy by prescribing the control parameters locally, allowing them to vary spatially across the domain. This capability is a direct consequence of the local and self-organizing nature of the NCA. When the control parameters are prescribed as smoothly varying spatial fields, each cell follows its local steering conditions while remaining compatible with its neighbors through local information exchange. As a result, the generated morphology transitions smoothly across the domain without introducing discontinuities or requiring any post-processing. In Fig. 14 (a), we first demonstrate such spatially varying growth for the individual control mechanisms. In each case, the corresponding control is prescribed as a smoothly varying spatial field from top to bottom. The same NCA then grows the microstructure under these local steering conditions, producing continuous morphological variations across the domain.

Beyond spatially varying steering within one microstructure type, the framework can also generate transitions between distinct microstructure types. The ring-shaped example in Fig. 14 (b) demonstrates smooth transitions among biofilm, leaf, bone, and spider silk microstructures. We achieve this by introducing a spatial interpolation mask during growth to blend the local update responses of multiple NCA models trained on different parent microstructure types. Because this blending is integrated into the iterative growth process rather than a direct interpolation of final generated structures, as in most existing methods, the smooth transition regions emerge naturally by design and inherit morphological features from both parent microstructures. In the transition region between leaf and bone, for example, the generated morphology combines leaf-like network junctions with bone-like porous features, forming a hybrid microstructure prototype that does not exist in either parent class. This could enable new design possibilities beyond naturally occurring or individually learned microstructures.

To stress test this transition capability, we further consider an extreme case with irregularly mixed microstructure types, as shown in Fig. 14 (c) and Vid. [vid:type transition]. The result shows that the model remains stable under such a complex spatial layout with a smooth transition. Together, these demonstrations extend the framework from global microstructure generation to spatially programmed heterogeneous design.

3.2 Transfer across geometric domains and discretizations↩︎

All previous cases assume growth on a regular pixelated grid. We extend the framework to grow on irregular geometric domains with irregular discretization as shown in Fig. 15, broadening its applicability across design cases. For example, implant design often requires the design to adapt to irregular, patient-specific bone shapes [63]. The key idea of such an extension builds on the PDE interpretation of the NCA dynamics we introduced in Sec. 2.1.1. We note that the same PDE can be solved and converge to the same result regardless of the domain or discretization, as long as the local operators are adapted. For example, a PDE yields the same solution whether solved on a triangular mesh, a hexagonal mesh, or even a point cloud, provided the local operators are redefined in the corresponding form that makes the solution compatible across domains. The governing PDE itself does not change, and only the local operators adapt to different instantiations. We generalize this idea to the NCA. NCA uses a neural network to map a set of local perception operators to the update. The neural network can thus be regarded as the universal PDE form. And to adapt the NCA to a new domain or discretization, we simply adapt the local perception operators, without retraining of the neural network. We take the irregular 3D surface mesh as an example to illustrate the transfer procedure.

Transferring the trained NCA to an irregular surface mesh faces two challenges that do not arise on a regular grid. First, local perception operators such as the \(x\)-gradient are directional, which requires a well-defined coordinate system to provide direction information. Unlike a consistent global coordinate system that generalizes for the whole grid, there is no such global system on the curved surface mesh because of the lack of global axis alignment. To address this issue, a consistent local frame must be established at every vertex so that the directions carry the same meaning across the surface. Second, the fixed weights of the grid kernels are no longer valid on an irregular mesh. Each vertex has a different number of neighbors at different relative positions, so the weights of local perception operators must be redefined intrinsically from the local mesh geometry. We address these challenges following well-established differential geometry tools to redefine local perception operators in local frames at each vertex [64]–[66].

To establish the local frame, we prescribe a global flow direction \(\mathbf{d} \in \mathbb{R}^3\) that defines a reference orientation across the surface [64], [65], as shown in Fig. 15 (a). At each vertex \(i\) with position \(\mathbf{v}_i\) and surface normal \(\mathbf{n}_i\), we project \(\mathbf{d}\) onto the local tangent plane to obtain the local \(x\)-axis: \[\mathbf{f}_i = \frac{(\mathbf{n}_i \times \mathbf{d}) \times \mathbf{n}_i}{\|(\mathbf{n}_i \times \mathbf{d}) \times \mathbf{n}_i\|}, \label{eq:tangent95flow}\tag{10}\] and define the local \(y\)-axis as \(\mathbf{t}_i = \mathbf{n}_i \times \mathbf{f}_i\). Together, \((\mathbf{f}_i, \mathbf{t}_i, \mathbf{n}_i)\) form an orthonormal local frame at the vertex \(i\), providing fixed local coordinate systems for directional perception operators across the surface.

With the local coordinate system of each vertex, we define an intrinsic, angularly-weighted stencil that generalizes the applied kernels to the local tangent frame [66]. Each neighbor \(j\) of the vertex \(i\), with position \(\mathbf{v}_j\), is projected into this local tangent frame to obtain its local 2D coordinates \((u_{ij}, v_{ij})\) relative to the vertex \(i\): \[u_{ij} = (\mathbf{v}_j - \mathbf{v}_i) \cdot \mathbf{f}_i, \qquad v_{ij} = (\mathbf{v}_j - \mathbf{v}_i) \cdot \mathbf{t}_i, \label{eq:neighbor95projection}\tag{11}\] where \(u_{ij}\) and \(v_{ij}\) are the local \(x\)- and \(y\)-offsets of neighbor \(j\) from the center vertex \(i\). On the uniform grid, the Sobel and Laplacian filters assign fixed weights according to the axial and diagonal positions of neighboring pixels. On the mesh, we replace these fixed weights with continuous weights computed from the projected neighbor offsets: \[w^{x}_{ij} = \frac{2\,u_{ij}}{|u_{ij}| + |v_{ij}|}, \qquad w^{y}_{ij} = \frac{2\,v_{ij}}{|u_{ij}| + |v_{ij}|}, \qquad w^{\mathrm{lap}}_{ij} = \frac{2}{|u_{ij}| + |v_{ij}|}, \label{eq:mesh95weights}\tag{12}\] and assign the center vertex the negative sum of its neighbor weights, ensuring that the operators vanish on constant fields. These intrinsic weights produce mesh analogs of the grid gradient and Laplacian kernels, and degenerate exactly to the 2D grid kernels, making the operators compatible with both discretizations. The resulting perception tensor is arranged in the same channel order as the grid perception, yielding a tensor of shape \([B, 4C, V]\) that is passed directly to the trained update network.

As shown in Fig. 15 (b), we demonstrate the transferred growth on a complex 3D surface. Starting from a random initial state, the target morphology begins to emerge coherently and finally reaches a stable state. This result shows that the learned local update rule is not restricted to the uniform grid discretization on which it was trained, provided that local perception supplies compatible local differential information on the target discretization.

3.3 Programmable asynchronous growth↩︎

Beyond controlling the morphology itself, we introduce an additional level of control to NCA by prescribing a growth sequence to control where and when growth occurs. Within the activated region, the NCA dynamics evolve as usual, while cells outside remain frozen. When those frozen regions are activated later, the self-organizing dynamics allow them to interact with the previously activated regions, merging multiple evolutions into one to ensure compatible connections and achieve asynchronous growth. This mirrors developmental processes in nature, where tissues grow at different times and rates yet still merge seamlessly, much as a healing wound regenerates and reintegrates with the surrounding tissue. We demonstrate this capability through two representative examples.

3.3.1 Two-source growth in complex geometry↩︎

The first example demonstrates asynchronous growth inside an irregular domain, as shown in Fig. 16. We represent the domain by a binary mask \(\mathbf{M} \in \{0,1\}\), where cells outside the mask remain fixed at zero throughout growth. We then prescribe two source points, \(\mathbf{m}_1\) and \(\mathbf{m}_2\), inside the domain to define where growth begins and how the active region expands. The growth sequence is determined by a normalized nearest-source distance field: \[Dist(\mathbf{p}) = \frac{\min_k \|\mathbf{p} - \mathbf{m}_k\|_2}{\max_{\mathbf{p}' \in \mathbf{M}} \min_k \|\mathbf{p}' - \mathbf{m}_k\|_2}, \label{eq:distance95field}\tag{13}\] where cells closer to either source have lower values of \(Dist(\mathbf{p})\) and start to grow earlier. The sequence allows the growth to expand outward from both sources simultaneously, as shown in Fig. 16 (a). The source(s) and growth sequence can always be customized for each specific case as required, serving as a flexible interface for programmable growth.

As shown in Fig. 16, the two growth fronts expand from the prescribed sources, merge where their influence regions meet, and eventually fill the irregular domain. We illustrate the whole growth process in Vid. [vid:asynchronous growth]. The result remains coherent and consistent while following the prescribed growth sequence. This example shows that customized growth sources and sequence can program the growth spatially and temporally in arbitrary domain shapes. The asynchronous growth allows parallel generation starting from multiple sources in large-scale structures, without retraining or modifying the learned NCA update network. By enabling only the active portion to grow, i.e., the growth front of each source, it reduces peak memory consumption and enables straightforward parallelization. In contrast to generating the entire domain all at once, which may not be scalable, this provides an alternative path that helps the generation scale to large and complex domains.

3.3.2 In-situ microstructure repair↩︎

Beyond controlling the sequence in irregular domains, asynchronous growth can also be used to regrow a prescribed region surrounded by preexisting microstructures, enabling the in-situ repair of a missing region in a biological microstructure, as shown in Fig. 17. This is achieved by restricting growth to the missing region while keeping the surrounding intact area fixed to serve as the boundary condition of growth. The NCA is thus directed to grow only within the missing region, respecting the existing material boundaries. We demonstrate this on a cross-sectional slice of trabecular bone, from which a circular region is removed to simulate a localized structural defect. Rather than requiring a reference structure in the missing region, we extract a small adjacent intact patch as the single training template for NCA used for the repair.

During repair, the NCA applies the stochastic asynchronous update in the missing region while the exterior remains fixed. The fixed region acts as a boundary condition, allowing the generated structure inside the defect to interact with and connect to the existing trabecular network at the repair interface. We illustrate the two-stage process from training to repairing in Vid. [vid:bone repair]. As a result, the generated microstructure inherits the local morphology of the neighboring template and embeds smoothly with the surrounding bone. This example highlights two key advantages of the framework. First, the one-shot training feature enables the model to learn from a single small neighboring patch, without requiring large patches or any reference structure in the missing region. Second, asynchronous growth allows the new material to grow into the defect and connect to the existing trabecular network rather than forming a separate, mismatched patch. Together, these properties give the framework strong potential for patient-specific implant design, and can serve as a robust, self-restoring representation for other engineering systems where resilience against damage or attack is critical.

4 Multiscale optimization↩︎

In this section, we integrate the proposed NCA generative framework into a multiscale optimization pipeline. We demonstrate this pipeline through a mechanical cloaking design, which features irregular domain geometries and complex heterogeneous property requirements. By following the growth cues provided by the optimization process, the NCA model steers its self-organizing dynamics to generate a heterogeneous multiscale structure with smooth transitions, accommodating complex geometries, and growing naturally across different resolutions without requiring a redesign. This capability enables high-performance cloaking designs without the complicated post-processing or scalability limitations of existing methods, demonstrating the generality and scalability of the proposed framework.

4.1 Optimization framework and problem setup↩︎

4.1.1 General optimization problem↩︎

To translate the controllable microstructure generator into a structural design tool, we couple the NCA model with multiscale optimization. The goal is to generate spatially varying microstructures to achieve a target structural response or maximize the system-level performance. We discretize the macrostructure into a finite-element mesh of \(n_e\) four-node quadrilateral elements. Each design element \(e\) is assigned an effective stiffness matrix \(\mathbf{C}_e\), constrained within the explored property space of a trained NCA in Sec. 2.3 enabled by the proposed hybrid control strategy. The macroscale response is governed by \[\mathbf{K}(\mathbf{C}_e)\,\mathbf{U} = \mathbf{F}, \label{eq:fe95equilibrium}\tag{14}\] where \(\mathbf{K}\) is the global stiffness matrix assembled from the constitutive matrix in each element, \(\mathbf{U}\) is the nodal displacement vector, and \(\mathbf{F}\) is the applied load vector. The general multiscale optimization problem is defined as: \[\min_{\mathbf{C}_e} \quad f\!\left(\mathbf{U}(\mathbf{C}_e)\right) \label{eq:general95opt}\tag{15}\] \[\text{s.t.} \quad \mathbf{K}(\mathbf{C}_e)\,\mathbf{U} = \mathbf{F},\] \[\mathbf{C}_e \in \mathcal{M}, \quad e \in \Omega_d,\] where \(\Omega_d\) is the set of design elements, \(f(\cdot)\) is a given objective, and \(\mathcal{M}\) denotes the explored property space. The elemental stiffness matrix \(\mathbf{C}_e\) in the 2D setting is fully characterized by six independent stiffness components in Voigt notation \([C_{11}, C_{12}, C_{22}, C_{16}, C_{26}, C_{66}]\).

For a given microstructure template, we can obtain the reachable property space of the trained NCA by varying the control parameters, as shown in Fig. 12. A straightforward approach is to directly optimize the stiffness of each element within this explored property space. However, while feasible, this strategy is not ideal and can be challenging in practice. The reason is that the feasible property space is highly skewed, as shown in Fig. 12, resulting in a high-dimensional irregular boundary. This makes it difficult to formulate effective constraints on the stiffness components. Even with such constraints, the optimizer may frequently violate the feasible region or become trapped during optimization.

To address this issue, we introduce a variational autoencoder (VAE) [44], [45] as a dimensionality reduction and regularization tool to construct a bidirectional mapping between the stiffness components and a lower-dimensional latent space \(\mathbf{z}_e \in \mathbb{R}^{d_z}\). The neural network encoder maps the stiffness matrix \(\mathbf{C}_e\) to the latent representation \(\mathbf{z}_e\), while the paired decoder reconstructs the stiffness matrix from the latent variables, i.e. \(\mathbf{C}_e = D(\mathbf{z}_e)\). Instead of optimizing and constraining the original stiffness components \(\mathbf{C}_e\), we can now reformulate the optimization problem in the latent space, which is lower dimensional and exhibits a more regular distribution (ideally Gaussian). This transformation greatly simplifies constraint construction and improves optimization efficiency. Specifically, the optimization is reformulated in terms of the latent variables \(\mathbf{z} = \{\mathbf{z}_e\}_{e \in \Omega_d}\): \[\min_{\mathbf{z}} \quad f\!\left(\mathbf{U}(\mathbf{z})\right) \label{eq:latent95opt}\tag{16}\] \[\text{s.t.} \quad \mathbf{K}(\mathbf{z})\,\mathbf{U} = \mathbf{F},\] \[\mathbf{g}(z)\leq0,\]

where \(\mathbf{g}(\mathbf{z})\) is a feasibility constraint in the latent space to ensure that the optimized design remains within the explored property space. Specifically, the feasible properties obtained in Sec. 2.3 by varying the control parameters serve as anchor points that span the explored property space. We require the optimized latent variables to remain close to these anchors and formulate \(\mathbf{g}(\mathbf{z})\) as a latent tube constraint. For each design element, we compute the normalized distance \(\mathrm{dist}_e\) to the nearest feasible latent anchor and define the tube radius as \(r_0 = \sqrt{d_z} \cdot r_{\mathrm{frac}}\), where \(r_{\mathrm{frac}}\) is a given fractional value to decide the radius of feasible range. A smooth activation function is then used to penalize latent variables outside the tube: \[h_e = \frac{1}{2}\left(\tanh\!\left(\beta \left(\mathrm{dist}_e - r_0\right)\right) + 1\right), \label{eq:tube95activation}\tag{17}\] where \(\beta\) increases progressively during optimization to sharpen the constraint boundary. The activation \(h_e\) varies smoothly between 0 and 1: it approaches 0 when \(\mathrm{dist}_e < r_0\), indicating that the latent variable lies safely inside the tube, and approaches 1 when \(\mathrm{dist}_e > r_0\), indicating that it has moved outside the feasible neighborhood. Therefore, an ideal design requires \(h_e\) to remain close to 0 for all design elements. Aggregating over the domain, the feasibility constraint is expressed as \[g(\mathbf{z}) = \frac{1}{|\Omega_d|}\sum_{e \in \Omega_d} h_e - \gamma \leq 0, \label{eq:tube95constraint}\tag{18}\] where \(\gamma\) is the allowable fraction of elements outside the tube. In addition to this smooth constraint, as a further measure to ensure feasible designs, any element with latent variables that falls outside the tube after each optimization iteration is projected to its nearest feasible anchor points. We emphasize that the VAE is applied to the explored property space rather than the microstructure geometry. Its purpose is solely to facilitate the formulation of constraints in multiscale optimization. The microstructure geometry is learned by the NCA framework and still requires only a single template as the training data.

The sensitivity of the objective with respect to the latent variables is obtained by chaining the adjoint sensitivity of \(f\) with respect to the stiffness components through the Jacobian of the decoder: \[\frac{\partial f}{\partial \mathbf{z}_e} = \mathbf{J}_e^{\top} \frac{\partial f}{\partial \mathbf{C}_e}, \qquad \mathbf{J}_e = \frac{\partial D(\mathbf{z}_e)}{\partial \mathbf{z}_e}, \label{eq:chain95rule}\tag{19}\] where \(\mathbf{J}_e \in \mathbb{R}^{6 \times d_z}\) is the Jacobian of the decoder at \(\mathbf{z}_e\), obtained by automatic differentiation, and \(\partial f / \partial \mathbf{C}_e\) is the adjoint sensitivity. The sensitivity of the aggregated constraint can be obtained in a similar way. The sensitivities of the objective function and constraints enable the use of gradient-based optimizers to navigate the latent space iteratively, such as the Method of Moving Asymptotes (MMA) adopted in this study.

4.1.2 Case study: mechanical cloak design↩︎

While the proposed optimization is general and can accommodate different multiscale design requirements, we choose to demonstrate it on a mechanical cloak design [67], due to its demanding requirement for spatially varying microstructures with anisotropic mechanical properties in different regions, making it an ideal stress test of the framework’s ability to translate optimized property fields into a heterogeneous structure with a smooth transition.

In this case, the macrostructure is a \(50 \times 50\) finite-element grid divided into three regions. As shown in Fig. 18 (a), we design the mechanical cloak region \(\Omega_c\) around a void \(\Omega_v\) within the structure. The goal is to make the displacement field \(\mathbf{U}_{\mathrm{ref}}\) in the surrounding reference region \(\Omega_r\) retain the values it would have in a homogeneous structure without the void. The void region is fixed to a near-zero stiffness, and the reference region is fixed to the stiffness of the reference microstructure \(\mathbf{C}^{\mathrm{ref}}\). The properties of the cloaking region \(\Omega_c\) are optimized to minimize the void-induced distortion of the displacement field. We define the objective as: \[f(\mathbf{U}) = \frac{\left\|\mathbf{\gamma} \odot \left(\mathbf{U} - \mathbf{U}_{\mathrm{ref}}\right)\right\|_2^2}{\left\|\mathbf{\gamma} \odot \mathbf{U}_{\mathrm{ref}}\right\|_2^2 + \varepsilon}, \label{eq:cloaking95obj}\tag{20}\] where \(\mathbf{\gamma}\) is a binary mask to distinguish the reference region and \(\varepsilon\) is a regularization constant. We apply constant compressive pressure on the left and right edges to induce deformation, while rigid-body motion is suppressed, as shown in Fig. 18 (a). All design elements are initialized with the latent vector of reference property \(\mathbf{C}^{\mathrm{ref}}\). The sensitivity of the objective with respect to the stiffness of each design element is computed using the adjoint method, which avoids solving the forward sensitivity equation separately for every design variable. The adjoint variable \(\boldsymbol{\lambda}\) is obtained from the global system \(\mathbf{K}^{\top}\boldsymbol{\lambda} = \partial f / \partial \mathbf{U}\). Substituting this adjoint solution and the explicit form of \(\partial f/\partial \mathbf{U}\) derived from the objective in Eq. 20 , the full sensitivity with respect to stiffness component \(C_{e,ij}\) of design element \(e\) is: \[\frac{\partial f}{\partial C_{e,ij}} = -\left[ \frac{2\,\boldsymbol{\gamma}\odot(\mathbf{U}-\mathbf{U}_{\mathrm{ref}})}{\left\|\boldsymbol{\gamma} \odot \mathbf{U}_{\mathrm{ref}}\right\|_2^2 + \varepsilon} \right]^{\top} \mathbf{K}^{-1} \frac{\partial \mathbf{K}}{\partial C_{e,ij}} \mathbf{U}. \label{eq:sensitivity}\tag{21}\] The complete optimization loop is summarized in Fig. 18 (b). In iteration \(t\), the current latent field \(\mathbf{z}^t\) is decoded into an element-wise stiffness matrix \(\mathbf{C}^t\), which is assembled into the global stiffness matrix for the finite element simulation. The resulting displacement field \(\mathbf{U}^t\) is compared with the homogeneous reference field. The adjoint solve provides the sensitivity of the objective with respect to the element stiffness components, and the decoder Jacobian maps this sensitivity to the latent variables. The MMA optimizer then updates the latent field using the objective gradient together with the feasibility constraint. After each update, latent points that leave the feasible region are projected back to their nearest samples in the dataset. The loop continues until convergence, progressively tailoring the stiffness distribution in \(\Omega_c\) to minimize the displacement mismatch in the reference region. Detailed setups are shown in Table 3.

4.2 Optimized results↩︎

The optimized stiffness component distributions and the corresponding NCA control parameter fields within the cloaking region are shown in Fig. 19 (a). The stiffness components in Fig. 19 (a) indicate a directional redistribution of the effective stiffness around the void, which manipulates the elastic response to minimize the void-induced distortion of the displacement field and recover the homogeneous displacement response in the reference region. Fig. 19 (b) shows the NCA control parameter fields after projection onto the explored property space, where \(s_1\) and \(s_2\) control spatially varying anisotropy and \(s_x\) and \(s_y\) control directional thickness at each element. The projected control parameter fields will guide the growth of microstructures to realize the optimized stiffness distribution. In the homogenized simulation, it reduces the relative displacement error to \(8.10\%\) compared with \(129.41\%\) for the uncloaked structure, as reported in Table 1.

Using these spatially varying control parameter fields, we then generate full-size cloaking structures with the trained NCA model. The control parameters act as growth cues that steer the self-organizing dynamics of the trained NCA, analogous to how local environmental cues guide natural growth. Because optimization is performed on properties rather than on any specific microstructure, the optimized design can be realized at any target resolution and on arbitrary discretizations, from regular grids to irregular meshes, without assembling discrete building blocks or retraining. As a demonstrative case, Figs. 20 and 21 show the generated structures and their displacement fields at \(500 \times 500\) and \(1000 \times 1000\) resolutions, respectively. In each case, we compare three configurations: the homogeneous reference structure, the uncloaked but voided structure, and the cloaked structure generated by the NCA. The void in the uncloaked structure strongly distorts both \(U_x\) and \(U_y\), leading to the largest relative error in Table 1. In contrast, the NCA-generated cloaked structure substantially restores the smooth displacement distribution of the homogeneous reference in the reference region. It can be noted that there are smooth morphological transitions, which are inherently guaranteed by the local self-organizing dynamics of the NCA, without requiring any post-processing steps. This provides a major advantage over existing building block assembly methods [44], which suffer from scalability and compatibility issues and cannot easily adapt to different discretizations or meshes.

The resolution study further demonstrates the flexibility of the generation framework. The \(1000 \times 1000\) structure produces finer microstructural detail and smoother displacement fields than the \(500 \times 500\) structure, with the mean relative displacement error over 20 independent samples decreasing from \(26.12\% \pm 1.71\%\) to \(16.86\% \pm 0.63\%\), indicating that higher resolution improves both the average cloaking performance and the consistency across independently generated samples. This improvement reflects the statistical nature of the disordered microstructures: the optimized effective properties represent the average response of a representative volume, which is perfectly recovered only when each design region contains enough features. Finer microstructures place more features within each region, reducing the property fluctuations and thereby improving both the accuracy and the consistency of the generated structure. Even at \(500 \times 500\) resolution, the NCA-generated structure still significantly reduces the displacement mismatch relative to the uncloaked case. These results confirm that optimized control parameter fields can be freely deployed across different scales, offering the flexibility to adapt the target resolution to practical manufacturing constraints.

5 Conclusion↩︎

In this work, we developed an NCA-based generative framework for disordered microstructures. Instead of directly learning the geometry of microstructures, the framework learns the underlying local rules that grow them. This makes it highly data-efficient, requiring only a single template to enable microstructure growth across different domains and scales. The locality of the NCA update rule and its self-organizing nature ensure good connectivity and smooth transitions in spatially varying microstructures by design. This capability mirrors natural growth processes, where spatially varying cues enable heterogeneous growth while maintaining seamless integration across regions. A byproduct of this capability is that the NCA framework enables controlled asynchronous growth, allowing flexible control over where and when microstructures evolve. Similar to biological tissues, different regions can now grow asynchronously while maintaining seamless connectivity. This capability facilitates efficient parallel generation of microstructures over large and irregular domains and opens opportunities for adaptive and self-repairing material systems. These features distinguish the framework from existing building-block assembly methods, which require interface compatibility treatments and often suffer from scalability limitations.

An important insight from this study is that the learned NCA can be viewed as a generalized PDE system, allowing mature computational tools and analytical principles from the PDE domain to be leveraged for microstructure generation. Beyond simply replicating the template morphology, we enable continuous and interpretable control over the orientation, anisotropy, scale, and directional thickness of the generated microstructure by modifying the local perception operator, without additional training. This allows a single trained model to cover a broad morphological and mechanical property space that extends well beyond the training template, a capability that conventional deep generative models cannot easily achieve. Furthermore, by adapting only the local perception operator based on established principles from PDE discretization, the learned growth rule can be transferred across different domain geometries and representations, including planar and curved surfaces as well as regular grids and irregular meshes. This representation independence provides substantially greater flexibility than existing generative models and assembly-based approaches that are tied to predefined geometries and discretizations.

With these new capabilities enabled by the NCA, we extend beyond microstructure generation and integrate the trained NCA model into multiscale structural design. Demonstrated through a mechanical cloaking problem, this approach establishes an extensible pathway from a single microstructure template to spatially programmed heterogeneous structures, highlighting the potential of local self-organizing generative principles for next-generation material design.

The proposed framework lays the foundation for several promising future extensions. First, while current controls steer the morphology and thereby influence the mechanical properties, conditioning target mechanical properties directly into the NCA model would extend the current framework to on-demand inverse generation of microstructures with target properties. Second, the framework can be extended beyond linear elastic properties to consider nonlinear mechanical responses, which will broaden the applicability of the framework to a wide range of mechanical behaviors that motivate disordered microstructure design. Third, the current 2D NCA can be extended to generate 3D microstructures by integrating 3D perception operators, such as 3D convolution kernels on voxel grids or point clouds.

We envision that the proposed framework will serve as a general generative engine for mechanical metamaterial design. Its one-shot efficiency, controlling flexibility and geometry adaptivity can help design in engineering domains that demand disordered microstructures with prescribed mechanical properties, including biomedical implant design [68]–[71] and soft robotics [72]–[74].

6 Model setups↩︎

Table. 2 provides the definition of hyperparameters in detail and the choice to train the NCA model. Table. 3 provides the detailed setup for the optimization process.

7 Supplementary videos↩︎

[]{#vid:nca growth label=“vid:nca growth”}Video S: The growth process of leaf vein microstructure in Fig. 2.

[]{#vid:individual control label=“vid:individual control”}Video S: Growth process under individual controls in Fig. 6, 8, 9, 11.

[]{#vid:type transition label=“vid:type transition”}Video S: Growth transitioning across multiple types of microstructure with irregular boundaries in Fig. 14 (c).

[]{#vid:asynchronous growth label=“vid:asynchronous growth”}Video S: Asynchronous growth in an irregular domain with customized growth sources and sequence in Fig. 16.

[]{#vid:bone repair label=“vid:bone repair”}Video S: In-situ bone repair process in Fig. 17.

Acknowledgments↩︎

The authors acknowledge the support from the Department of Mechanical Engineering at Carnegie Mellon University.

References↩︎